Quantum Simulation with Interacting Photons

Enhancing optical nonlinearities so that they become appreciable on the single photon level and lead to nonclassical light fields has been a central objective in quantum optics for many years. After this has been achieved in individual micro-cavities representing an effectively zero-dimensional volume, this line of research has shifted its focus towards engineering devices where such strong optical nonlinearities simultaneously occur in extended volumes of multiple nodes of a network. Recent technological progress in several experimental platforms now opens the possibility to employ the systems of strongly interacting photons these give rise to as quantum simulators. Here we review the recent development and current status of this research direction for theory and experiment. Addressing both, optical photons interacting with atoms and microwave photons in networks of superconducting circuits, we focus on analogue quantum simulations in scenarios where effective photon-photon interactions exceed dissipative processes in the considered platforms.

Abstract. Enhancing optical nonlinearities so that they become appreciable on the single photon level and lead to nonclassical light fields has been a central objective in quantum optics for many years. After this has been achieved in individual micro-cavities representing an effectively zero-dimensional volume, this line of research has shifted its focus towards engineering devices where such strong optical nonlinearities simultaneously occur in extended volumes of multiple nodes of a network. Recent technological progress in several experimental platforms now opens the possibility to employ the systems of strongly interacting photons these give rise to as quantum simulators. Here we review the recent development and current status of this research direction for theory and experiment. Addressing both, optical photons interacting with atoms and microwave photons in networks of superconducting circuits, we focus on analogue quantum simulations in scenarios where effective photon-photon interactions exceed dissipative processes in the considered platforms. Quantum simulation [36,50,89] is a useful concept in several situations. One application considers the emulation of physical phenomena that are in their original form not accessible with existing experimental technology. This could for example be because the required energy or length scales exceed anything that is realisable.
Another, possibly more prominent application of quantum simulations, which is the subject of this review is quantum many-body physics.
The dimension of the Hilbert space for a quantum many-body system grows exponentially in the number of its constituents (the number of particles). Hence, to specify an arbitrary quantum state of such a system one needs to store an exponentially growing amount of complex numbers, the probability amplitudes of the state's components. This task quickly becomes impossible on available classical computers. With numerical simulations of such systems being out of reach, one could hope to nonetheless explore them in detail via experiments. Yet, unfortunately it is in many situations not possible to resolve their microscopic properties.
Here quantum simulation promises a way forward. Complex quantum many-body systems that defy experimental access are emulated by other systems which allow for much better experimental control and measurement resolution. The latter is often due to different temperature, energy, length or time scales at which quantum simulators work as compared to their simulation targets. This idea to simulate a complex quantum system with another well controllable one was originally proposed by Feynman [82] by suggesting that a 'computer' built of quantum mechanical elements is needed to simulate highly complex quantum systems.
The observation that an exponentially growing amount of data is needed to specify a quantum state on the other hand indicates that quantum dynamics can execute computations with massive parallelism. Indeed there are known examples where a quantum computer could solve a problem exponentially faster than any known classical algorithm [74]. In the mid 1990s, Lloyd showed that a universal quantum computer would also be able to simulate any quantum system efficiently [169].
With its conceptual use being well understood, the concept of quantum simulation received a boost in terms of practical implementations through the successful emulations of Bose-Hubbard models in equilibrium with ultra-cold atoms in optical lattices [94]. Quantum simulations with ultra-cold atoms have since then developed into a very active and successful field of research, see [29] and [28] for recent reviews. Concepts for quantum simulation have also been successfully pursued with trapped ions [27] and in quantum chemistry [134,172].
Quantum simulators can operate in two conceptually different ways. In digital quantum simulation, the time evolution to be simulated is decomposed into a sequence of short evolution steps using the Trotter formula. The individual steps of such a decomposition are often very similar to quantum gates, so that a digital quantum simulation can be viewed as a special purpose quantum computation. In analogue quantum simulation, in contrast, the device is operated in such a way that it's dynamics can, for suitable initial conditions, be approximated by an effective timeindependent Hamiltonian, the physics of which one aims to simulate.
In this review we will summarise the developments towards quantum simulations with photons throughout recent years. In doing so, we will concentrate on analogue simulations of quantum many-body physics with photons that are made to interact with each other via optical nonlinearities. We focus on regimes where the effective photon-photon interactions are a dominant process and in particular exceed the dissipation processes for individual photons.
We will thus only very briefly touch developments in digital quantum simulations but not cover the progress in simulations with linear optics devices and refer interested readers to the recent review by Aspuru-Guzik and Walter [11].
Generating strong optical nonlinearities or effective photon-photon interactions has been a central goal of research in quantum optics for many years. Since more than a decade, optical nonlinearities on the single photon level have now been realised in single cavities that confine the trapped light fields as strongly as possible [254]. Due to technological progress in the micro-fabrication of high finesse resonators it has now become feasible to couple several resonators coherently to form an extended network [141]. Alternatively onedimensional waveguides doped with optically highly nonlinear media can be considered. These developments give rise to networks or extended volumes where strong light-matter coupling and hence appreciable effective photon-photon interactions take place simultaneously in multiple locations.
The term "interacting photons" used in the title calls for some explanation. Pure photons do only exist in infinitely extended vacuum. Any matter that light enters into or any boundary conditions, e.g. in the form of a waveguide, it is subjected to will cause the emergence of a polarisation in these media. The elementary excitations of such fields are then polaritons, combinations of photons and excitations of the polarisation-fields, and these do interact if the medium has an appreciable non-linearity. It is photons in this sense that are considered in this review.
The remainder of this review is organised as follows. In section 2, we review the development towards quantum simulation with optical photons that are subject to strong effective interactions mediated by strongly polarisable atomic media or quantum well structures. We cover the research on cavity arrays leading to effective Bose-Hubbard or spin models and Jaynes-Cummings-Hubbard models, see section 2.1, as well as work on continuum models leading to Lieb-Lininger models and generalisations thereof, see section 2.2. Finally we briefly touch on approach employing atomic Rydberg media in section 2.3. In section 3, we then describe more recent developments with microwave photons hosted in networks of superconducting circuits. After briefly introducing the quantum theory for electrical circuits, c.f. section 3.1, we cover developments towards Bose-Hubbard, c.f. section 3.2, and Jaynes-Cummings-Hubbard models, c.f. section 3.3, in these devices. We then point out some developments in digital quantum simulation, c.f. section 3.4, and discuss recent experimental progress, c.f. section 3.5. In section 4, we discuss the quantum many-body dynamics that can be explored in the considered quantum simulators. We first review equilibrium calculations of phase diagrams, c.f. section 4.1, then discuss non-equilibrium and driven-dissipative regimes, c.f. section 4.2, and finally comment on employed and recently developed calculation techniques, c.f. section 4.3, and the experimental signatures of the predicted phenomena, c.f. section 4.4. We then discuss the more recent work on quantum simulation of many-body systems subject to artificial gauge fields, see section 5, which has already seen some important experimental advances with photons and conclude with a summary in section 6.

Optical photons interacting with atoms
First theory developments towards quantum simulations with interacting photons considered optical photons that couple to atoms. In order to generate sufficiently strong optical nonlinearities for this aim, strongly confined light-fields with appreciable interactions between individual photons and atoms where considered. Such strong atom-light interactions had been achieved in high finesse optical cavities [26] in the years CONTENTS 4 before and so multiple cavities coupled via tunnelling of photons between them were considered initially, see also the reviews by Hartmann et al. [112], Tomadin et al. [242] and Noh et al. [190]. Yet, for optical frequencies, the strength of achievable light matter couplings is about six orders of magnitude smaller than the frequency of the employed photons and it is thus quite demanding to build multiple cavities with sufficiently low disorder so that photons can tunnel efficiently between them. As a consequence setups with thinly tapered optical fibres have been considered subsequently to avoid these challenges, see also a recent review by Roy et al. [218]. We here first review the work on cavity arrays and then turn to discuss continuum models in tapered optical fibres.
2.1. Lattice models in cavity arrays 2.1.1. Coupling of cavities 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, ( r) = ( r + R n), where r is the three dimensional position vector, R the lattice constant (i.e. the distance between the centres of adjacent cavities) and n a vector of integers. Expanding the electromagnetic field in terms of Wannier functions w R , localised in the cavities at locations R = R n, the tunnelling rate of photons between neighbouring cavities can be expressed as where | R − R | = R, ω C is the resonance frequency of the considered cavity mode and the dielectric function R ( r) describes a single cavity surrounded by bulk material only. Introducing creation and annihilation operators a † R and a R for the Wannier modes, the Hamiltonian of the field can thus be written as where < R, R > is a sum over nearest neighbours. Since the tunnelling rate is typically much less than the photon frequency, J phot ω 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 cavitycavity interactions. In practise there will always be some disorder in the array and the resonance frequencies, ω C ( R), and tunnelling rates, J phot ( R, R ) 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.

Bose-Hubbard Models
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, where p † R creates a boson at site 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 > 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 [94,29,28].
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.

Bose-Hubbard model with dark state polaritons
A Bose-Hubbard model can be generated in an array of cavities that are filled with atoms of a specific An array of cavities as described by our model. Photon hopping occurs due to the overlap (shaded green) of the light modes (green lines) of adjacent cavities. Atoms in each cavity (brown), which are driven by external lasers (blue) give rise to an on site potential. Reprinted with permission from [21].
Although such an interaction naturally appears in some media as a result of a non-zero third order electric susceptibility, its effect is negligible on the level of individual quanta, which is a reason for the difficulty of realising nonlinear optics for individual quantum systems. Using the enhanced light-matter interaction in cavity QED, it is possible to engineer much stronger nonlinearities. Intuitively the strong interactions of the light mode with atoms inside the cavity, under particular circumstances, mediates strong nonlinear interaction among the photons of the cavity mode. The strength of the nonlinearity can be increased even further if instead of considering photons as the bosonic particles of the model, one considers polaritons, joint photonicatomic excitations. In the next section we review the proposal for realising the Bose-Hubbard Hamiltonian using such polaritons.

Polaritonic Bose-Hubbard model
The polaritonic Bose-Hubbard model is hosted in an array of cavities that are filled with atoms of a particular 4 level structure, which are driven with an external laser, see figures 3 and 4. Thereby the laser drives the atoms in the same manner as in Electromagnetically Induced Transparency (EIT) [76]: The transitions between levels 2 and 3 are coupled to the laser field and the transitions between levels 2-4 and 1-3 couple via dipole moments to the cavity resonance mode. Levels 1 and 2 are assumed to be metastable and their spontaneous emission rates are thus negligible.
It has been shown by Imamȏglu and co-workers, that this atom cavity system can exhibit a very large nonlinearity [77,78,79,80], and a similar nonlinearity has been observed experimentally [67].
Considering the level structure of figure 4, in a rotating frame with respect to H 0 = ω C a † a + 1 2 + N j=1 ω C σ j 22 + ω C σ j 33 + 2ω C σ j 44 , the Hamiltonian of The level structure and the possible transitions of one atom, ωC is the frequency of the cavity mode, Ω is the Rabi frequency of the driving by the laser, g13 and g24 are the parameters of the respective dipole couplings and δ, ∆ and ε are detunings.
the atoms in the cavity reads, where σ j kl = |k j ⟩⟨l j | projects level l of atom j to level k of the same atom, ω C is the frequency of the cavity mode, Ω is the Rabi frequency of the driving by the laser and g 13 and g 24 are the parameters of the dipole coupling of the cavity mode to the respective atomic transitions.
In a cavity array, N atoms in each cavity couple to the cavity mode via the interaction H I and photons tunnel between neighbouring cavities as described in equation (5). Hence the full Hamiltonian describing this system reads and can emulate a Bose-Hubbard model for polaritons as we will see in the following. Assuming that all atoms interact in the same way with the cavity mode, the description can be restricted to Dicke type dressed states, in which the atomic excitations are delocalised among all the atoms. In the case where g 24 = 0 and ε = 0 , level 4 of the atoms decouples from the dressed-state excitation manifolds [79]. The Hamiltonian (8) can then be expressed in terms of the following creation (and annihilation) operators: Copyright line will be provided by the publisher Figure 1. An array of nonlinear cavities. Photon hopping occurs due to the overlap (shaded green) of the light modes (green lines) of adjacent cavities. Atoms in each cavity (brown), which are driven by external lasers (blue) give rise to an on site potential. Reprinted from [108]. Figure 2. The level structure and the possible transitions of one atom, ω C is the frequency of the cavity mode, Ω is the Rabi frequency of the driving by the laser, g 13 and g 24 are the parameters of the respective dipole couplings and δ, ∆ and ε are detunings.
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 and |3 couple to the laser field whereas the transitions |2 ↔ |4 and |1 ↔ |3 couple to the cavity resonance mode. Levels |1 and |2 are assumed to be metastable with negligible spontaneous emission rates.
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,77,207,75,32,152] 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 in a rotating frame with respect to Here σ kl j = |k j l j | is the transition operator between levels |l and |k of atom j, ω C is the frequency of the cavity mode, Ω is the Rabi frequency of the laser drive and g 13 (g 24 ) are the coupling strengths of the cavity mode to the atomic transitions |2 ↔ |4 (|1 ↔ |3 ).
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 lightmatter 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 where g = √ N g 13 is a collective coupling rate, S † 12 = 1 √ N N j=1 σ 21 j creates a spin wave in the metastable atomic levels and B = g 2 + Ω 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 of the atoms induces an effective interaction between all dark state polaritons in one cavity. For |g 24 gΩ/B 2 | |∆| the strength of this interaction can be calculated in a perturbative manner [257], Similarly, the two photon detuning ε leads to an energy shift of µ BH = 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 CONTENTS 6 rate [108] 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 µ = µ BH . Note that both, J BH and U BH can be tuned by varying the intensity of the laser driving Ω 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.

Bose-Hubbard model with exciton polaritons
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, where a † (a) creates (annihilates) a cavity photon and b † (b) creates (annihilates) an exciton. The excitons have transition frequencies ω X , interact with each other at a strength U X and couple to photons at a strength g. An anharmonic excitonphoton 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 and frequencies with sin(θ) = g/ g 2 +∆ 2 , cos(θ) =∆/ g 2 +∆ 2 and Since g U X the dynamics of both polariton species decouples from each other and the interaction between excitons leads to an interaction 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.
2.1.5. Jaynes-Cummings-Hubbard model 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], H JC = ω C a † a + ω 0 |e e| + g(a † |g e| + a |e g|), (11) where ω C and ω 0 are the frequencies of the cavity mode and the atomic transition, g is Jaynes-Cummings lightmatter coupling, a † is the creation operator of a photon in the resonant cavity mode, and |g (|e ) are the ground (excited) states of the two level system. Since the energy of a photon ω C and the atomic transition energy ω 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 = nω C + ∆ 2 ± ng 2 + ∆ 2 4 for n ≥ 1, where ∆ = ω 0 − ω 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 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 [77,207,75] and a circuit quantum electrodynamics (circuit QED) setup [32,152], see also section 3 for further details.

7
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 ∆ = 0 one finds (13) so that the interaction energy per excitation decreases with growing excitation number n, U JC;n /n → g/n for n 1. The degrading of the anharmonicity of the spectrum of the Jaynes-Cummings model at high excitation numbers was observed in a selftrapping experiment [204], see Sec. 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 bimodality 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, where 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 manybody 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 [212,157], see section 4 for further discussions.

Spin Models
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,109,112,92,8]. 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 [109,112,262,92,8].
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.

Experimental requirements for lattice models
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 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 2 13 κγ and g 2 24 κγ, for the Bose-Hubbard model, or a strong coupling regime, g κ, γ for the Jaynes-Cummings-Hubbard model. For the approach to realise a Bose-Hubbard model as explained in section 2.1.2 the condition ∆ > γ 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 κ [167] but requires low rates of pure dephasing γ * [80].

8
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 → |2, 0 , destructively interferes with the excitation path where the first excitation tunnels to the second cavity, |1, 0 → |0, 1 , a second excitation is generated in cavity one via the drive, |0, 1 → |1, 1 , 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 γ * . Despite these requirements for an experimentally suitable technology, there has already been significant progress as we discuss next.
2.1.8. Experimental progress 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 [175,126]. 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 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 U n(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 figure 3a-c.
For interactions as present in the Jaynes-Cummings model, the interaction energy per excitation degrades as the number of excitations grows, see Eq. (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 figure 3d-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.

Continuum Models in Optical Fibers
For optical frequencies, building mutually resonant cavities of sufficient finesse is very challenging. This can be appreciated by observing that the largest achievable atom-photon couplings reach 1 -10 GHz [254]. Hence disorder in the resonance frequencies of the cavities needs to be suppressed to 10 9 Hz or below, which corresponds to a disorder in cavity dimensions below 10 −6 times the wavelength of the trapped photons. A possible alternative to cavity arrays are therefore one-dimensional waveguides, which avoid the need to build mutually resonant cavities but nonetheless feature a large light-matter coupling due to a tight confinement of the light modes in transverse directions. Moreover, in contrast to lattice structures, such devices emulate a different class of quantum many-body models. The probably most prominent |1i |2i Figure 4. One dimensional waveguide, here a hollow core fibre doped with N atoms that have the internal level structure depicted in figure 2. Blue arrows represent classical control fields whereas green arrows represent quantum probe fields.
representative of this class is the Lieb-Liniger model, which describes bosons of effective mass m eff in one dimension which interact via a contact interaction of strengthg ( is Planck's constant). In equation (16), L is the length of the waveguide and ψ the field describing the polaritons. All ground state features of the Lieb-Liniger model [166] are characterised by a single, dimensionless parameter G = m effg /( 2 N p /L) [N p is the particle number], that quantifies the effective interaction strength between the particles. For weak interactions, i.e. |G| < 1, the bosons are in a superfluid state. In contrast, in the strongly correlated regime |G| 1, they form a Tonks-Girardeau gas [91] of impenetrable hard-core particles. A characteristic feature of this regime is that the density-density correlations withn(z) = ψ † (z)ψ(z) vanish for z = z and exhibit Friedel oscillations [87] for |z − z | finite. The Friedel oscillations indicate a crystallisation of the particles by showing that they prefer to occur at specific separations from one another. For further information about the physics of interacting bosons in one dimension we here refer the reader to the review by Cazalilla et al. [44] and now proceed to discussing realisations of this physics with interacting photons.

Lieb-Liniger Model with Dark-State Polaritons
A realisation of the model (16) with dark-state or slowlight 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 Fig 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]. 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 where k c is the momentum of the control fields, sin θ = √ N g 13 / N g 2 13 + 2Ω 2 , cos θ = √ 2Ω c / N g 2 13 + 2Ω 2 and N −1/2 N µ=1 σ 12 µ e −ikzµ describes a spin coherence. The dynamics of the dark state polaritons can be described in terms of their density matrix ρ D . Neglecting single particle dissipation which can be sufficiently suppressed for suitable parameters, ρ D obeys the equation of motion [138,139], where H eff = H LL with effective mass m eff = − (N g 2 13 + 2Ω 2 )/(2δc 2 cos 2 θ) and a complex interaction constantg = 2 Lg 2 24 cos 2 θ/(∆ − cos 2 θ∆ω + iγ 42 /2). Here c is the speed of light. The term 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 and is maximal for purely dissipative (∆ = 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 densitydensity correlations, one finds g (2) for z = z which vanishes in the limit |G| → ∞. 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 ψ 1 and ψ 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 with strength V 12 [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].

Polariton correlations and dynamics
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.

Experimental Progress
There has been substantial progress towards realising quantum manybody systems of strongly interacting photons with optical fibres, although the simulation of a Lieb-Liniger model appears to still pose challenges. Laser-cooled 87 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 ↔ 4, see figure 2, the initial transmission of a beam coupled to the transition 1 ↔ 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.

Polaritons in ensembles of Rydberg atoms
For generating stronger interactions between polaritons, ensembles of Rydberg atoms have more recently received increasing attention. Here the large dipole moment of Rydberg atoms leads to large van der Waals interactions between two atoms that scale proportional to the sixth power of the principal quantum number. These interactions lead to an effect called Rydberg blockade which describes a situation where a driving field cannot generate a second Rydberg excitation in the vicinity of a Rydberg excitation due to these strong forces. Reviewing the development in this branch of research is beyond the scope of this article and we thus refer the interested reader to a series of excellent recent reviews [221,171,115,46,218] and references therein.

Microwave photons in superconducting circuits
In the quest for realising strong effective photonphoton interactions, superconducting circuits supporting elementary excitations in the form of microwave photons have received increasing attention as they offer very favourable properties for this task. Importantly the wavelength of microwave photons are about 10 mm. Therefore resonators that trap them are of a similar size and the accuracy of available micro-fabrication techniques is sufficient for producing multiple resonators of the same resonance frequency on the same chip.
More precisely, any residual disorder in resonance frequencies is well below typical values for their mutual coupling [246]. In this sense, building large, coherently coupled resonator arrays of sufficient finesse is for microwave photons significantly less challenging than in the optical domain.
Free, that is non-interacting photons are in this technology supported by various forms of LC circuits composed of an inductance and a capacitance. In these, current and voltage oscillations occur together with associated oscillations of electric and magnetic fields.
For the appropriate dimensions of these circuits, their elementary excitations are microwave photons with frequencies that remain below the superconducting gap but are large enough to keep thermal excitation numbers vanishingly small at cryogenic temperatures.
The employed LC circuits are thereby built as so called lumped element versions, with dimensions smaller than the wavelength of the explored photons, or as so called coplanar waveguide resonators. The latter can be viewed as a flattened version of a coaxial cable with three conductors patterned in parallel on a chip. Whereas the two outer conductors are grounded, the central conductor carries an electric signal that thus leads to a time and space dependent voltage drop -an electromagnetic field -between the central and outer conductors, see figure 5. The structure thus acts as a waveguide for microwave photons. A discrete spectrum of resonances can then be engineered by cutting the central conductor at two points so that capacitances form and enforce nodes of the current profile leading to anti-nodes of the voltage modes.
V ac Figure 5. A coplanar waveguide resonator formed by three superconducting lines on a chip. The resonator can be excited via a microwave tone applied to the central conductor, which is intersected by two capacitances that play a similar role as the mirrors in a Fabry-Pérot cavity. The voltage profile of a resonance mode is indicated by the sinusoidal lines.
As LC circuits are harmonic oscillators, their excitations do not interact with each other. For applications as quantum simulators, nonlinear elements in the circuit are thus highly desirable. These are provided by Josephson junctions and the circuit elements featuring one or multiple Josephson junctions are usually called a superconducting qubits [176,65]. There are two physical processes that determine the physics of a Josephson junction and hence a superconducting qubit, the Coulomb interaction between Cooper pairs at both sides of the junction that is determined by the junction's capacitance (including a possible shunt capacitance) and the tunnelling of Cooper pairs through the junction as quantified by the Josephson energy.
In the so called charging regime, as generated by a low capacitance junction, the energy eigenstates of the qubit are characterised by the difference in the number of Cooper pairs at both sides of the junction. Here the qubit with its strongly anharmonic spectrum can be interpreted as a matter component, playing a similar role as a two-level atom coupled to optical photons in the approaches discussed in section 2. Such scenarios have been investigated intensively following an influential experiment in 2004 by Wallraff et al. [252] that showed coupling between a superconducting qubit in the charging regime and a coplanar waveguide resonator, with a very high ratio of coupling strength to dissipation rates. Due to the analogies to cavity quantum electrodynamics (QED), this line of research was coined circuit-QED [231] and initially explored scenarios with direct analogies in optical cavity QED [119,84,22,32,152]. The strong coupling coupling strength between qubit and resonator, that is achieved in these setups, results from the large dimensions of the superconducting qubits (∼ 1 µm) leading to large dipole moments and from the strong confinement of the electromagnetic field between the superconducting wires of the coplanar waveguide resonators.
Superconducting qubits have been considered in various forms. Besides the regime where their eigenstates are mostly determined by the charge degree of freedom, so called phase qubits and flux qubits have been investigated intensively. Rather than discussing all these types in detail we here refer the reader to excellent reviews of this matter [176,64,51,197].
To improve robustness against dephasing noise induced by fluctuating background charges on the chip, novel designs for superconducting qubits have been developed within the last decade. The currently most prominent version is the transmon qubit [149], where the charging energy due to Coulomb interaction of Cooper pairs is reduced by shunting the junction with a large capacitance. Further improvement of the coherence times for transmon qubits coupled to superconducting resonators has recently been achieved by building three-dimensional resonators [196,142]. This design is based on a qubit only made out of two superconducting islands connected by a Josephson junction that is kept inside a threedimensional cavity machined out of aluminium which becomes superconducting at the employed cryogenic temperatures. Experiments have already successfully explored quantum many-body physics with multiple qubits in one three-dimensional resonator, see Sec. 3.5.1. Before reviewing quantum simulator designs in this technology we briefly discuss the quantum description of superconducting circuits.

Quantum Theory of Circuits
An excellent introduction to the quantum theory of superconducting circuits has been written by Devoret [63], so that we here merely summarise the main aspects for completeness.
The dynamics of a superconducting electronic circuit can be described in terms of so called node variables [63]. To this end, the description is reduced to a set of independent variables by eliminating the remaining variables via Kirchhoff's laws, which state that the sum of all voltages around a loop and the sum of all currents into a node should be zeros as long as the flux through the loop and the charge at the node remain constant. A convenient choice is then to associate to each node of the network a variable φ(t) = t −∞ dt V (t ), where V (t) is the voltage drop between the considered node and the ground plane. φ has dimensions of a magnetic flux and for superconducting circuits can be linked to the phase of the Cooper pair "condensate" ϕ via the relation ϕ = φ/ϕ 0 , where ϕ 0 = /(2e) is the reduced quantum of flux with e the elementary charge.
Using this language, the energy of an element between nodes j and j + 1 of a circuit network is given by the expressions, 2 for an inductance L and −E J cos(ϕ j − ϕ j+1 ) for a Josephson junction with Josephson energy E J (A dot denotes a time derivative.). If two identical Josephson junctions are included in a closed ring, such as in a superconducting quantum interference device (SQUID), they behave like a single effective junction for currents across the ring, where the Josephson energy of the effective junction can be tuned via a magnetic flux threaded through the ring [176].
A quantum theory for a circuit can be derived in the canonical way. One first sets up a Lagrangian L such that its Euler-Lagrange equations are identical to the classical equations of motion for the currents in the circuit. This Lagrangian is then Legendre transformed into a Hamiltonian by introducing canonical momenta π j for each node variable via π j = ∂L ∂φj [63,191] and the theory is quantised by imposing canonical commutation relations [ϕ j , π l ] = i δ j,l .
Rather than in a chronological order, we here review the work on employing networks of superconducting circuits for the simulation of quantum many-body physics by starting with those models that have been of central interest in the optical domain as well, in particular the Bose-Hubbard and Jaynes-Cummings-Hubbard models.

Bose-Hubbard models
For a lattice of coplanar waveguide resonators, that each couple to a transmon qubit, see figure 6, it can be shown that the dynamics of polaritons, formed by a superposition of resonator and qubit excitations, is described by a Bose-Hubbard Hamiltonian for suitable parameters [160,159]. Since transmon qubits feature a moderate non-linearity due to their large ratio of E J /E C , where E C = e 2 /(2C Σ ) is the charging energy and C Σ the total capacitance of the qubit to ground, they can be modelled by a nonlinear oscillator. A single lattice site of the envisioned architecture is thus described by the Hamiltonian where a † (a) and q † (q) are creation (annihilation) operators of the resonator and qubit and ω q = −1 √ 8E C E J is the transition frequency of the qubit. The coupling g between resonator and qubit typically greatly exceeds the dissipation rates for both modes. If moreover, this coupling is stronger than the nonlinearity of the transmon qubits, g > E C the qubit and resonator modes hybridise and two species of polaritons with annihilation operators p + and p − of the same form as in equation (9) become the elementary excitations of the system [160], see also section 2.1.4.
Neighbouring resonators can for example be coupled via a mutual capacitance C J , see figure 6, leading to the energy, where ∆V j,j+1 is the voltage difference across the coupling capacitance. In a rotating wave approximation, this coupling leads to frequency shifts for both coupled resonators and a tunnelling of photons between both resonators at a rate J a . In terms of the polaritons p + and p − , c.f. equation (9), the Hamiltonian describing the circuit reads where H p± are Bose-Hubbard Hamiltonians as in equa- j describes a density-density interaction between both species with U +− = E C sin 2 (θ) cos 2 (θ). In the derivation of equation (25) the difference between the frequencies of both polariton species has been assumed to greatly exceed all interaction strength and tunnelling rates so that all processes that would lead to a mixing of both species are strongly suppressed, and a rotating wave approximation has been applied.
Another approach to Bose-Hubbard physics in superconducting circuits [158] considers coplanar waveguide resonators that have been made nonlinear due to Josephson junctions or dc-SQUIDs inserted in their central conductors at the location of a current node of the bare resonator [31,158], see figure 7 for a sketch. This approach leads to a Bose-Hubbard model for a single excitation species. The approach can also be interpreted as the SQUIDs and resonators being ultra-strongly coupled [62,31,158,30] so that the splitting between the two normal modes in each resonator becomes comparable to their frequencies and mixing processes between excitations in both normal modes are truly absent.
The concept of making a coplanar waveguide resonator nonlinear by inserting a Josephson junction or SQUID into its central conductor at the location of a voltage node or current anti-node can be taken a step further by inserting multiple junctions into the central conductor. When applied to a waveguide or very long resonator, this approach gives rise to an extended "artificial" medium, which is optically nonlinear at the single photon level. The idea has been considered by Leib et al. [161], where a synchronised switching of the architecture's normal modes from a low excitation quantum regime to a highly excited classical regime has been found.

Chains of transmission line resonators
Circuit QED o↵ers ample possibilities to couple several resonators to a network as there are basically no limitations on the topology and geometry of networks of transmission line resonators [20,21], except for constraints imposed by the detection circuitry and by space on the chip. In figure 1 we display two specific configurations of one-dimensional networks (chains) of transmission line resonators. In figure 1 a), the resonators are coupled at both ends [25], which form a small capacitance and thereby enable photon hopping between adjacent resonators. In this way, also two-dimensional lattices can be formed by coupling more than two resonators at their ends [21]. The advantage of this coupling scheme is scalability whereas a drawback results from the fact that one can only probe the output fields at the borders of the entire network [25]. However, transmission line resonators can also be coupled by reducing the distance between their central conductors for a certain length as shown in figure 1 b). Depending on the exact position of the convergence, and the specific mode under consideration, the coupling between the two resonators is in general both capacitive and inductive [15,14]. The advantage of this coupling scheme is that each transmission line resonator can be probed individually. This advantage however only holds for one-dimensional networks.
For both designs the coupling constants are typically smaller than the frequencies of the field modes and therefore a rotating wave approximation is applicable. The coupling Figure 7. Sketch of array of coplanar waveguide resonators (only central conductors are drawn), that are nonlinear due to a Josephson junction (shunted by a capacitance) inserted into their central conductor. All resonators couple to in-and output lines on the top or bottom. The mutual coupling between the resonators is generated by bringing them into close proximity at suitable points so that a mutual capacitance and inductance forms [23]. Figure reproduced from [158].

Non-local interactions
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 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 (2) , and an expansion up to 4-th order in ϕ 1 − ϕ 2 leads to with C the capacitance in each resonator and C J the capacitance that shunts the dc-SQUID. The coefficient α 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. Lattice elements coupled by Josephson junctions have furthermore been considered for the simulation of Anderson and Kondo lattices [88] and tunable coupling elements [23].

Relation to Josephson junction arrays
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 manybody models, such as the Jaynes-Cummings-Hubbard model.

Jaynes-Cummings-Hubbard models
If the nonlinearities of the employed qubits are larger than all other interaction and tunnelling processes in the circuit network, its dynamics can no longer be approximated by a Bose-Hubbard model for polaritonic excitations, but is described by a Jaynes-Cummings-Hubbard model as given in equation (14) [229,120]. In this regime, hg E C and the qubit is approximated by a two-level system, b → σ − and b † → σ + in equation (23), which gives rise to a Jaynes-Cummings model describing the individual lattice site. Together with the tunnelling of microwave photons between adjacent resonators as described in equation (24), this leads to a Jaynes-Cummings Hubbard model for the entire lattice.
The approximation of a circuit QED system consisting of a transmon qubit and a resonator by a Jaynes-Cummings Hamiltonian has been investigated experimentally, where it was possible to show the characteristic scaling of its nonlinearity with √ n, where n denotes the number of excitations in the system [84]. In contrast to atoms, the approximation of a superconducting qubit as a two level system needs to be considered with much more care, in particular for transmon qubits, where the increased robustness against charge noise comes at the expense of a slightly reduced nonlinearity as compared to charge qubits. Indeed the experiment [152] also found corrections to the approximations by a two-level system due to the finite nonlinearity of the qubit. A Kagomé lattice of coupled coplanar waveguide resonators that each interact with a superconducting qubit such that individual lattice sites can be described by a Jaynes-Cummings model has been investigated by Koch et al. [147,120]. The approach approach considered three-port coupling elements that are capacitively connected to the three coplanar waveguide resonators that meet at each vertex. These couplers lead to time-reversal symmetry breaking, which is a prerequisite for accessing certain classes of quantum many-body states such as fractional quantum Hall states. We discuss quantum simulators for this class of systems in more detail in section 5.
A dimer of two capacitively coupled coplanar waveguide resonators that are each capacitively coupled to a superconducting qubit has been investigated by Schmidt et al. [229] to explore a localisationdelocalisation transition from a self-trapped phase to an oscillating phase as the interaction strength between qubits and cavity modes crosses a critical value. Discrepancies between the classical and quantum prediction for this critical value have been resolved in [233]. This transition has later been observed in an experiment by Raftery et al. [204], c.f. Sec. 2.1.8.
Moreover approaches to simulating spin systems in coupled arrays of circuit cavities and superconducting qubits have been put forward recently [150,9].

Digital Quantum Simulations
In this review we focus our attention on quantum simulations where an experimentally well controllable system is tuned such that it emulates a specific Hamiltonian. This approach to quantum simulation is termed analog quantum simulation. In contrast, one can also employ a device where specific unitary operations can be performed with high precision. A sequence of such operations can then lead to the same dynamics as the target Hamiltonian of a quantum simulation as can be seen via Trotter's formula. This approach of digital quantum simulation is well amenable to devices that have been designed to implement quantum gates and was recently demonstrated in superconducting circuits to simulate spin [222] and fermionic [20] Hamiltonians in one dimension as well as molecular energies [192]. The theory for these approaches was worked out in [113] and [10], see also [183] for a related theoretical and [21] for related experimental work. Other approaches considered the quantum simulation of the regime of s so called ultra-strong light-matter coupling [62,208,209,238] in driven systems [16].

Experimental Progress
The experimental progress towards assembling and controlling larger and larger networks of superconducting circuits has been remarkable in recent years. Whereas most of the efforts are aiming at implementations of quantum information processing tasks, quantum simulation applications have received increasing interested lately.
For quantum computation applications, threequbit gates have been demonstrated in 2012 [206,173,79]. More recently, larger networks of coherently coupled qubits (up to 9 at the time of writing) with performance fidelities suitable for surface code computation [19] and state preservation by error detection in the repetition code [136] have been shown. To boost the scalability of such networks, multilayer structures are now being built [34].
For quantum simulation of quantum manybody systems, very low disorder (below 10 −4 ) has been shown for 12 coupled coplanar waveguide resonators on a Kagomé lattice [246] and large resonator lattices (more than 100 resonators) have been built [120].
Weak localisation has been simulated in a multiple-element superconducting quantum circuit [48] and topological phases together with transitions between them have been measured via the deflection of quantum trajectories with two interacting qubits [216]. Moreover digital quantum simulations of spin models including their mapping to non-interacting fermions [222,20] have been shown and a novel quantum simulation concept employing an experimental generation of Matrix Product States [193] has been demonstrated for the Lieb-Liniger model [71]. There also have been first few-site realisations of analogue quantum simulation devices for Bose-Hubbard chains which may form building blocks of larger scale quantum simulators with interacting microwave photons.
3.5.1. Few-site Bose-Hubbard chains 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 Sec. 4.1.1. An extension to more qubits in a common resonator for exploring dipolar spin models has been considered in [56].

Quantum many-body dynamics and phase diagrams
Besides their application as quantum simulators, the approaches to generate effective quantum many-body Hamiltonians discussed in this review, also call for investigations of the many-body phenomena they give rise to. An important question is whether the phase diagrams for the approaches as discussed in sections 2 and 3 show any deviations from those of the target models or what the phase diagrams of the novel models motivated by these approaches are? We discuss these questions in this section, see also related reviews by Tomadin et al. [242], Schmidt et al. [230] and Le Hur et al. [155]. As most of the interest was initially focussed on equilibrium regimes, we begin our discussion with these.

Equilibrium studies
Initial investigations of the phase diagrams of effective Hamiltonians for interacting photons or polaritons in arrays of coupled cavities [93,212,148] considered equilibrium scenarios by introducing a chemical potential, the physical realisation of which remains an open question, see section 4.1.1 for further discussion of this aspect. For the Jaynes-Cummings-Hubbard model, these works addressed the immediate question whether the phase diagram would show any differences to the Bose-Hubbard case due to the microscopic differences related to the two component nature of the model with atomic excitations and photons. Before discussing the phase-diagrams of the most prominent models, we here first elaborate further on the absence of a chemical potential for photons and approaches to create such a potential.

Engineering a chemical potential for photons
Photons can typically only be trapped for limited times in optical or microwave resonators.
The achievable trapping times are moreover comparable or even below the experimental equilibration timescales. 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 [133,135,143,57,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 coulping 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 [226,99]. 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.

Mean-field calculations for three dimensions
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 ψ = a j as a superfluid order parameter, a mean-field decoupling was performed in the Hamiltonian (14) and a chemical potential µ was added to get where z denotes the number of neighbouring lattice sites, i.e. z = 6 for three dimensions. Depending on whether the value of ψ that minimises the expectation value ofH 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 ω 0 equals the cavity resonance frequency ω C , ∆ = ω 0 − ω 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 field-theoretic 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].
happens for ͑ − ͒ Ͼ 0. nstable: adding more pors the total energy. This n additional mechanism riton number, is unphysivalues of ͓−ϱ , ͔ in similar instability in the it of large photon hopimit: Õ g š 1 1, the Hamiltonian H hop ing, and the latter may be approximation, we concoupling is dropped and ecouple completely. The f the atomic contribution e ground states͒ and the onic tight-binding model As usual, the tight-binding model may be diagonalized in terms of single-particle Bloch waves. For the relevant cases of a two-dimensional ͑2D͒ cubic lattice and a 2D honeycomb lattice, the resulting energy dispersions are given by for the cubic lattice ͓25͔, and for the honeycomb lattice ͑a denotes the lattice constant͒ ͓26͔. In the energy dispersion of the honeycomb lattice, the two signs refer to the lower and upper ͑ and ‫ء‬ ͒ bands. Independent of the specific lattice type, the bosonic ground state is obtained by N-fold occupation of the k = 0 ͑͒ state with corresponding energy the JC lattice model in the resonant case ⌬ =0. ͑a͒ Ground states of the JC lattice system in the atomic The ground state of the system is given by a product state of Jaynes-Cummings eigenstates on each cal potential, the system assumes either the state ͉0͘ j or one of the antisymmetric states ͉n−͘ j . −͘ mark the onset of superfluidity, occurring for finite photon hopping Ͼ 0. The onset points, for − ͱ n + 1, become dense as approaches . For Ͼ , the system becomes unstable. Degeneracies negative and positive detunings ⌬ ͑solid curves͒, except for the lowest degeneracy between ͉0͘ j and the dashed curve. ͑b͒ Mean-field phase diagram of the resonant JC lattice system as a function of nd photon-hopping strength . The color/gray scale shows the magnitude of the order parameter insulating phases ͑denoted "MI"͒ with = 0 and fixed number n = 0 , 1 , . . . of polaritons per site, and The phase boundary can be obtained analytically ͓cf. Eq. ͑20͔͒ and is marked by a black curve. For th, the system becomes unstable with respect to the addition of polaritons. A crude estimate of the z c depicted by the white dashed curve. In the hatched region close to − = 0, numerical results are r the maximum photon number. ͑c͒ Improved numerical results near − =0 ͓range of − and ͔ can be obtained by employing a "sliding" truncation ͑see text͒ centered at the photon occupation RANSITION OF… PHYSICAL REVIEW A 80, 023811 ͑2009͒ 023811-3 Figure 9.
Ground state phase diagram of the Jaynes-Cummings-Hubbard model for ∆ = 0 as obtained from equation (28) for a three dimensional lattice with z = 6. Note that the axes labels use a different notation with β instead of g and κ instead of J JC . Figure reproduced from [148] 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].

4.1.3.
Calculations for one and two dimensions Whereas mean-field approaches are expected to become increasingly accurate for higher and higher lattice dimensions, ground states and dynamics of onedimensional 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.
j3i j ! j2i j is driven by a classical coupling field with Rabi frequency ; the cavity mode of frequency ! cav couples the j1i j ! j3i j and j2i j ! j4i j transitions with coupling constants g 1 and g 2 ; the parameters and account for the detunings of levels 3 and 4, respectively. The atomic part of the system wavefunction for the ith cavity can be fully characterized by the number of atoms in each of the four possible states: fjn 1 ; n 2 ; n 3 ; n 4 ig, with P 4 i1 n i N. The total number of photons plus the number of atomic excitations in the whole system (where states j2i j , j3i j count for one excitation, while j4i j counts for two excitations), is a conserved quantity. Hereafter we assume g 1 ' g 2 g and define the relative atomic detuning w ÿ .
Mott insulator.-The phase diagram of the coupled cavity system is characterized by two distinct phases [3][4][5]: the Mott Insulator (MI) is surrounded by the superfluid (SF) phase. In the MI polaritons are localized on each site, with a uniform density n pol =L, where n pol is the total number of polaritons in a system of L cavities; there is a gap in the spectrum, and the compressibility @=@ vanishes. A finite hopping renormalizes this gap, which eventually vanishes at t . The phase boundaries between the two phases can thus be determined by evaluating, as a function of the hopping, the critical values of at which the gap vanishes. Our data have been obtained by means of the density matrix renormalization group (DMRG) algorithm with open boundary conditions [16]. In numerical calculations, the Hilbert space for the on-site Hamiltonian is fixed by a maximum number of admitted photons n phot max . We chose n phot max 6 for model I and n phot max 4 for model II; we also retained up to m 120 states in the DMRG procedure, such to guarantee accurate results, and checked that our data are not affected by increasing n phot max . We simulated systems with up to L 128, and up to N 5 atoms per cavity [17]; the asymptotic values in the thermodynamic limit have been extracted by performing a linear fit in 1=L. By combining these results with strong coupling perturbation theory [18] we were able to locate the phase boundaries for all values of N. Most of this Letter is devoted to the case ! for model I and 0 for model II. These regimes could not be accessed by the perturbative approach of Ref. [3].
Let us start with zero photon hopping (t 0). For model I, at fixed N, there exists a value I of the detuning I ! ÿ such that, for I > I , the width of the lobe with a polaritonic density N is greatly enhanced with respect to the other lobes. We estimate I numerically and find a scaling I N p . For model II, at a given relative atomic detuning ! > 0 the situation is similar to model I, where the resonating lobe with N is much larger than the other lobes, if < . In the opposite case, ! < 0, some of the lobes disappear.
For model I, numerical data at finite photon hopping for features emerge in the structure of the lobes. In both models, for fixed N, contrary to Bose-Hubbard model, the critical values t of the hopping strength at which the various lobes shrink in a point are not proportional to the lobe width at t 0. Furthermore, the ratio between the upper and the lower slopes of the lobes at small hopping is greater than the one predicted in Ref. [18]; this discrepancy disappears on increasing the number of atoms inside the cavity. In terms of an effective Bose-Hubbard model, this may be understood as a correlated hopping of the polaritons, i.e., the hopping depends on the occupation of the cavity.

H Y S I C A L R E V I E W L E T T E R S
week ending 2 NOVEMBER 2007 Figure 10.
Ground state phase diagram of the Jaynes-Cummings-Hubbard model for ∆ = 0 for a one-dimensional lattice. Note that the axes labels use a different notation with β instead of g and t instead of J JC . Left (right) column for N = 1 (N = 2) atoms per cavity. The bottom row shows the compressibility κ = ∂n/∂µ, where n is the average number of excitations per lattice site. Figure reproduced from [212].
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 [199,117,116], 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 strongcoupling 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].

Non-equilibrium explorations
In contrast to ultra-cold atoms, photons or polaritons are only trapped for considerably shorter times in the samples. On the other hand photons can be produced and injected into the device at low 'cost' and in large numbers. As moreover a chemical potential for photons does not appear in nature but needs to be carefully engineered, c.f. section 4.1.1, it appears to be much more natural and feasible to explore quantum many-body systems of interacting photons in driven scenarios where continuous or pulsed inputs counteract dissipation. We first consider pulsed input drives and discuss continuous driving schemes afterwards.

Pulsed driving fields
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 nonequilibrium dynamics, it was found that the rescaled superfluid order parameter ψ/ √ 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.

Continuous driving fields and driven-dissipative regimes
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 for 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 ω d of the input fields. In this frame, the frequency of a photon in each cavity (resonator), ω r , is replaced by the detuning between ω r and the frequency of the drive ω r → ∆ r = ω r − ω d and similarly the transition frequency of emitters in the resonators is replaced by their detuning from the drive frequency, ω e → ∆ e = ω e − ω d . A field that continuously drives the cavity modes is, in this rotating frame, described by an additional term in the Hamiltonian, where Ω j is the amplitude and ϕ 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 ughly be Þ is phase t cavities, alyze the mode is mode can ies U. For e state of coherent tions. We k;ð=2Þ þ plex numl quantum equation the steady field, n ¼ ch has study its particle statistics, which can be calculated via its characteristic function [15]. For the first order coherence between modes k and p we obtain hB y k B p i ¼ k;ð=2Þ p;ð=2Þ Nn þ k;p m, whereas the density-density 113601-2 Figure 11. Density of photons in the classical background (left) and ratio of photon densities in the quantum fluctuations and classical background (right), for relative phase of π/2 between adjacent driving fields. A semiclassical description is justified for small values in the right plot. Reproduced from [107].
described by master equations of the froṁ 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.

Driven-dissipative 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 tunnelcoupled 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.
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, where j and l label lattice sites, have been investigated by Hartmann [107] via a numerical integration of equation (30) with Matrix Product Operators [232,111]. A spatial 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 (2) j,j 1, slight bunching for neighbouring sites, g (2) j,j+1 > 1, and anti bunching for further separated sites, g (2) j,l < 1 for |j − l| > 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].
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 10, and andñð8Þ g ð2Þ r ð8; jÞ The crysng J only gly pror findings rities, U,     [1]). We see in Fig. 1 that there are regions with 1 or 2 stable solutions, but also regions with no stable homogeneous solution [shown in red (labeled ''0'')]. Notice however that, within our mean field approach, we have direct access only to spatially uniform solutions, where all the cavity sites are equivalent. Interestingly, we find that the bistability induced by the coupling between the cavities also appears when the pump 233601-2 Figure 13. Number of mean-field solutions for Ω/∆ = 0.4, γ/∆ = 0.2 and ∆ > 0. Note that the notation ∆ω = ∆ is used in the axes labels. Reproduced from [153].
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 [256,174].
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 meanfield 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) < g (2) (j, j), see figure 14.
In the studies discussed so far, the excitation dissipation has been assumed to be caused by the electromagnetic vacuum throughout. As has been Reproduced from [128].
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.

Driven-dissipative Jaynes-Cummings-Hubbard model
The coherence and fluorescence properties of a coherently pumped driven-dissipative Jaynes-Cummings-Hubbard model were explored by Nissen at 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. A comparative study of the features of drivendissipative 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 The 2LS approximation is appropriate as the nonlinearity Þ of the JCM prevents higher states from being excited ( 1;2 are the lowest energies in the Hilbert space sector with 1 or 2 excitations, respectively).  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.

Regimes of incoherent pumping
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 Sec. 4.1.1 about engineering effective chemical potentials.
For further work on non-equilibrium photon condensation, see [144,236,3], 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. 4.2.6. Transport studies Transport of photons in a waveguide that are scattered at a localised emitter has been theoretically explored via scattering theory [234,165,187], 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 [221,171,115,46,218] 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 wavefunction 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 manybody 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.

Josephson interferomenter
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 antibunching 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.

Calculation Techniques
The modelling of quantum many-body systems is a formidable challenge since the dimension of their Hilbert space grows exponentially in the number of constituents. This scaling renders exact descriptions, including exact numerical approaches infeasible, even for moderate system sizes. For dissipative quantum many-body systems the computational effort is even more dramatic as mixed quantum states need to be considered.
Exceptions to this intractability are quantum systems that do not explore their entire Hilbert space, where numerical optimisation approaches such as DMRG [232] become efficient. For calculating equilibrium phase diagrams DMRG has thus been used in a number of works considering manybody systems of interacting photons [212,213,215]. In turn for calculating dynamics, Matrix Product State representations of DMRG [193] in form of the Time-Evolving Block Decimation (TEBD) [250,251,263] have frequently been applied [110,107,95,128,198,225,202,25].
Alternatively to approximations with states of limited entanglement, one may aim for only obtaining the information of interest about the quantum state of the entire system and try to find accurate and efficient approximations for the sought quantities. Mean-field approaches [131] can be understood as representatives of this strategy as they only predict properties of a single constituent of the many-body system [220,85]. Such approaches, which calculate local quantities but ignore all correlations between subsystems, have been applied in both, the calculation of equilibrium states [93] as well as dynamics [243] and stationary states of driven-dissipative systems [189,153]. In some cases, e.g. some scenarios of strongly interacting Rydberg gases as discussed in section 2.3, only states with low excitation numbers contribute so that full numerics in a strongly truncated Hilbert space provides a good approximation [171].
Yet, as DMRG approaches are limited to onedimensional systems and mean-field techniques are only expected to become accurate in very high lattice dimensions, which often do not correspond to the physical realisations, there is a need for further efficient methods for accurately calculating stationary states of quantum many-body systems. As a consequence a substantial amount of research has recently been dedicated to the development of such methods.
Keldysh path integral methods have been used to explore long-range properties [240,235] and dynamical mean-field theory has been generalised to nonequilibrium scenarios [7]. For solving Lindbald type master equations for their stationary states, del Valle et al. have expanded the resulting equations for correlators in powers of the inter-site coupling [61]. Li et al. have developed a perturbation theory for the Lindbladian including a resummation technique for the perturbations [163,164]. Degenfeld-Schonburg et al. have generalised open quantum system techniques to take dynamical environments into account so that they can describe the interacting constituents of a manybody system in a consistent Mori projector theory (c-MoP) [58,59].
For directly finding the stationary states of a master equation without doing a time integration, variational approaches have been developed. These include a variational expansion around product states [256] and variational Matrix Product Operator approaches [55,181]. Alternative approximations that focus on a restricted part of the Hilbert space in a similar way as Matrix Product State representations by keeping only the dominant eigenvalues of reduced density matrices but can be applied in two-dimensional lattices have been introduced by Finazzi et al. as a "corner-space renormalization method" [83]. Moreover, a dynamical polaron ansatz has been introduced for treating very strong light-matter couplings [66,151].
For experimental investigations of the predicted phase diagrams and phenomena, a crucial question is whether their experimental signatures are accessible in measurements. We therefore discuss some of these signatures in the next section.

Experimental signatures
To clarify whether the phase diagrams and transitions discussed here can eventually be observed in experiments, it is important to determine their signatures in measurable observables. In this context it is natural to consider the photons emitted from the individual resonators. These output fluxes are related to intraresonator quantities via input-output relations [253]. In this way the number of polaritons in each resonator, n l = p † l p l for site l, and the number fluctuations, can be measured. These quantities provide information about the localisation and delocalisation of the polaritons, e.g. in a Mott insulator to superfluid transition. The visibility of interference fringes of photons emitted from the resonators [212,121] moreover provides information about coherences between resonators that may have built up in condensation, thus signalling a superfluid phase, c.f. [29]. The visibility V = max(S) − min(S) max(S) + min(S) (33) of the interference pattern can be expressed in terms of the polariton number distribution in momentum space, where L is the number of sites in a one dimensional array.
In the superfluid regime, V approaches unity whereas it is very small for a Mott insulator. Remarkably, equation (34) also allows to bound the entanglement in a many-body system from below without further assumptions [53,52].
Moreover, coincidence counting measurements of the photons emitted from one or multiple resonators [253] allow to reconstruct g (2) -functions of the polaritons in the resonator for both, coincidences from one resonator leading to g (2) (j, j) and coincidences from separate resonators leading to g (2) (j, l) with j = l. These quantities would reveal the density wave ordering predicted for driven-dissipative Bose-Hubbard models with phase off-set in the driving fields [107] or with cross-Kerr interactions [128]. Higher order correlation-functions [60] could then be measured to reveal further properties of the investigated quantum many-body states.
In the next section we turn to review a recent development that has been largely triggered by the ample possibilities for engineering the band structure for photons in many devices. This is the realisation of artificial gauge fields for the quantum simulation of many-body systems under their influence.

Artificial gauge fields
Charged particles moving in magnetic fields are an important paradigm in quantum mechanics and give rise to intriguing phenomena including the celebrated quantum Hall effect [260,239].
The increasing understanding of the geometric foundations of the quantum Hall effect also led to the discovery of topological insulators and topological superconductors which show exotic properties routed in the topological structure of their electron bands [201].
The first approaches to exploring photons that are subject to an artificial gauge field were put forward by Haldane and Raghu [105,205], who considered a hexagonal array of dielectric rods leading to two dimensional photonic bands that show an appreciable Faraday effect leading to a breaking of time-reversal symmetry. The photonic bands can then be characterised by their Chern number [241], which is a topological invariant. Moreover, when two "materials" with different Chern numbers are joined, a chiral state emerges, that propagates along the interface in one unique direction only. These so called edge states are robust with respect to scattering at impurities provided the associated interaction energy is smaller than the gap to neighbouring bands. For photons these properties allow to engineer one-way waveguides which are free of backscattering.
For microwave photons, an approach to engineer a time reversal symmetry breaking in networks of superconducting circuits has been introduced by Koch et al. [147], see also [191]. In this scheme, a circulator element formed by a superconducting ring intersected by three Josephson junctions that is threaded by a flux bias couples three coplanar waveguide resonators.
The absence of backscattering at impurities for edge modes has been considered by Hafezi et al. [103] for engineering robust optical delay lines. In their approach, toroidal micro-cavities coupled by loops of tapered optical fibre with different path length for different directions of propagation are considered to emulate the effect of an artificial gauge field similar to a perpendicular magnetic field for the photon dynamics. An alternative approach to generating effective gauge fields via a dynamical modulation of the photon tunnelling rate at difference between the oscillation frequencies of the two adjacent resonators has been presented by Fang et al. [76]. The emerging gauge fields do in both concepts not need to be spatially uniform and can thus be employed to guide photon propagation and implement waveguides [168].
Experimentally, the absence of backscattering in chiral edge modes of photons has first been investigated in the microwave regime in photonic crystals of a square lattice geometry [255]. Profiles of edge modes have then been imaged in silicon photonics devices [101]. More recently, topological invariants in such systems have been measured via the shift of the spectrum as a response to a quantum of flux inserted at the edge [185]. This technique allows to access the winding number which is directly related to the Chern number via the bulk-boundary correspondence. An alternative approach to engineering an artificial gauge field has been explored with a continuous drive on one site of a Bose-Hubbard dimer [211]. For further details on the physics of artificial gauge fields for propagating photons, we refer the interested reader to the review by Hafezi [97].
Whereas experimental progress has so far been mostly made with samples where photon-photon interactions via optical nonlinearities can be neglected, theory research has also addressed the intriguing regime where strong artificial gauge fields and strong effective photon-photon interactions coexist. Most notably this leads to regimes for fractional quantum Hall physics.
A first proposal for generating a fractional quantum Hall regime in coupled photonic resonators was put forward by Cho et al. [49]. The approach considers optical cavities doped with single atoms with a lambda shaped level structure featuring two metastable states. By driving two photon transitions with external lasers with nonuniform relative phases, an artificial gauge field is engineered in this strongly nonlinear system. This driving pattern leads to a bosonic version of a fractional quantum Hall regime, with even number of magnetic fluxes per excitation.
In circuit-QED systems in turn, three-body interactions in the presence of artificial gauge fields have been explored, which lead to Pfaffian states [98]. The artificial gauge field is in circuit-QED architectures implemented via externally modulated coupling SQUIDs as discussed in section 3.2.1. The adjacent building blocks have different transition frequencies and the coupling SQUID is driven by an external flux which oscillates at a frequency that equals the difference of the transition frequencies of its two neighbouring lattice sites. Relative phases between the drives at several coupling SQUIDs then encode an applied gauge field. If in turn, the photon tunnelling is dynamically modulated at the sum of the transition frequencies of the adjacent lattice sites, a Kitaev spin chain showing Majorana zero modes can emerge [18].
The dynamical modulation of a coupling circuit for engineering an artificial gauge field has recently been demonstrated in an experiment by Roushan et al. [217]. Here three superconducting transmon qubits were coupled in a ring via tunable couplers [47] that were modualted by oscillating magnetic fluxes to simulate the effect of a perpendicular gauge field. In this setup, that thus combines an artificial gauge field with local interactions provided by the qubit nonlinearities, single excitations and excitation pairs were circulated around the ring in a controlled way.
The emergence of fractional quantum Hall physics in driven-dissipative lattices of coupled photonic resonators has been investigated by Umucalilar et al. [245] by showing that the stationary states of such driven dissipative systems can, when projected on a specific excitation number, approximate Laughlin wave-functions with even number of magnetic fluxes per excitation. The overlap of the stationary states of such cavity lattices with a Laughlin state was analysed with an efficient approximation for low excitation numbers in [104].
The potential of engineered dissipation for stabilising topological states by coupling the cavity lattice to two-level systems with fast dissipation, has been explored by Kapit et al. [132]. Moreover a quantum simulator for topological order with superconducting circuits has been proposed in [223] The theory research on models with artificial gauge fields has recently also been pushed further to explore the generation of non-Abelian [184] and dynamical gauge fields in lattice gauge filed theories [178,180]. For the latter, the gauge fields are not set by an external current or voltage source by formed by dynamical degrees of freedom of the network.

Summary and Outlook
After optical nonlinearities at the single photon level, i.e. effective interactions between individual photons, have been realised in single micro-cavities, coupling several optical or microwave resonators to form a network has become a new research goal. In parallel avenues to generate effective interactions between individual photons in extended one-dimensional volumes have been considered. Whereas the initial work in the research field was mostly of theoretical nature, technological advances in recent years have now matured the experimental platforms to such an extent that an increasing number of experimental investigations are to be seen in the coming years. This development has two exciting perspectives.
From a scientific perspective, quantum many-body systems of interacting photons can be expected to exhibit a wealth of new quantum many-body phenomena as they naturally operate under driven, non-equilibrium conditions which are different to the equilibrium scenarios usually explored in quantum many-body physics. From a technology perspective, photons are the most suitable carrier for transmitting information over long distances as they are largely immune to environmental perturbations. As optical nonlinearities in conventional media are weak it is therefore of great importance to conceive means of making individual photons interact with each other at multiple nodes of a network to make them suitable for information processing. The coupled micropillars used in our studies (Fig. 1a) are obtained by dry etching of a planar semiconductor microcavity with a Rabi splitting of 15 meV at 10 K, the temperature of our experiments (see Methods). Each individual pillar has a diameter of 4 µm and presents a series of confined polaritonic states 21-23 with a lifetime ⌧ ⇠ 33 ps. The centre-to-centre separation of 3.7 µm results in a tunnel coupling of the lowest energy confined polaritonic states (ground state) of J = 0.1 meV. This double potential well system can be described by equations (1a) and (1b) with the addition of a phenomenological decay term 9,24 , i(h/2⌧ ) L(R), accounting for the polariton losses due to the escape of photons out of the cavity, and with a positive value of U coming from the polariton-polariton repulsive interactions 25 .
The Madelung transformation L,R(t ) = p NL,R(t )e i✓L,R(t ) allows us to rewrite equations (1) in their dynamical form 4 : where z(t ) = NL NR NT is the population imbalance between the two micropillars (with NT = NL + NR the total population), and (t ) = ✓L(t ) ✓R(t ) is the phase difference. E 0 L E 0 R is the energy difference between the ground states of the left and right micropillars in the absence of coupling, negligible in our case as the two pillars are nominally identical. The second term on the right hand side of equation (2b) contains all the features due to interactions, and it is the only one that is affected by the polariton finite lifetime.
To study the different Josephson regimes we excite the system with a 1.7 ps pulsed laser resonant with the ground state energy of the micropillars (⇠780 nm). The spectral width of the laser (0.4 meV > J ) allows us to initialize the system in a linear combination of bonding and antibonding states 26 (B and AB in the figures). We use a 10 µm wide Gaussian excitation spot, covering the whole molecule. In this geometry, (0) = 0 and the initial population imbalance z(0) can be tuned by shifting the spot with respect to the centre of the molecule.
The light emitted from a single polaritonic molecule is collected with a microscope objective in reflection geometry and analysed in energy and time by means of a spectrometer coupled to a streak camera. To avoid the strong reflection of the laser beam in our detectors, we select the linear polarization of emission that is perpendicular to that of the excitation (parallel to the molecule long axis). The molecules present an intrinsic linear polarization splitting along a non-trivial direction which slowly rotates the polarization direction of the injected polaritons 27 . This allows us to measure in cross-polarized detection the energy, population imbalance z(t ) and phase difference (t ) between the two sites (see Methods). The polarization splitting is of the order of ⇠40 µeV, resulting in a polarization rotation that is much slower (>70 ps) than the Josephson oscillation timescales in our system (up to 21 ps).
At low excitation density (2J UNT ⇡ 0), interaction effects are negligible and the dynamics of the system is entirely dominated by the tunnel coupling. This is the situation presented in Fig. 1a The coupled micropillars used in our studies (Fig. 1a) are obtained by dry etching of a planar semiconductor microcavity with a Rabi splitting of 15 meV at 10 K, the temperature of our experiments (see Methods). Each individual pillar has a diameter of 4 µm and presents a series of confined polaritonic states 21-23 with a lifetime ⌧ ⇠ 33 ps. The centre-to-centre separation of 3.7 µm results in a tunnel coupling of the lowest energy confined polaritonic states (ground state) of J = 0.1 meV. This double potential well system can be described by equations (1a) and (1b) with the addition of a phenomenological decay term 9,24 , i(h/2⌧ ) L(R), accounting for the polariton losses due to the escape of photons out of the cavity, and with a positive value of U coming from the polariton-polariton repulsive interactions 25 .
The Madelung transformation L,R(t ) = p NL,R(t )e i✓L,R(t ) allows us to rewrite equations (1) in their dynamical form 4 : where z(t ) = NL NR NT is the population imbalance between the two micropillars (with NT = NL + NR the total population), R is the energy difference between the ground states of the left and right micropillars in the absence of coupling, negligible in our case as the two pillars are nominally identical. The second term on the right hand side of equation (2b) contains all the features due to interactions, and it is the only one that is affected by the polariton finite lifetime.
To study the different Josephson regimes we excite the system with a 1.7 ps pulsed laser resonant with the ground state energy of the micropillars (⇠780 nm). The spectral width of the laser (0.4 meV > J ) allows us to initialize the system in a linear combination of bonding and antibonding states 26 (B and AB in the figures). We use a 10 µm wide Gaussian excitation spot, covering the whole molecule. In this geometry, (0) = 0 and the initial population imbalance z(0) can be tuned by shifting the spot with respect to the centre of the molecule.
The light emitted from a single polaritonic molecule is collected with a microscope objective in reflection geometry and analysed in energy and time by means of a spectrometer coupled to a streak camera. To avoid the strong reflection of the laser beam in our detectors, we select the linear polarization of emission that is perpendicular to that of the excitation (parallel to the molecule long axis). The molecules present an intrinsic linear polarization splitting along a non-trivial direction which slowly rotates the polarization direction of the injected polaritons 27 . This allows us to measure in cross-polarized detection the energy, population imbalance z(t ) and phase difference (t ) between the two sites (see Methods). The polarization splitting is of the order of ⇠40 µeV, resulting in a polarization rotation that is much slower (>70 ps) than the Josephson oscillation timescales in our system (up to 21 ps).
At low excitation density (2J UNT ⇡ 0), interaction effects are negligible and the dynamics of the system is entirely dominated by the tunnel coupling. This is the situation presented in Fig. 1a of ⇡), a situation also reproduced by simulations of equations (2) (see Supplementary Information). The observed oscillations are probably damped owing to intensity fluctuations of the pump laser, resulting in fluctuations in the recovery time of the oscillations. This regime is different to the case depicted in Fig. 2, where anharmonic oscillations take place around = 0. Here, ⇡-oscillations appear naturally from the untrapping of the condensate from the quasiantibonding mode, which is characterized by = ⇡. The dynamical transition from self-trapping to ⇡-oscillations can be obtained only thanks to the finite particle lifetime and cannot be easily realized in systems such as Josephson junctions with atomic condensates.
In our experiments we have taken advantage of the large particleparticle interactions originating from the excitonic component of polaritons to demonstrate macroscopic self-trapping in a photonic system. This result opens the way to highly nonlinear photonic phenomena such as polarization chaos 13 , or spontaneous symmetry breaking and pitchfork bifurcations 17,28,29 , expected to give rise to highly squeezed macroscopic states 14 . By reducing the size of the system, the many-particle interactions evidenced here could be taken into the regime of single-photon nonlinearities 30,31 . Furthermore, the pair of coupled cavities we have used could be extended to arrays with minimal frequency dispersion, where photon fermionization 32,33 and the emergence of Bose-Hubbard physics [34][35][36][37] have been predicted.

Methods
Sample and experimental set-up. The planar sample was grown by molecular beam epitaxy and consists of a ⌦/2 cavity sandwiched between two distributed Bragg reflectors containing 26 and 30 pairs, respectively, of Al0.95Ga0.01As, Al0.20Ga0.80As ⌦/4 layers. The structure contains 3 groups of 4 GaAs quantum wells of 70 Å width, located at the maxima of the electromagnetic field in the structure, resulting in a Rabi splitting of 15 meV. The quality factor of the etched structure is 16,000, resulting in a photon lifetime of 15 ps. As we perform the experiments in a polaritonic molecule with zero exciton-photon detuning, the polariton lifetime is extended up to 30 ps.
The polaritonic molecules are fabricated by dry etching of the planar structure. The basic building block of the molecule is a single round micropillar, which presents a set of discrete polariton levels arising from three-dimensional confinement. When two such micropillars spatially overlap, polaritons can tunnel from one micropillar to the other. The tunnel coupling is proportional to the overlap and can thus be engineered 8 .
Excitation of the sample is performed with a 1.7 ps pulsed laser (repetition rate 82 MHz) resonant with the bonding and anti-bonding states. Its spectral width is 0.4 meV, which is larger than the bonding-antibonding splitting, thus allowing the preparation of a linear combination of both states 26 . Simultaneously, it is smaller than the energy separation between the ground state and the first excited state in a single micropillar 22 , and also smaller than the energy distance to the polariton reservoir (located 7.5 meV above). In our analysis we can thus ignore the excited states and only consider the ground state coupled modes.
The photoluminescence dynamics are recorded using a streak camera synchronized with the excitation laser, with a time resolution of 4 ps. Recorded images are the result of integration of several million realizations. The excitation pulse is linearly polarized along the long axis of the molecule and detection is performed in the orthogonal polarization. Despite the cross-polarized detection, stray light from the laser prevents access to the population imbalance and phase difference during the first few picoseconds.

Measurement of the phase difference .
To measure the phase difference between the emissions from each micropillar, we use a Michelson interferometer. The real space emission of the molecule interferes with its mirror image at the entrance slit of a streak camera. In this way, we monitor interferences between the emission of one of the micropillars of the molecule and that of the other micropillar. Constructive/destructive interference appears depending on the phase difference between the emitted light from each micropillar, and also from the delay between the two arms of the interferometer. By varying the delay of one of the arms with a piezoelectric stage by 6⇡, we obtain an oscillating interference pattern from which we extract the phase difference (t ) at different delay times 20 t .

Energy-resolved emission.
To measure the energy of the emission, we image the micropillars on the entrance slit of a spectrometer placed in front of the streak camera. In this way we can measure the time evolution of the energy of the bonding (E ⇤ B ) and antibonding (E ⇤ AB ) states. with the associated single-photon nonlinearities. We develop a characterization technique useful for systems with low dissipation that relies on the Jaynes-Cummings nonlinearity, which leads to bistability with a sharp transition to a bright state as an applied continuous microwave tone is swept in power [57][58][59]. The threshold for this transition (above which the bright state behaves linearly) is sensitive to the frequency difference between the uncoupled mode being monitored and the  BECs [21]. The exponential decay of the oscillations later gives way to a superexponential drop in homodyne signal, a signature of the crossover from delocalized to localized behavior: photon escape is a stochastic process, and for a given trial the photon number falls below N c at a random time, with an average time dependent on the initial photon number. When approaching this point, the Josephson oscillations become nonlinear, exhibiting a critical slowing-down [22]. Oscillations of different trials within an ensemble dephase with respect to each other, and individual trials once localized exhibit very rapid collapse; thus, ensemble averages ofÎ andQ die out faster than exponentially, and only trials where NðtÞ ≫ N c continue to contribute to the homodyne signal. Figure 6(b) shows the observed homodyne dynamics for various initial photon numbers, revealing the logarithmic dependence of the critical time to reach the transition on the initial photon number, t c ∼ 1 κ 0 logðN i =N c Þ.
Unlike the homodyne signal, photon number measurement is insensitive to the coherence of field in the monitored cavity. Additionally, the dispersion between individual trials arising from critical slowing-down does not cause the photon number to decay superexponentially. Figure 6(c) compares homodyne measurement to photon number for the same initial condition. For short times, the two signals match, demonstrating a high degree of coherence within the ensemble. In contrast to the superexponential decay of the homodyne signal, we see an exponential decay of the photon number.
The homodyne observation maps out a dynamical phase diagram as a function of initial photon number N i and time, as displayed in Fig. 6(d). An applied drive during an initialization phase where photon-photon interactions are off initiates Josephson oscillations. An undriven and noninteracting region lasting a constant time τ follows, continuing the oscillations, at the end of which the  Fig. 5. With interactions off, the system undergoes oscillations and exponential decay at rate κ, which can be observed for several microseconds. Significantly, the presence of interactions causes superexponential decay of the homodyne signal as N approaches Nc, a signature of crossover into the localized regime. (b) The time of onset of superexponential decay tc shifts with initial photon number, as shown here for initialization pulses with varying drive power lasting 11.5 μs (100=J). The use of a long initialization pulse makes it possible to drive the system to very large initial photon numbers (top to bottom: Ni ≈ 12 000, 3800, 1100, 550, 40), but introduces complications (see Supplemental Material for more details [45]). (c) Directly measuring the photon number (green) reveals that incoherent photons remain in the system after the homodyne signal (blue) has undergone superexponential decay. Oscillations in the photon number can also be observed to die out, as critical slowingdown constrains the envelope of oscillations, finally leaving only exponential decay. Here, Vdrive is 1.15 μs (10=J) and τ ¼ 1 μs.
Background voltages leading to distortion of the signal were removed from the photon number measurement.   [1]. Polaritons are initially generated in the left pillar. (b) Collective oscillations of the polaritons between both pillars as measured by the emitted intensities for moderate initial polariton imbalance. (c) Self trapped regime where the initial polariton imbalance is strong enough to suppress collective oscillations. (d) Set-up of two coupled superconducting resonators that each interact with a transmon qubit [204]. (e) Phase diagram as observed in the experiment [204] with collective oscillations for large initial imbalance that are suppressed as the imbalance decreases. Plots a,b and c adapted from [1] and plots d and e adapted from [204] Figure 6. Chain of coplanar waveguide resonators (central and grounded conductors indicated by grey lines), that each couple to a transmon qubit (balck boxes). Neighbouring resonators are coupled via a capacitance (indicated by the small gaps in the central conductor).