Quantum simulations and many-body physics with light

In the following we describe in detail the initial proposal by Angelakis et al [5] (first appeared in June 2006 and published in 2007), to implement the JCH model using two-level-system-doped CRAs in order to observe the Mott-superfluid transition of photons.

Consider a series of optical resonators coupled through direct photon hopping. The hopping could be generated through an overlap in the photonic modes as we assumed here or using mediating elements such as fibres [48]. The system is described by the Hamiltonian

$$\begin{eqnarray}&&H=\underset{k=1}{\overset{N}{\sum}}\,{{\omega}_{d}}a_{k}^{\dagger}{{a}_{k}}-\underset{k=1}{\overset{N}{\sum}}\,J\left(a_{k}^{\dagger}{{a}_{k+1}}+\text{H}.\text{c}.\right)\end{eqnarray} \tag{ 1 }$$

which corresponds to a series of coupled quantum harmonic oscillators. The photon frequency and hopping rate is ${{\omega}_{d}}$ and J respectively and no nonlinearity is present yet. In order to simulate strongly correlated effects, interaction is necessary and the most natural way to induce it is to dope each resonator with a two level system. The latter could be an atom or a quantum dot or a superconducting qubit. Let us denote $|g{{\rangle}_{k}}$ and $|e{{\rangle}_{k}}$ as its ground and excited states at site k, respectively. The Hamiltonian describing the total system can be divided into three terms: H_free the Hamiltonian for the free resonators and dopants, H_int the Hamiltonian describing the local coupling of the photon and the dopant, and H_hop for the photon hopping between cavities. Explicitly, these read

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{\text{free}}}\,=\,{{\omega}_{d}}\underset{k=1}{\overset{N}{\sum}}\,a_{k}^{\dagger}{{a}_{k}}+{{\omega}_{0}}\underset{k}{{{\sum}^{}}}\,|e{{\rangle}_{k}}\langle e{{|}_{k}}, \end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{\text{int}}}\,=\,g\underset{k=1}{\overset{N}{\sum}}\,{{(a_{k}^{\dagger}|g\rangle}_{k}}\langle e{{|}_{k}}+\,\text{H}.\text{c}.), \end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \,{{H}_{\text{hop}}}\,=\,-J\underset{k=1}{\overset{N}{\sum}}\,\left(a_{k}^{\dagger}{{a}_{k+1}}+\text{H}.\text{c}.\right), \end{array}\end{eqnarray} \tag{ 4 }$$

where g is the atom-photon coupling strength. The local part of the Hamiltonian ${{H}_{\text{free}}}+{{H}_{\text{int}}}$ can be diagonalised in a basis of mixed photonic and atomic excitations. These excitations, known as dressed states or polaritons in quantum optics and condensed-matter communities respectively, involve a mixture of photonic and atomic excitations and are defined by annihilation operators $P_{k}^{(\pm,n)}=|g,0{{\rangle}_{k}}\langle n\pm {{|}_{k}}$. The polaritonic states of the kth atom-resonator system are given by $|n+{{\rangle}_{k}}={{(\sin {{\theta}_{n}}|g,n\rangle}_{k}}+\cos {{\theta}_{n}}|e,n-1{{\rangle}_{k}})/\sqrt{2}$ and $|n-{{\rangle}_{k}}={{(\cos {{\theta}_{n}}|g,n\rangle}_{k}}-\sin {{\theta}_{n}}|e,n-1{{\rangle}_{k}})/\sqrt{2}$, where $|n{{\rangle}_{k}}$ denotes the n-photon Fock state in the kth resonator, $\tan 2{{\theta}_{n}}=2g\sqrt{n}/ \Delta$, and $\Delta ={{\omega}_{0}}-{{\omega}_{d}}$ is the atom-light detuning. These polaritons have the energies $E_{n}^{\pm }=n{{\omega}_{d}}+ \Delta /2\pm \sqrt{n{{g}^{2}}+{{(\Delta /2)}^{2}}}$ and are also eigenstates of the sum of the photonic and atomic excitation operators ${{\mathcal{N}}_{k}}=a_{k}^{\dagger}{{a}_{k}}\,+\,|e\rangle \langle e{{|}_{k}}$ with the eigenvalue n (figure 1).

Zoom In Zoom Out Reset image size

Studying the energy eigenstates of the system consistent with a given number of net excitations per site (or filling factor), it was shown that the ground state of the system reaches a Mott state of the polaritons for integer values of the filling factor. In contrast to the optical lattice systems where the lattice depth is tuned, here the detuning $\Delta ={{\omega}_{0}}-{{\omega}_{d}}$ is varied. The latter results in further shifts of the dressed states from the bare resonances, providing an effective control over effective nonlinearity. The hopping strength J is usually fixed by the fabrication process and is harder to change experimentally. To see how the Mott state arises, it is convenient to rewrite the Hamiltonian in terms of the polaritonic operators (assuming $\Delta =0$):

$$\begin{eqnarray}\begin{array}{*{35}{l}} H = & \underset{k=1}{\overset{N}{\sum}}\,\underset{n=1}{\overset{\infty}{\sum}}\,\left[n\left({{\omega}_{d}}-g\right)P_{k}^{(-,n)\dagger}P_{k}^{(-,n)} \right. \\ &+n\left({{\omega}_{d}}+g\right)P_{k}^{(+,n)\dagger}P_{k}^{(+,n)}+g\left(n-\sqrt{n}\right)P_{k}^{(-,n)\dagger}P_{k}^{(-,n)} \\ &\left. +g\left(\sqrt{n}-n\right)P_{k}^{(+,n)\dagger}P_{k}^{(+,n)}\right] \\ &-J\underset{k=1}{\overset{N}{\sum}}\,\left(a_{k}^{\dagger}{{a}_{k+1}}+\text{H}.\text{c}.\right). \end{array}\end{eqnarray} \tag{ 5 }$$

Remember that the lower of the two energy eigenstates of the local system with a given number of excitations η is the '-' dressed state $|\eta -\rangle$. Therefore, assuming the regime $Jn\ll g\sqrt{n}\ll {{\omega}_{d}}$, one needs only consider the first, third and last terms of the above Hamiltonian for determining the low-lying energy eigenstates of the whole system. The first line corresponds to a linear spectrum, equivalent to that of a harmonic oscillator of frequency ${{\omega}_{d}}-g$. If only that part were present in the Hamiltonian, then it would not cost any extra energy to add an excitation (of frequency ${{\omega}_{d}}-g$) to a site already filled with one or more excitations, as opposed to an empty site.

The third term, on the contrary, provides an effective 'on-site' photonic repulsion, leading to the formation of the Mott state of polaritons as the ground state of the system for a commensurate filling. Reducing the strength of the effective nonlinearity (or the photonic repulsion) through detuning, one can then drive the system to the superfluid (SF) regime. This could be done, for example, by Stark shifting the atomic transitions from the resonator by an external field. The shift being inversely proportional to the detuning, the new detuned polaritons are not as well separated as before, and in the case of large detuning, it costs little extra energy to add excitations (excite transitions to higher polaritons) in a single site. As a result, the system moves towards the SF regime with increasing detuning. Note here that the mixed nature of the polaritons could in principle allow for excitations that are mostly photonic, resulting in a photonic Mott state [74].

Beyond the qualitative description reproduced above, an exact numerical simulation was performed in [5], for a different number of sites (from three to seven) using the Hamiltonian $H={{H}_{\text{free}}}+{{H}_{\text{int}}}+{{H}_{\text{hop}}}$. Note that such finite-size calculations are often employed to study the physics of the transition between Mott and superfluid phases [75] because mean-field theories in low dimensions are not guaranteed to work. We also note here that numerical simulations of the JCH Hamiltonian are more computationally costly than the ones for the BH Hamiltonian due to the existence of both bosonic (photon) and atomic operators in every site.

As discussed in [5], an appropriate order parameter in the JCH lattice is the variance of the polariton number per site, as the usual order parameter employed in other cases—like the expectation value of the annihilation operator in the BH model—is identically zero in a closed system. The variance $\text{var}\left({{\mathcal{N}}_{k}}\right)$ is plotted in figure 2 as a function of ${{\log}_{10}}(\Delta /g)$ for a filling factor of one net excitation per site. For this plot, we have taken the parameter ratio g/J  =  10² (g/J  =  10¹ gives quantitatively similar results), with $\Delta$ varying from $\sim {{10}^{-3}}g$ to $\sim g$ and ${{\omega}_{d}},{{\omega}_{0}}\sim {{10}^{4}}g$. The order parameter for dissipationless atoms and cavities is displayed as solid curves and corresponding curves in the presence of a resonator dissipation and atomic spontaneous emission of $\sim 0.01g$ are depicted with dotted lines.

Zoom In Zoom Out Reset image size

In figure 3 the quantum dynamics characterising the excitation probability of the middle site in a chain of three cavities for different regimes are calculated assuming the middle resonator is initially empty. We observe the blockade of higher than one excitation in the limit of resonant interaction (where the Mott phase is expected) as discussed in figure 2.

Zoom In Zoom Out Reset image size

The phase diagram of the JCH model has first been studied in 2D using mean field theory [7] and in 1D using a numerical method based on DMRG techniques [8]. Figure 4(top) shows slices of the superfluid order parameter $\psi =\langle {{a}_{i}}\rangle$ as a function of the photon hopping and the chemical potential for different values of atom-photon detuning. A triangular lattice was assumed and the usual decoupling approximation $a_{i}^{\dagger}{{a}_{j}}=\langle a_{i}^{\dagger}\rangle {{a}_{j}}+a_{i}^{\dagger}\langle {{a}_{j}}\rangle +\langle a_{i}^{\dagger}\rangle \langle {{a}_{j}}\rangle$ was used [7]. The bottom panels of figure 4 show the numerically-calculated phase diagram (top row) and the compressibility ($\kappa =\partial \rho /\partial \mu$, ρ  =  filling factor) (bottom row) as functions of the hopping strength for a 1D chain.

Zoom In Zoom Out Reset image size

In both figures the corresponding Mott and SF phases are explored by changing the hopping strength and the effective photonic chemical potential (the free energy term in the JCH Hamiltonian), following the standard treatment of quantum phase transitions. We note that in CRA experiments the hopping strength is in general not tunable but fixed by the fabrication process. In these systems, it is easier to control the effective repulsion (the photon blockade), as discussed in the previous section. Nevertheless, the phase diagrams are useful in showing the existence of Mott lobes in this photonic system and how they depend on system parameters. These works have been supplemented using the variational cluster approach, where on top of the phase diagram, single-particle spectra and finite temperature effects have been investigated [35]. Quantum Monte-Carlo method has also been employed to find various properties of the JCH-type systems [36, 76–78].

An application of the Mott regime for implementing quantum information processing (QIP): The typical implementation of cluster state quantum computing requires initialising all qubits in a 2D lattice in the $|+\rangle =(|0\rangle +\,|1\rangle )/\sqrt{2}$ state and then performing controlled-phase gates (CP) between all nearest-neighbours [79]. In the resonator array system the first proposal to create cluster states utilised the Mott state and the resulting XY dynamics [16]. The ground state $|g,0\rangle$ and the first excited state $|1-\rangle$ of the combined atom-photon (polaritonic) system in each site was used as a qubit as depicted in figure 5.

Zoom In Zoom Out Reset image size

The need to disable the evolution was realised by combining the system's natural dynamics with a protocol where some of the available physical qubits are used as gate 'mediators' and the rest as the logical qubits. The mediators can be turned on and off on demand by tuning the atoms in the corresponding sites on and off resonance from their cavities. This is done by Stark shifting them through the application of an external field (gates A, B, C D in figure 5), effectively shifting the site frequencies and inhibiting photon hopping, thereby isolating each logical qubit. Soon after there have also been other approaches focusing on the generation of the Ising interaction [80] and also towards the simulation of a family of effective spin models in CRAs [15, 41, 81–84]. These will be analysed in the next section.

Quantum spin models have played a crucial role as basic models accounting for magnetic and thermodynamic properties of many-body systems. In implementations involving arrays of Josephson junctions [85] or quantum dots [86], the spin-chain Hamiltonian naturally emerges from the spin-like coupling between qubits, albeit with limited control over the coupling constants. In optical lattice simulators, perturbative dynamics with respect to the Mott-insulator state allows for an effective spin-chain Hamiltonian [72, 87]. This direction has the advantage that the spin-coupling constants can be optically controlled to a great extent. In CRAs, a spin is represented by either polariton states or hyperfine ground levels of the embedded emitter and the system was shown to generate an effective Heisenberg model with optically tunable interactions

$$\begin{eqnarray}&&{{H}_{\text{eff}}}=-\underset{\langle \,j,k\rangle}{\sum}\,{{\lambda}_{z}}\sigma _{j}^{Z}\sigma _{k}^{Z}+{{\lambda}_{x}}\sigma _{j}^{X}\sigma _{k}^{X}+{{\lambda}_{y}}\sigma _{j}^{y}\sigma _{k}^{y}+h\underset{j}{\sum}\,\sigma _{j}^{Z},\end{eqnarray} \tag{ 6 }$$

where the angled brackets denote the sum over nearest neighbours. We will use the same notation throughout this review.

In the work [81], an array of cavities was assumed on the vertices of a 2D lattice (note that any array is in principle possible) as depicted in figure 6(a). Each resonator is doped with a single system (which we refer to as an atom), whose energy level structure comprises of a ground state $|g\rangle$ and two degenerate excited states $|A\rangle$ and $|B\rangle$. Within each lattice site, the Hamiltonian takes the form

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{0}}={{\omega}_{0}}\left(a_{i}^{\dagger}{{a}_{i}}+b_{i}^{\dagger}{{b}_{i}}+|A\rangle \langle A{{|}_{i}}+|B\rangle \langle B{{|}_{i}}\right) \\ +g\left(|A\rangle \langle g{{|}_{i}}\otimes {{a}_{i}}+\,|g\rangle \langle A{{|}_{i}}\otimes a_{i}^{\dagger}\right) \\ +g\left(|B\rangle \langle g{{|}_{i}}\otimes {{b}_{i}}+\,|g\rangle \langle B{{|}_{i}}\otimes b_{i}^{\dagger}\right), \end{array}\end{eqnarray}$$

where $a_{i}^{\dagger}$ and $b_{i}^{\dagger}$ create photons of orthogonal polarisations at site i and the atomic levels and the resonator modes are assumed to be on resonance. The strength of the coupling between the resonator and the atom is denoted by g. Coupling between nearest-neighbour resonators is described by the Hamiltonian

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{\text{hop}}}=\underset{\langle i,j\rangle}{\sum}\,{{J}_{a}}\left(a_{i}^{\dagger}{{a}_{j}}+\text{H}.\text{c}.\right)+{{J}_{b}}\left(b_{i}^{\dagger}{{b}_{j}}+\text{H}.\text{c}.\right). \end{array}\end{eqnarray} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

Assuming strong atom-photon interaction regime, $g\gg {{J}_{a,b}}$, and restricting to the unit filling fraction, such that one excitation per site is expected, one obtains an effective Hamiltonian (6) with

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} h\,=\,\frac{5}{8g}\left(J_{a}^{2}-J_{b}^{2}\right),\;\;\;\;{{\lambda}_{z}}=\frac{9}{32g}\left(J_{a}^{2}+J_{b}^{2}\right), \\ {{\lambda}_{x}}\,=\,{{\lambda}_{y}}=\frac{9{{J}_{a}}{{J}_{b}}}{16g}. \end{array}\end{eqnarray}$$

Note that we have ignored the constant term $\kappa =\frac{31}{32g}\left(J_{a}^{2}+J_{b}^{2}\right)$. The local magnetic fields can be manipulated by applying local Stark fields, thereby leaving an XXZ Hamiltonian where the coefficients ${{\lambda}_{z}}$ and ${{\lambda}_{x}}$ are independently tunable (at device-manufacture stage). This approach allows a stronger spin–spin coupling than an alternative one in [15], but lacks the optical control of the coupling. The latter exhibits greater optical controllability, but relies on rapid switching of optical pulses and the consequent Trotter expansion, which unavoidably introduces additional errors and makes the implementation more difficult.

The question of simulating chains of higher spins, which may have a completely different phase diagram, was dealt with in [88]. In this proposal, a fixed number of atoms are confined in each resonator with collectively applied constant laser fields, and is in a regime where both atomic and resonator excitations are suppressed (figure 6(b)). It was thus shown that as well as optically controlling the effective spin-chain Hamiltonian, it is also possible to engineer the magnitude of the spin.

Since its first discovery in a semiconductor heterostructure in the early 1980s [89], systems exhibiting FQHE has now become routinely available in laboratories. In a different context, the achievement of trapping ultracold atomic gases in a strongly correlated regime has prompted an interest in mimicking various condensed matter systems with trapped atoms, thereby allowing one to tackle such complex systems in unprecedented ways [90]. The great advantage of this approach is that such a tailored system is highly controllable at the microscopic level and is close to idealised quantum models used in condensed matter systems. In the context of FQHE, a 2D atomic gas can be either trapped in a harmonic potential [91, 92] or in an optical lattice [93, 94]. Since the atoms in consideration have no real charge, the magnetic field is simulated artificially. For instance, it can be done by rapidly rotating the harmonic trap [91], by exploiting electromagnetically induced transparency [92], or by modulating the optical lattice potential [93, 94].

The first proposal to realise the fractional quantum Hall system with strongly correlated photons was based on an optical manipulation of atomic internal states and inter-resonator hopping of virtually excited photons in a 2D array of coupled cavities [18]. It was shown that the local addressability of coupled resonator systems allows one to simulate any system of hard-core bosons on a lattice in the presence of an arbitrary Abelian vector potential. The proposal is based on the realisation of the Heisenberg spin model introduced earlier and is depicted in figure 7. To facilitate a spatially varying gauge potential, an asymmetry in the resonator modes is assumed, where two orthogonal resonator modes along the x and y directions have different resonant frequencies (this asymmetry is not an essential requirement as explained in [18]). Realising this geometry would be viable in several promising platforms, such as photonic bandgap microcavities and superconducting microwave cavities. The frequency difference between the two modes is assumed to be much larger than the atom-resonator coupling rates. In this way, either direction of the spin exchange can be accessed individually.

Zoom In Zoom Out Reset image size

Similar works with different approaches and platforms within CRA have appeared since. Ways to introduce synthetic gauge fields—allowing one to break the time-reversal symmetry which is crucial in quantum Hall physics—have been devised in circuit QED [95, 96] and solid state cavities [97]. A photonic version of the quantum spin Hall system, which does not break the time-reversal symmetry, has also been proposed [98] and experimentally demonstrated [99]. Interacting versions realising fractional quantum Hall-type physics of photons were proposed in Kerr-nonlinear cavities [100] and the JCH cavities [101] and signatures of fractional quantum Hall physics in a nonequilibrium driven-dissipative scenario have been investigated [102]. Experimentally, a building block consisting of three superconducting qubits has been demonstrated recently, where both the signatures of a synthetic magnetic field and strong interaction have been observed [103]. Synthetic Landau levels were demonstrated for photons in a single twisted optical cavity also, opening up a new route towards photonic fractional quantum Hall physics [104].

Here we will take a closer look at the first proposal. But first, to set the stage, let us very briefly describe the model Hamiltonian and its ground state wave function that captures the physics of FQHE in its fermionic counterpart. A system of bosonic particles confined in a 2D square lattice in the presence of an artificial perpendicular magnetic field B is described by the Hamiltonian

$$\begin{eqnarray}&&{{H}_{0}}=-t\underset{p,q}{\sum}\,\left(c_{p+1,q}^{\dagger}{{c}_{p,q}}{{\text{e}}^{-\text{i}2\pi \alpha q}}+c_{p,q+1}^{\dagger}{{c}_{p,q}}+\text{H}.\text{c}.\right),\end{eqnarray} \tag{ 8 }$$

where the positions of lattice sites are represented by $a\left(\,p\hat{x}+q\hat{y}\right)$ with p and q being integers and α is the spacing of the lattice. Here, $\alpha \equiv B{{a}^{2}}/{{ \Phi }_{0}}$, the number of magnetic flux quanta through a lattice cell, plays a crucial role in characterising the energy spectrum, whose self-similar structure is known as the Hofstadter butterfly. In the hardcore (maximum one boson per site) and continuum limit ($\alpha \ll 1$) the ground state is accurately described by the Laughlin state

$$\begin{eqnarray}&&{{ \Psi }_{m}}\left(\left\{{{z}_{j}}\right\}\right)={{\text{e}}^{-\frac{1}{4}\underset{j}{\sum}\,|{{z}_{j}}{{|}^{2}}}}\underset{j<k}{\prod}\,{{\left({{z}_{j}}-{{z}_{k}}\right)}^{m}},\end{eqnarray} \tag{ 9 }$$

where ${{z}_{j}}={{x}_{j}}+\text{i}{{y}_{j}}$ and the filling factor $\nu =1/m$, with m even, corresponds to the ratio of the number of bosons to the number of magnetic flux quanta.

Assume now a 2D array of coupled cavities, each confining a single atom with two ground levels as in figure 7(b). Each resonator has two modes along the x and y directions, whose annihilation operators will be denoted by a^x and a^y. The atom interacts with these resonator modes with coupling rates g^x and g^y, and with detunings ${{ \Delta }^{x}}$ and ${{ \Delta }^{y}}$, respectively. Furthermore, to induce site-dependent hopping, two classical fields with (complex) Rabi frequencies ${{ \Omega }^{x}}{{\text{e}}^{-\text{i}\theta _{(\,p,q)}^{x}}}$ and ${{ \Omega }^{y}}{{\text{e}}^{-\text{i}\theta _{(\,p,q)}^{y}}}$ are applied. In the rotating frame, the Hamiltonian for the system reads

$$\begin{eqnarray}\begin{array}{*{35}{l}} H = & \underset{\begin{array}{c} \mu =x,y \\ j=(\,p,q) \end{array}}{\sum}\,\left[{{g}^{\mu}}{{\text{e}}^{-\text{i}{{ \Delta }^{\mu}}t}}a_{j}^{\mu}(|e\rangle {{\langle 0|)}_{j}} \right. \\ {} & \left. +\,{{ \Omega }^{\mu}}{{\text{e}}^{-\text{i}\theta _{j}^{\mu}}}{{\text{e}}^{-\text{i}{{ \Delta }^{\mu}}t}}(|e\rangle {{\langle 1|)}_{j}}+\text{H}.\text{c}.\right] \\& -\underset{p,q}{{\,{\sum}^{}}}\,\left({{J}^{x}}a_{p+1,q}^{x\dagger}a_{p,q}^{x}+{{J}^{y}}a_{p,q+1}^{y\dagger}a_{p,q}^{y}+\text{H}.\text{c}.\right), \end{array}\end{eqnarray} \tag{ 10 }$$

where J^x (J^y) denotes the inter-resonator hopping rate of the photon along the x (y) direction. We assume ${{ \Delta }^{x}}-{{ \Delta }^{y}}\gg {{g}^{x}},{{g}^{y}}$, as mentioned above, and furthermore that ${{ \Delta }^{\mu}}\gg {{g}^{\mu}}\gg {{ \Omega }^{\mu}},{{J}^{\mu}}$.

In this strong coupling regime, the atomic excitation is suppressed and one can perform adiabatic elimination to obtain an effective Hamiltonian

$$\begin{eqnarray}\begin{array}{*{35}{l}} H = & \underset{\begin{array}{c} \mu =x,y \\ j=(\,p,q) \end{array}}{\sum}\,\left[{{\delta}^{\mu}}a_{j}^{\mu \dagger}a_{j}^{\mu}(|0\rangle {{\langle 0|)}_{j}}+{{\omega}^{\mu}}\left({{\text{e}}^{\text{i}\theta _{j}^{\mu}}}a_{j}^{\mu}\sigma _{j}^{+}+\text{H}.\text{c}.\right)\right] \\& -\underset{p,q}{\sum}\,\left({{J}^{x}}a_{p+1,q}^{x\dagger}a_{p,q}^{x}+{{J}^{y}}a_{p,q+1}^{y\dagger}a_{p,q}^{y}+\text{H}.\text{c}.\right), \end{array}\end{eqnarray} \tag{ 11 }$$

where ${{\delta}^{\mu}}={{\left({{g}^{\mu}}\right)}^{2}}/{{ \Delta }^{\mu}}$, ${{\omega}^{\mu}}={{g}^{\mu}}{{ \Omega }^{\mu}}/{{ \Delta }^{\mu}}$, and ${{\sigma}^{+}}=|1\rangle \langle 0|$. Here, we have ignored the ac Stark shift induced by the classical fields, which is negligible in the assumed regime (it may also be compensated by other lasers). We further assume ${{\delta}^{\mu}}\gg {{J}^{\mu}}\gg {{\omega}^{\mu}}$, which can be satisfied, along with the above conditions, when

$$\begin{eqnarray}&&{{g}^{\mu}}/{{ \Delta }^{\mu}}\gg {{J}^{\mu}}/{{g}^{\mu}}\gg {{ \Omega }^{\mu}}/{{ \Delta }^{\mu}}.\end{eqnarray} \tag{ 12 }$$

Now the resonator photon is suppressed and an adiabatic elimination can be applied once more. Keeping the terms up to the third order and taking only the sub-manifold with no photons, the effective Hamiltonian further simplifies to

$$\begin{eqnarray}\begin{array}{*{35}{l}} H = & -t\underset{p,q}{\sum}\,\left[\sigma _{p+1,q}^{+}\sigma _{p,q}^{-}{{\text{e}}^{\text{i}\left(\theta _{p+1,q}^{x}-\theta _{p,q}^{x}\right)}} \right. \\ &\left. \,+\,\sigma _{p,q+1}^{+}\sigma _{p,q}^{-}{{\text{e}}^{\text{i}\left(\theta _{p,q+1}^{y}-\theta _{p,q}^{y}\right)}}+\text{H}.\text{c}.\right], \end{array}\end{eqnarray} \tag{ 13 }$$

where ${{\delta}^{\mu}}$, ${{\epsilon}^{\mu}}$, ${{\omega}^{\mu}}$, and ${{J}^{\mu}}$ are chosen to yield

$$\begin{eqnarray}&&t={{J}^{x}}{{\left(\frac{{{\omega}^{x}}}{{{\delta}^{x}}}\right)}^{2}}={{J}^{y}}{{\left(\frac{{{\omega}^{y}}}{{{\delta}^{y}}}\right)}^{2}}.\end{eqnarray} \tag{ 14 }$$

It is easy to see that this Hamiltonian reduces to the Hamiltonian H₀ in (8) in the hardcore bosonic limit, if we adjust the phases of the classical fields such that $\theta _{p,q}^{x}=-pq\centerdot 2\pi \alpha$ and $\theta _{p,q}^{y}=0$.

The ground state of the Hamiltonian (13) is thus expected to be described faithfully by the Laughlin wavefunction (9), modulo changes due to the difference in gauge and boundary conditions, which can be calculated analytically [105]. For a 4 by 4 lattice with $\alpha =1/4$, m  =  2, the fidelity between the numerically calculated ground state of the ideal Hamiltonian and this modified Laughlin state (with m  =  2) was shown to be excellent indeed ($\approx 0.989$) [18]. Actually this number is a bit misleading because $\alpha =1/4$ is quite a large value for which the continuum limit is not guaranteed to be a good approximation. Indeed, for 4 bosons, the bosonic Laughlin state ceases to be a good variational ground state for $\alpha >0.2$ [106]. In such a case, one can compute topological invariants such as the Chern number [106], or look at the edge properties [99], in order to verify the topological nature of the ground state.

The numerical ground state of the effective Hamiltonian (11) for ${{\delta}^{\nu}}/10=10{{\omega}^{\mu}}={{J}^{\mu}}$, corresponding to ${{ \Delta }^{\mu}}/1000=$ ${{g}^{\mu}}/100={{ \Omega }^{\mu}}={{J}^{\mu}}$, has the fidelity of 0.976 [18]. The ground state can be prepared adiabatically by exciting two atoms at different sites initially without the laser fields (t  =  0) and then delocalising these excitations via adiabatically tuning up the effective hopping rate. The quasi-excitations can also be generated via local tuning of the laser fields to mimic an infinitely thin solenoid, around which the excitations form. More details on these issues can be found in [18, 107]. Lastly, it should be mentioned that one needs a higher number of bosons in order to perform braiding with anyonic excitations [108, 109].

Let us begin by describing how to model the driven-dissipative systems. Firstly, coherent driving by laser fields can be modelled by introducing a driving term in the Hamiltonian, ${{H}_{\text{drive}}}={\sum}_{i}\left[{{ \Omega }_{i}}(t)a_{i}^{\dagger}+ \Omega _{i}^{\ast}(t){{a}_{i}}\right]$, where ${{ \Omega }_{i}}(t)={{ \Omega }_{i}}\exp \left(-\text{i}{{\omega}_{L}}t\right)$ is the time-dependent driving field and a_i is the driven resonator mode. As we have noted earlier, we will encounter two types of model Hamiltonians. One is the Bose–Hubbard Hamiltonian which, including the driving term, can be written as

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{H}_{\text{BH}}} = & \underset{i}{{{\sum}^{}\,}}\,\left[{{\omega}_{c}}a_{i}^{\dagger}{{a}_{i}}+\frac{U}{2}a_{i}^{\dagger}a_{i}^{\dagger}{{a}_{i}}{{a}_{i}}\right]\,-\underset{\langle i,j\rangle}{\sum}\,Ja_{i}^{\dagger}{{a}_{j}} \\& +\underset{i}{\sum}\,\left[{{ \Omega }_{i}}(t)a_{i}^{\dagger}+ \Omega _{i}^{\ast}(t){{a}_{i}}\right], \end{array}\end{eqnarray} \tag{ 15 }$$

while the other is the JCH Hamiltonian

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{H}_{\text{JCH}}} = & \underset{i}{{{\sum}^{}\,}}\,\left[{{\omega}_{c}}a_{i}^{\dagger}{{a}_{i}}+\left({{\omega}_{c}}- \Delta \right)\sigma _{i}^{+}\sigma _{i}^{-} \right. \\& \left. +g\left(a_{i}^{\dagger}\sigma _{i}^{-}+{{a}_{i}}\sigma _{i}^{+}\right)\right] \\& -\underset{\langle i,j\rangle}{\sum}\,Ja_{i}^{\dagger}{{a}_{j}}+\underset{i}{\sum}\,\left[{{ \Omega }_{i}}(t)a_{i}^{\dagger}+ \Omega _{i}^{\ast}(t){{a}_{i}}\right]. \end{array}\end{eqnarray} \tag{ 16 }$$

In both models, ${{\omega}_{c}}$ is the bare resonator frequency and J is the photon hopping rate between the nearest neighbour resonators—denoted by the angled brackets. $\Delta$ is the detuning between the resonator and atom and g is the atom-photon coupling strength. $\sigma _{i}^{+}$ and $\sigma _{i}^{-}$ are the atomic raising and lowering operators at site i.

Secondly, photon losses are modelled by assuming that the system is (weakly) interacting with a bath of harmonic oscillators. Taking the Born–Markov approximation, one obtains a quantum master equation describing the dynamics of the reduced density matrix of the system ρ [123]. In this review, we will mostly deal with master equations of the form

$$\begin{eqnarray}\begin{array}{*{35}{l}} \overset{\centerdot}{{\rho}}\, & = & -\text{i}\left[{{H}_{X}},\rho \right]+\frac{\gamma}{2}\underset{i=1}{\overset{N}{\sum}}\,\left(2{{a}_{i}}\rho a_{i}^{\dagger}-a_{i}^{\dagger}{{a}_{i}}\rho -\rho a_{i}^{\dagger}{{a}_{i}}\right) \\ {} & {} & +\,\,\frac{\kappa}{2}\underset{i=1}{\overset{N}{\sum}}\,\left(2\sigma _{i}^{-}\rho \sigma _{i}^{\dagger}-\sigma _{i}^{\dagger}{{\sigma}_{i}}\rho -\rho \sigma _{i}^{\dagger}{{\sigma}_{i}}\right), \end{array}\end{eqnarray} \tag{ 17 }$$

where X can be $\text{JCH}$ or $\text{BH}$ (in the latter case, $\kappa =0$ is assumed); γ and κ are the dissipation rates of photons and atoms, respectively, assumed here to be identical for all sites.

Now we are ready to compare the two models. Let us start from the simplest case: a single-site system. The eigenstates of the Jaynes–Cummings Hamiltonian are the 'dressed' states $|n,\pm \rangle$, where n is the total excitation number and  ±  labels two polaritonic (hybrid atom-photon) modes with energies $\omega _{n}^{\pm }=n{{\omega}_{c}}- \Delta /2\pm \sqrt{{{(\Delta /2)}^{2}}+n{{g}^{2}}}$. Because of the dependence on n, the energy of $|2,-\rangle$ is not twice that of $|1,-\rangle$, as illustrated in figure 8(a) for $\Delta =0$. Comparing with the energy levels of the Kerr-nonlinear model shown in figure 8(b), one thus concludes that if '−' polaritons are driven and if the driving is sufficiently weak such that only the lowest states are occupied, an effective interaction strength U_eff can be defined to formally connect the two models. That is, one can assign an effective BH parameter

$$\begin{eqnarray}&&\frac{{{U}_{\text{eff}}}}{g}=\frac{\omega _{2}^{-}-2\omega _{1}^{-}}{g}=\frac{ \Delta }{2g}+2\sqrt{{{\left(\frac{ \Delta }{2g}\right)}^{2}}+1}-\sqrt{{{\left(\frac{ \Delta }{2g}\right)}^{2}}+2},\end{eqnarray} \tag{ 18 }$$

to the JCH model.

Zoom In Zoom Out Reset image size

Deep in the strongly interacting regime, we would thus expect the two models to yield similar qualitative behaviour: the polaritons will delocalise to minimise the kinetic energy, but the double occupation of any site is inhibited at the same time. On the other hand, we also know that the two models should converge in the opposite limit of linear CRAs. These two observations give a heuristic argument as to why we can see the Mott-superfluid transition in a JCH system. However, this does not mean that the two types of systems are similar in all aspects. After all, there are additional atomic degrees of freedom in the JCH systems. To illustrate similarities and differences in the two models, figure 9 plots various expectation values of the nonequilibrium steady state (NESS) for a homogeneously-driven periodic three-site CRA.

Zoom In Zoom Out Reset image size

Upon comparing the JCH (figure 9(a)) and the BH (figure 9(b)) results, we immediately see that the JCH has a more complex structure. Broadly, the spectrum is divided into two parts. The left (right) wing is mostly composed of the '−' ('+') polaritons for the small value of J used here (U  =  20J). The vertical dash-dotted line indicates the driven single-excitation eigenstate while the dotted lines indicate the eigenstates of a single resonator (J  =  0 case). A BH system by comparison has a relatively simple spectrum with peaks near (but not exactly at) the single-site n-photon energies indicated by the dotted vertical lines. Note that the spectrum is quite broad as we have assumed $\gamma =J$. Despite these differences, there is a notable similarity when the systems are driven at the single-particle resonance (vertical dash-dotted lines)—enhanced steady-state photon number is accompanied by a dip in the same-time intensity–intensity correlation function ${{g}^{(2)}}={{g}^{(2)}}\left({{i}_{1}},{{i}_{1}}\right)$. The latter is defined as

$$\begin{eqnarray}&&{{g}^{(2)}}\left({{i}_{1}},{{i}_{2}}\right)=\langle a_{{{i}_{1}}}^{\dagger}a_{{{i}_{2}}}^{\dagger}{{a}_{{{i}_{2}}}}{{a}_{{{i}_{1}}}}\rangle /\langle a_{{{i}_{1}}}^{\dagger}{{a}_{{{i}_{1}}}}\rangle \langle a_{{{i}_{2}}}^{\dagger}{{a}_{{{i}_{2}}}}\rangle,\end{eqnarray} \tag{ 19 }$$

where the expectation values are taken with respect to the steady state density operator. The dip in g⁽²⁾ signifies that two excitations do not want to occupy the same site simultaneously. This quantifies the photon blockade effect mentioned earlier, arising due to the strong effective repulsion between photons.

At this point, the reader might be wondering whether it is possible at all to achieve a good quantitative agreement between the two models. It is indeed possible in the 'photonic' regime, in which the atom-photon detuning is so large that the polaritons are mostly photonic. In this limit, however, the effect of atom-photon interaction is to merely shift the resonance (figure 9(c)). To counteract the reduced nonlinear effect, one needs a very large value of g. In figure 9(d), $\Delta /g=-10$ and $g=6\times {{10}^{4}}\gamma$ were used. Since these kind of numbers lie well outside the experimental realm, we conclude that differences between the two types of systems are there in realistic systems and one must be careful in choosing the correct model to describe a given CRA. In particular, one should check whether a phenomenon found in a BH system can also be observed in a JCH system in a reasonable parameter regime. We will come across such examples below.

Initial investigations on many-body physics in driven-dissipative CRAs employed the 1D BH model. Signatures of strong on-site interaction leading to fermionisation and crystallisation of the excitations (much as in Tonks–Girardeau gases) were shown to be visible in the NESS of driven-dissipative CRAs systems [124, 125]. JCH versions of these models have been studied by Grujic et al [122], confirming that the predicted phenomena also occur in the JCH versions. The effect of geometric frustration on the steady-state has also been studied more recently, where it was shown that in a homogeneously-driven flat-band lattice, the steady state exhibits an incompressible state of photons characterised by a short-ranged crystalline order [126].

3.2.1. Fermionisation of photons.

We begin with the work of Carusotto and coworkers [124]. The system is composed of Kerr nonlinear cavities coupled in a cyclic fashion as illustrated in figure 10. The master equation describing the system is given by (17) with $X=\text{BH}$, with homogeneous driving fields such that ${{ \Omega }_{i}}= \Omega$. As shown on the right panel of figure 10, the 'transmission' spectrum of the total 'transmitted' intensity ${{n}_{T}}={\sum}_{i}\text{Tr}\left[\rho a_{i}^{\dagger}{{a}_{i}}\right]\equiv \langle a_{i}^{\dagger}{{a}_{i}}\rangle$ show the traces of strong interaction between photons. For the case of impenetrable bosons ($U\gg J$), the Bose–Fermi mapping can be employed [127], which translates an interacting Boson model into a free Fermion model to enable calculations of the energy eigenvalues. These eigenvalues with their eigenstates are shown in the figure, showing excellent agreement between the spectra of the driven-dissipative system and the eigenvalues of the closed system.

Zoom In Zoom Out Reset image size

To further illustrate the fermionisation of photons, figure 11 plots the on-site (${{g}^{(2)}}\left({{i}_{1}}={{i}_{2}}\right)$) and nearest-neighbour (${{g}^{(2)}}\left({{i}_{1}}\ne {{i}_{2}}\right)$) correlations as functions of U/J, when the pump lasers are kept on resonance with the lowest two-photon transition. Clearly, both the JCH and BH systems yield qualitatively similar results when the effective nonlinearity is matched via equation (18). In the weak interaction limit, all the levels are harmonic and the emitted light inherits the Poissonian nature of the pump. In the intermediate regime, $U/J\approx 1$, the (driven) two photon peak starts to separate from the single photon peak leading to bunching. Finally, in the large U/J limit, the two photons fermionise, resulting in strong antibunching with enhanced cross-correlations.

Zoom In Zoom Out Reset image size

A detailed investigation on the fluorescence spectrum and the spectrum of the second-order correlation function under coherent and incoherent driving for a two coupled cavities of JCH type was carried out by Knap et al [128]. Similar to the work described above, the resonance peaks of the measured spectra of the dissipative system were matched to the eigenenergies of the closed system.

3.2.2. Crystallisation of polaritons.

In an independent work, Hartmann has studied the effects of strong interaction between polaritons in coupled resonators (described by the BH model in which a_j is the polariton annihilation operator) [125]. In that setup, instead of driving the system homogeneously, a flow of polaritons is created via modulation of the phase of the driving laser such that ${{ \Omega }_{j}}= \Omega {{\text{e}}^{-\text{i}{{\phi}_{j}}}}$, where ${{\phi}_{j}}=j\pi /2$. The interplay between this flow and the nonlinearity then results in polariton crystallisation, where excitations are predominantly found within a specific distance from one another. To describe this phenomenon in detail, we will use the spatial correlation function g⁽²⁾(i, j) as well the momentum correlation function defined by

$$\begin{eqnarray}&&g_{\text{BM}}^{(2)}(\,p,q)=\langle B_{p}^{\dagger}B_{q}^{\dagger}{{B}_{q}}{{B}_{p}}\rangle /\langle B_{p}^{\dagger}{{B}_{p}}\rangle \langle B_{q}^{\dagger}{{B}_{q}}\rangle,\end{eqnarray} \tag{ 20 }$$

where the momentum eigenmode operator is defined as ${{B}_{p}}=1/\sqrt{N}{\sum}_{j}\exp \left(\text{i}pj\right){{a}_{j}}$, with N the total number of cavities and $p=2\pi l/N$ for l  =  −N/1  +  1, −  N/2  +  2,..., N/2.

There are two interesting regimes determined by the relative strength between the nonlinearity U and the driving strength $\Omega$: (1) Weakly nonlinear regime, $\Omega \gg U$, and (2) Strongly nonlinear regime, $U\gg \Omega$. Let us start with the weakly nonlinear regime. In this strongly-driven regime, a semiclassical treatment is valid, i.e. ${{B}_{p}}=\sqrt{N}\beta {{\delta}_{p,\pi /2}}+{{b}_{p}}$, where β is a mean amplitude for the driven mode at $p=\pi /2$ and b_p is an annihilation operator for the quantum of fluctuations. To second order in fluctuations the Hamiltonian, neglecting the driving term for now, can be written as

$$\begin{eqnarray}&&\tilde{H}=\underset{p}{\sum}\,\left[\left({{\omega}_{p}}+2Un\right)b_{p}^{\dagger}{{b}_{p}}+\left(\frac{U{{\beta}^{2}}}{2}b_{p}^{\dagger}b_{{\bar{p}}}^{\dagger}+\text{H}.\text{c}.\right)\right],\end{eqnarray} \tag{ 21 }$$

where $\bar{p}=\pi p/|\,p|-p$ and $n=|\beta {{|}^{2}}$ is the mean number of polaritons. This Hamiltonian is known to exhibit two-mode squeezing between the modes p and $\bar{p}$, arising due to the scattering between p and $\bar{p}$ modes induced by the nonlinear term. The latter also induces density–density correlations between the modes, which can be seen from the intensity correlation function between the Bloch modes evaluated to the leading order in 1/(Nn):

$$\begin{eqnarray}\begin{array}{*{35}{l}} g_{\text{BM}}^{(2)}(\,p,q) \approx & 1+{{\delta}_{p\ne \pi /2}}\left({{\delta}_{p,q}}+{{\delta}_{p+q,\pm \pi}}\frac{|g{{|}^{2}}}{{{m}^{2}}}\right) \\& -{{\delta}_{p,\pi /2}}{{\delta}_{q,\pi /2}}\frac{2m}{Nn}, \end{array}\end{eqnarray} \tag{ 22 }$$

where $m=2{{U}^{2}}{{n}^{2}}/\left(12{{U}^{2}}{{n}^{2}}+{{\gamma}^{2}}\right)$ is the contribution to the mean photon number from quantum fluctuations and $g=\langle {{B}_{p}}{{B}_{{\bar{p}}}}\rangle =-{{\beta}^{2}}\left(4{{U}^{2}}n+\text{i}U\gamma \right)/\left(12{{U}^{2}}{{n}^{2}}+{{\gamma}^{2}}\right)$. The term ${{\delta}_{p+q,\pm \pi}}\frac{|g{{|}^{2}}}{{{m}^{2}}}$ clearly shows the correlations between the modes p and $\bar{p}$ for non-driven modes, while the driven mode $p=\pi /2$ exhibits slight antibunching $g_{m}^{(2)}(\pi /2,\pi /2)\approx 1-2m/Nn$.

In the strongly nonlinear regime, the semiclassical description is no longer valid. Instead, the time evolving block decimation (TEBD) algorithm (see [129, 130] for reviews) was used to calculate the steady state for an array of 16 cavities (for other methods to deal with larger systems, see e.g. [131, 132]). Contrary to the weakly nonlinear regime, there is now strong correlations between the same Bloch modes ($g_{m}^{(2)}(\,p,p)>1$), while different modes are strongly anticorrelated ($g_{m}^{(2)}(\,p,q)<1$, $p\ne q$), as shown in figure 12.

Zoom In Zoom Out Reset image size

The same figure also shows an evidence of photon crystallisation, by which we mean the following: if a polariton is found in one resonator, a second polariton is most likely to be found in an adjacent resonator and this happens at the reduced probability of finding a polariton at farther sites. This is clearly observed in the polariton correlation function between the cavities, where g⁽²⁾(j, j  +  1)  >  1 and g⁽²⁾(j, j  +  k)  <  1 for k  =  1, 2,.... The correlations fall below 1 only when the driving phase is $\ne 0$, indicating that the observed phenomenon emerges due to the interplay between the nonlinearity and the flow of polaritons. For comparison, figure 12 also plots (green lines) the density correlations of a Tonks–Girardeau (TG) gas $g_{\text{TG}}^{(2)}(i,j)=1-{{\left(\sin \left(\pi \tilde{n}(i-j)\right)/\pi \tilde{n}(i-j)\right)}^{2}}$, where $\tilde{n}$ is the number of polaritons per resonator. While the TG gas shows an oscillation in the density–density correlation, no such oscillation is observed in crystallised polaritons. Similar behaviour is observed for the JCH model when the polaritonic mode $|1,-{{\rangle}_{p=\pi /2}}$ is driven, as depicted in figure 13. It is interesting to note that the on-site repulsion is strengthened with increasing $\Delta /g$—this can be attributed to an increase in the atomic portion of the JC polariton with the detuning.

Zoom In Zoom Out Reset image size

In the previous section we have seen that the steady states of driven-dissipative CRAs are capable of exhibiting signatures of underlying many-body states in the equilibrium system and furthermore that interesting phases can be created due to driving. However, this is not the only nonequilibrium scenario that reveals interesting physics. One well-known example is quenching, in which the system is initially excited in a particular initial state and subsequent dynamics is observed. Here we review two early works taking this approach. We begin with a study of the nonequilibrium signatures of the superfluid-insulator transition by Tomadin et al [133] and briefly summarise a similar study on a hybrid system by Hummer et al [49]. Then we review the 'localisation–delocalisation' transition of photons [134], which has been verified experimentally in a superconducting circuit setup [135].

3.3.1. Superfluid-insulator transition.

In [133], Tomadin and coworkers investigated a dissipative BH system in 2 dimensions without any driving field, initially prepared in a Mott-like state. Much as in the quantum quench scenario, subsequent transient evolutions depend sensitively on the underlying equilibrium phase (without dissipation). A self-consistent cluster mean-field approach was employed with up to two sites per cluster, where the evolution within a cluster is exact but the interaction with the rest of the system is treated at the mean-field level.

Let us first consider the dynamics of the initial Mott state. With negligible dissipation, the dynamics obviously depends sensitively on the value of J/U. For small J/U, the system is expected to stay more or less in the same state, whereas for large J/U, nontrivial dynamics are expected. As in the equilibrium case, a useful observable to distinguish between these distinct cases is the second-order correlation function, since strong antibunching is expected to persist for small J/U. In fact, in the transient scenario, a more useful observable was shown to be the time-averaged second-order correlation function, which we denote ${{\langle g\rangle}_{t}}$. Figure 14 illustrates the usefulness of this quantity, where the averaging is over $2/\gamma <t<4/\gamma$. The inset shows the time-averaged mean-field order parameter ${{\langle |\psi |\rangle}_{t}}$ displaying similar behaviour as ${{\langle g\rangle}_{t}}$; z is the coordination number denoting the number of nearest neighbours. For small J/U, the average intensity correlation is 0 indicating that the system remains in the initial Mott state. However, for large J/U, the initial state becomes unstable and a crossover towards Poisson statistics (g⁽²⁾  =  1) is observed, the latter occurring over a narrow region of J/U near the critical point indicated by the vertical dashed line.

Zoom In Zoom Out Reset image size

Experimentally, the initial state is most likely to be prepared by a series of laser pulses. In this case, there will inevitably be some imperfections in preparing the Mott state. Tomadin et al have studied the effect of such imperfections by considering non-ideal filling factors in the initial state. Figure 15 shows the results. Clearly, a deviation from the ideal initial state does not destroy the stark difference between the superfluid (triangles) and Mott (circles) regimes as measured by ${{\langle {{g}^{(2)}}\rangle}_{t}}$. The same conclusion holds for a general initial condition, with up to 20% depletion of the average filling, for the whole range of initial coherences allowed by the manifold spanned by 0 and 1 Fock states in each resonator.

Zoom In Zoom Out Reset image size

A similar investigation on nonequilibrium signatures of the phase transition has been carried out by Hummer et al [49], in a hybrid system of flux qubits coupled with nitrogen-vacancy centres. The model describing this system is different from the BH system in that (i) each site comprises of qubits coupled to collective bosonic modes of NV centres surrounding it, and (ii) the coupling is between the qubits, not bosons. Despite these differences, the time averaged probability to find two excitations in a single site is shown to be a very good indicator of the underlying phase transition in the equilibrium case.

3.3.2. Localisation–delocalisation transition of photons.

The transient phenomenon to be reviewed here is not associated with the underlying phases of the equilibrium system, but with the dynamics of the system. The concept was initially introduced by Schmidt et al [134] essentially in a closed Jaynes–Cummings dimer, but the subsequent experimental demonstration [135] showed that dissipation can play an important role in achieving the transition dynamically. As this experiment provides the state-of-the-art realisation of a CRA, we discuss it in some detail. The essential physics of the localisation–delocalisation transition is similar to that of self-trapping: when two nonlinear systems are tunnel-coupled, the oscillation between the systems can be inhibited by the on-site nonlinearity. While previous studies in optical fibres [136], molecules [137], cold atoms [138–140], and polaritons [141] have concentrated on the Kerr nonlinearity, Schmidt et al showed that the same phenomenon also occurs in a Jaynes Cummings dimer. The effect has been generalised to more than two sites by Schetakis and coworkers [142].

As first explained in [134], the physics of the localisation–delocalisation transition are nicely captured in the semiclassical picture without dissipation. One resonator is assumed to be initially in a coherent state with N photons on average, while the other resonator is assumed to be empty. Then, the semiclassical analysis shows that the photon imbalance $Z(t)=\left({{n}_{L}}(t)-{{n}_{R}}(t)\right)/N$ undergoes an oscillation whose frequency decreases with increasing g and furthermore goes to 0 at the critical coupling strength $g_{\text{c}}^{cl}\approx 2.8\sqrt{N}J$. This behaviour survives in the quantum regime as shown in figure 16 for a lossless dimer with 7 photons initially placed in the left resonator. $\bar{Z}(t)$ is the time averaged imbalance up until time t. While the initial imbalance is washed out over time below the critical coupling strength, $g=0.3g_{\text{c}}^{cl}$ (figure 16(a)), the signature of photon localisation is evident above the critical coupling strength, $g=2.0g_{\text{c}}^{cl}$ (figure 16(b)). The localisation behaviour can be explained as follows. For a symmetric dimer, the two states with opposite imbalance 1 and  −1 are approximate eigenstates of the system when $g\gg J$. Using the degenerate perturbation theory the period of oscillation between these states can be calculated as $1/\left({{c}_{{{N}_{0}}}}J{{(J/g)}^{{{N}_{0}}-1}}\right)$ with ${{c}_{{{N}_{0}}}}$ a constant that depends on the total photon number N₀. Clearly, this number diverges rapidly with increasing N₀ when $J/g\ll 1$.

Zoom In Zoom Out Reset image size

A detailed analysis shows that there are features not captured by the semiclassical treatment such as a beating in the coherent oscillation that is observable in the weakly-nonlinear case, smoothening of the transition, and a shift in the critical coupling strength where $g_{\text{c}}^{qu}<g_{\text{c}}^{cl}$ [134]. The shift in the critical coupling strength has been explained recently using a Fock-state WKB approach, yielding an analytical expression for the classical phase boundary depending on the initial state of the spins [143]. Rather than delving into these features further, let us move on to the experimental realisation of the dissipative localisation–delocalisation transition, where dissipation plays a key role in driving the transition dynamically. The latter occurs due to the photon number dependence of the critical coupling strength. For an experimental setup with fixed g, the initial photon number N₀ decides whether the system is in the localised regime ($g>g_{\text{c}}^{qu}\left({{N}_{0}}\right)$) or in the delocalised regime ($g<g_{\text{c}}^{qu}\left({{N}_{0}}\right)$). In a closed system, the system has no choice but to remain in a given regime. However, in a dissipative system, the decrease in the total number of photons can take the system from a delocalised regime to a localised regime (because $g_{\text{c}}^{qu}$ decreases with N).

This effect has been experimentally demonstrated by Raftery and coworkers in a superconducting circuit platform [135]. In the experiment, the initial state was prepared in a coherent state instead of a Fock state, but the described physics remain qualitatively unchanged. Initially, the qubits and the resonators are far-detuned so they are practically uncoupled. Then the system is driven by a coherent pulse and left to evolve a specific period of time, which prepares the system with the imbalance 1. The nonlinearity is then rapidly increased by bringing the qubits into resonance with the resonator mode. The results for various initial photon numbers (N_i) are plotted in figure 17. It shows the homodyne signal in the initially unoccupied resonator as a function of the initial photon number and time. For a small enough N_i, the system starts off and remains in the localised regime, whereas for a larger N_i the system starts off in the delocalised regime and crosses over to the localised regime after enough photons are lost from the system so that the critical coupling strength is smaller than the experimental coupling strength.

Zoom In Zoom Out Reset image size

Coupled resonator arrays, especially in the superconducting circuit platform, allow one to engineer different types of interactions that are difficult to achieve in other systems. Here, we look at two proposals utilising this advantage to realise exotic phases in nonequilibrium scenarios. We begin by reviewing a proposal on a photonic supersolid phase in a driven-dissipative setup and then a proposal to achieve Majorana-like modes in a chain of driven cavities.

3.4.1. Photon supersolid phase.

Let us start with the proposal by Jin and coworkers, which investigated the consequences of adding a cross-Kerr term between nearest neighbours [144, 145]. Such a term naturally arises in electronic systems due to the Coulomb interaction between electrons, but are necessarily much smaller than the on-site interaction term. However, in CRAs the magnitude of the new term can be made arbitrarily large as compared to the on-site interaction term, giving rise to interesting new phases. The open extended Hubbard system studied by the authors is described by the master equation (17) with

$$\begin{eqnarray}&&{{H}_{X}}={{H}_{\text{BH}}}+V\underset{\langle i,j\rangle}{\sum}\,{{n}_{i}}{{n}_{j}},\end{eqnarray} \tag{ 23 }$$

which, as they show, can be realised using nonlinear coupling elements between the cavities.

Consider a 2D bipartite lattice—lattice that divides into two sublattices—with homogeneous coherent pumping (${{ \Omega }_{i}}= \Omega$). Due to the presence of the cross-site interaction, a non-zero photon number difference between the sublattices may develop. This broken symmetry phase is characterised by the order parameter $\Delta n=|\langle {{n}_{A}}\rangle -\langle {{n}_{B}}\rangle |$, where A and B are sublattice indices. Figure 18(a) shows the order parameter as a function of the hopping strength and cross-Kerr nonlinearity for U  =  1. The coloured region indicates the broken-symmetry (or crystalline) phase, where the photon number distribution exhibits a checkerboard pattern. We see that the interaction term favours the broken symmetry phase, whereas the hopping tends to restore the symmetry. Remarkably though, having $V\ne 0$ is not a prerequisite to obtain the checkerboard phase, unlike in the equilibrium counterpart [144]: it is possible, although much harder, for the nonequilibrium steady-state to have non-zero $\Delta n\ne 0$ even with V  =  0. The requirements are a non-zero initial imbalance and a finite hopping strength.

Zoom In Zoom Out Reset image size

In figure 18(b), the effects of driving field detuning, δ, is shown for U  =  0 and zV  =  0.6. Here the situation gets more interesting. A new phase appears (shaded green area) where the steady-state never becomes truly stationary in that there is a residual time dependence of $\langle a\rangle$ (on top of a trivial time dependence due to driving). A closer inspection on the reduced density matrix of the sublattice shows that there is a global dynamical phase coherence on top of the checkerboard pattern, which led Jin and coworkers to call this phase a photonic supersolid phase. Similar phases were also found in the equilibriums setting for an extended JCH model with the nearest neighbour interaction between the qubits [146].

The existence of the crystalline phase was corroborated by performing numerical simulations using the matrix product operator technique in 1D [125]. This approach takes into account quantum fluctuations ignored in the mean-field analysis and is therefore much more reliable. The results are shown in figure 19. The nature of the crystalline phase is displayed in the second-order intensity correlation function g⁽²⁾(i, j). Figure 19(a) plots this function for 20 cavities in the absence of hopping (J  =  0). In this regime, the staggered distribution of photons is manifest, especially for larger values of V. Figure 19(b) plots g⁽²⁾ for 21 cavities for various values of hopping strengths, displaying the suppression of crystalline phase with increasing hopping strength.

Zoom In Zoom Out Reset image size

3.4.2. Majorana-like modes.

As the last detailed review of this section we look at a proposal by Bardyn and Imamoglu on a photonic realisation of Majorana-like modes (MLMs) in a 1D CRA [147]. The nonequilibrium steady state of the driven-dissipative version of the model has been studied by Joshi and coworkers in a different context [148]. There the model was shown to be equivalent to a transverse field anisotropic XY model, whose NESS was calculated using a matrix product operator representation. Quantum entanglement and correlations near the critical point (of the equilibrium system) were found to reveal interesting similarities and differences to the equilibrium physics. For instance, singular behaviour was found in the isotropic limit where the range of entanglement diverges (with vanishing magnitude), which is very distinct from the ground state behaviour. We refer the interested reader to the original article.

The Majorana modes are exotic 'half-a-fermion' zero energy modes that can form, among others, in a 1D topological p-wave superconductor described by Kitaev's toy model [149]. The latter involves the so-called p-wave paring term, ${{c}_{i}}{{c}_{i+1}}+\text{H}.\text{c}.$, where c_i is the fermion annihilation operator at site i, and is crucial for the formation of Majorana modes. To realise Kitaev's model with photons, the authors employ the hard-core limit ($U\gg J$) of the BH model and a parametric driving of the inter-cavity field as shown in figure 20. In the hard-core limit ($U\gg | \Delta |$), the effect of the parametric driving (with an amplitude $\Delta =| \Delta |{{\text{e}}^{\text{i}\phi}}$) is to introduce the 'squeezing' term

$$\begin{eqnarray}&&{{H}_{\text{drive}}}=\,| \Delta |\underset{i}{\overset{N-1}{\sum}}\,\left({{\text{e}}^{\text{i}\left(2{{\omega}_{p}}t+\phi \right)}}{{\tilde{{a}}\,}_{i}}{{\tilde{{a}}\,}_{i+1}}+\text{H}.\text{c}.\right),\end{eqnarray} \tag{ 24 }$$

where ${{\tilde{a}}_{i}}$s are the hard-core photon operators. Using the Jordan–Wigner mapping ${{c}_{i}}={\prod}_{j=1}^{i-1}\left(-\sigma _{j}^{z}\right)\sigma _{i}^{-}$ with $\sigma _{i}^{-}={{\tilde{a}}_{i}}$, the total Hamiltonian is written as

$$\begin{eqnarray}\begin{array}{*{35}{l}} H & = & -J\underset{i=1}{\overset{N-1}{\sum}}\,\left(c_{i}^{\dagger}{{c}_{i+1}}+\text{H}.\text{c}.\right)\,+\,| \Delta |\underset{i=1}{\overset{N-1}{\sum}}\,\left({{\text{e}}^{\text{i}\phi}}{{c}_{i}}{{c}_{i+1}}+\text{H}.\text{c}.\right) \\ {} & {} & -\mu \underset{i=1}{\overset{N-1}{\sum}}\,c_{i}^{\dagger}{{c}_{i}}, \end{array}\end{eqnarray} \tag{ 25 }$$

where c_i's are fermionic operators and $\mu ={{\omega}_{p}}-{{\omega}_{c}}$ is the pump-resonator detuning.

Zoom In Zoom Out Reset image size

The above Hamiltonian is identical to Kitaev's model of 1D p-wave superconductor of spinless fermions and hosts localised zero-energy Majorana modes [149]. The latter can most readily be seen in the limit $J= \Delta >0$ and $\mu =0$, in which the Hamiltonian can be written as

$$\begin{eqnarray}&&H=\text{i}J\underset{i=1}{\overset{N-1}{\sum}}\,d_{2i}^{\dagger}{{d}_{2i+1}}=-J\underset{i=1}{\overset{N-1}{\sum}}\,\sigma _{i}^{x}\sigma _{i+1}^{x},\end{eqnarray} \tag{ 26 }$$

where the Majorana operators are defined as ${{d}_{2i-1}}={{c}_{i}}+c_{i}^{\dagger}$ and ${{d}_{2i}}=-\text{i}\left({{c}_{i}}-c_{i}^{\dagger}\right)$. The Majorana operators d₁ and d_2N are absent from the Hamiltonian and therefore the energy does not depend on whether or not the non-local fermion $f=\left({{d}_{1}}+{{d}_{2N}}\right)/2$ is present. This means that the ground state is two-fold degenerate, characterised by the Majorana qubit operator

$$\begin{eqnarray}&&\sigma _{M}^{z}=\text{i}{{d}_{1}}{{d}_{2N}}=\underset{j=1}{\overset{N}{\prod}}\,\left(-\sigma _{j}^{z}\right)\sigma _{1}^{x}\sigma _{N}^{x}.\end{eqnarray} \tag{ 27 }$$

The stringlike operator $P\equiv {\prod}_{j=1}^{N}\left(-\sigma _{j}^{z}\right)$ corresponds to the parity operator of the total number of (hardcore) photons. It commutes with the Hamiltonian and is therefore conserved, meaning that the Majorana qubit is associated with the end qubits only, thus (topologically) protected from local perturbations that preserves parity.

In the proposed photonic system, local dissipation breaks the parity and therefore also the topological protection of the Majorana qubit. However, the existence of localised Majorana modes and their characteristic non-Abelian braiding do not depend on the parity protection, and can therefore be simulated. The Majorana modes can be detected by attaching two probe cavities at the two ends of the array and measuring their second order cross-correlations. Due to the so-called Majorana mediated Cooper pair splitting [150, 151], the nonlocal Majorana qubit mediates a nonlocal coherent exchange of photons between the probe cavities, and the cross-correlation function reveals strong bunching in the topologically nontrivial regime as shown in figure 21. The braiding operation can also be simulated by using the 'tunnel-braid' operation [152], where the Majorana edge mode in one end of the array is tunnel-coupled, via an intermediate cavity, to the Majorana edge mode in another array. A sequence of operations such as turning on and off the tunneling amplitude to the intermediate cavity can be performed to simulate the braiding operation.

Zoom In Zoom Out Reset image size

There are many more works that address interesting out-of-equilibrium physics arising in CRAs, which we did not have the space to cover in detail. Here, we provide a brief survey of these works for the interested reader.

Photon statistics in driven-dissipative BH or JCH dimers have been investigated by many authors. One of the main themes is how the photon statistics, as measured by the local intensity–intensity correlation function g⁽²⁾, change with the nonlinearity-hopping ratio. Gerace and coworkers showed that photon correlations can be used to reveal strongly correlated physics in a photonic Josephson junction—a three-site array with two end cavities coherently driven and the middle site containing single-photon nonlinearity [153]. Ferretti et al have studied a homogeneously driven case for a BH dimer and found a crossover between coherent and anti-bunched regimes [154]. An independent study on the JCH dimer was carried out by Leib and Hartmann, by mapping the system to an attractive BH model [155]. Driving one of the cavities and measuring g⁽²⁾ on the other, they observed the same crossover. A more detailed study on photon 'transport' was carried out by Biella et al for up to 60 Kerr-nonlinear resonators, in which both the transmission spectrum and second-order correlation functions were evaluated [156]. Going beyond coherent driving, effects of quantised driving fields were investigated by sending two-photon input states into a BH dimer [157]. Effects of incoherent driving was also investigated by Ruiz-Rivas and coworkers, in which a spontaneous build-up of collective coherence was found [158].

The crossover from the weak to strong hopping for a dimer as well as the infinite lattice was investigated by Nissen and coworkers [159]. There, in the dimer case, the antibunching was shown to persist up to a large value of hopping ($J/g\approx 10$), for driving strengths comparable to or larger than the dissipation rate. However, a semiclassical analysis in the limit of large coordination number shows that the antibunching disappears with increasing system size. On the other hand, photon super-bunching was also shown occur in driven dissipative coupled resonator arrays, if the laser driving is tuned to the lowest-lying two-photon (polaritonic) mode in the BH (and JCH) system as shown by Grujic et al [160]. In fact, we have already observed the same behaviour while discussing fermionisation of photons. TEBD simulations of 2–7 coupled JC cavities with homogeneous driving shows that when the hopping strength J is larger than a critical value, the system exhibits photon (super) bunching, the magnitude of which decreases with increasing system size. Next, we would like to discuss an interesting effect called unconventional photon blockade. Normally, a large value of nonlinearity ($U,g\gg \gamma$) is necessary to induce photon antibunching. However, Liew and Savona discovered that antibunching can occur in a weakly nonlinear ($U<\gamma$) dimer upon making a judicious choice of parameters [161]. Soon after, the origin of this effect was identified as a destructive quantum interference between excitation pathways [162] and a similar effect was shown to arise in resonators with a ${{\chi}^{2}}$ nonlinearity [163].

Quantum entanglement in driven-dissipative CRAs has also been investigated. Angelakis et al have studied entanglement between polaritons formed at two cavities, where the latter are coupled via an empty resonator mode [164]. Furthermore, a generalisation to three sites that are interconnected via three coherently-driven waveguides modes allows the entanglement to be controlled [165, 166]. In an alternative scheme, Liew and Savona considered three-site BH system with the middle site driven and found steady-state entanglement between the spatially separated modes [167]. An interesting aspect of this work is that, just like in the authors' antibunching work, only a small amount of nonlinearity is required. Furthermore, by repeating the three site setup, a multimode entanglement can be generated [168].

Finally, we would like to mention studies in 2D systems. We have already come across the latter when we reviewed nonequilibrium signatures of the superfluid-insulator transition [133] and the photonic supersolid phase [144, 145]. The phase diagram of the driven-dissipative BH array (without the extra cross-Kerr interaction term) was studied by Le Boité and coworkers [169, 170]. While the Mott-lobe-like structures were found, the phase diagram in the nonequilibrium setting exhibits clear differences to its equilibrium counterpart. This points to rich many-body physics beyond the equilibrium regime, which arises due to the interplay between driving and dissipation. Verifications of these mean-field results are an important future direction, as predicted effects such as bistability may not survive quantum fluctuations. Issues of bistability has been studied in 1D lattices by comparing exact numerical results against mean-field results. Mendoza-Arenas et al [171] have used matrix product states to simulate up to 60 sites to observe vanishing bistability, while Wilson et al [172] showed that in a finite system, the bistability manifests itself as a switching between 'macroscopically' distinct states. Note that while 1D lattices can be numerically simulated, an efficient method for 2D arrays is still missing. Towards this direction, new numerical methods such as the 'corner-space renormalisation method' [173] and 'self-consistent projection operator theory' [132] are currently being developed.

An interesting possibility in 2D that is absent in 1D is to introduce artificial magnetic fields by inducing phase-dependent hopping terms. This allows one to implement quantum Hall-like models for photons as we have looked at earlier. Driven-dissipative signatures in such systems were studied by various authors. Umucalılar and Carusotto proposed a method to experimentally verify the fractional quantum Hall state hidden in the NESS of a weakly driven array [100]. In a similar setting, Hafezi and coworkers showed how g⁽²⁾ measurements can be used to infer the underlying fractional quantum Hall states [102]. Efforts to devise methods to experimentally measure topological invariants in integer quantum Hall systems have also been made. Ozawa and Carusotto showed that the Chern number and the Berry curvature can be measured by following the displacement of a localised excitation by a coherent driving field [174], while Berceanu et al showed that by adding a weak harmonic potential to the model, the 'momentum-space Landau levels' can be observed [175]. Alternatively, Hafezi [176] and Bardyn et al [177] showed that twisted boundary conditions can be implemented to experimentally measure the topological invariants, while Ma et al have proposed a way to measure the real space Chern number using spectroscopic signals [178]. Methods to stabilise strongly-correlated states of light using auxiliary two-level atoms have also been suggested [179, 180].

Before we finish this section, perhaps a few words on experimental realisations are in order. There are a number of promising experimental platforms to realise the JCH (or the BH) model: superconducting circuits, photonic crystal devices with quantum dots, microcavities connected by optical waveguides, polaritons in microstructures, and trapped ions. The trapped ion implementation uses phonons, the quantised vibrational motions of trapped ions, instead of photons and have been experimentally demonstrated for the dimer case [181]. 2D patterns of micropillars etched out of a semiconductor microcavity that can host exciton polaritons is a viable experimental platform especially if the polaritonic nonlinearity can be enhanced to the few-photon level. Recently, spin–orbit coupling has been demonstrated in this platform [182]. Microcavities plus optical waveguides is an innovative platform suggested by Lepert et al [48] that uses cavities formed by open microcavities closed by coated waveguide tips. An emitter placed in each resonator acts as a two-level atom, while evanescent coupling between waveguides provide the photon hopping mechanism. Photonic crystals are dielectric devices with periodically patterned holes, where missing holes called defects can act as cavities. When a semiconductor is used to make a crystal, a quantum dot (artificial atom) grown in each defect resonator can allow the resulting system to be described by the JCH model. Majumdar et al have experimentally implemented two cavities coupled to a single quantum dot in this system [183].

The platform that shows most promise is superconducting circuits. The JC dimer has been implemented by Raftery et al [135] (the BH dimer has also been implemented by Eichler et al [184]), and lattices up to 200 resonators in a Kagome geometry—without 'atoms' (superconducting qubits)—have been fabricated [46, 185]. Scanning defect microscopy has been used to experimentally map out the normal-mode structure of microwave photons in such a lattice [186]—a technique which may prove useful in studying nonequilibrium phase transitions and photonic quantum simulators. Right after the submission of this review, a report demonstrating a 1D lattice of 72 sites, each coupled to superconducting qubits, have appeared on a preprint server, demonstrating a dissipative phase transition in which the steady state changes drastically as the mean photon number is increased [187].

Quantum simulations of 1D quantum field theories with strongly interacting photons (or more precisely, polaritons) to be reviewed in this section rely on two main techniques. The first is the technique of stationary light, allowing one to trap pulses of light and describe its time evolution via an effective Hamiltonian. The second is creating strong nonlinearity, achieved by introducing specifically structured atoms to mediate interaction between photons. Both are closely tied to the phenomenon of electromagnetically induced transparency (EIT) [188–190] as we will explain shortly.

4.1.1. Stationary pulses of light.

The concept of stationary light has been proposed and demonstrated by Bajcsy et al [191]. The main idea is based on dark state polaritons [192] formed in EIT, which we will briefly explain here. Consider an electric field propagating through an atomic medium. Normally, if the electric field is resonant with an atomic transition, it will be strongly attenuated as it propagates through the medium. However, a laser field driving another transition can induce destructive quantum interference between the excitation pathways and cancel this attenuation, i.e. the extra laser field may induce transparency. In EIT, this is accompanied with a steep dispersion, resulting in a dramatic reduction of the group velocity and enhanced nonlinearity.

Dark state polaritons: Dark state polaritons are collective excitations of light and matter formed when a quantum field propagates under specific conditions that allow EIT. A typical atomic level scheme employed in realising EIT, called a $\Lambda$-scheme, is depicted in figure 22. It comprises of a dissipative excited state $|b\rangle$ and two steady states $|a\rangle$ and $|c\rangle$. One transition to the excited state is coupled to a quantum field $\hat{E}$ while the other is coupled to a classical control (laser) field of Rabi frequency $\Omega$. The Hamiltonian describing the evolution of the atom-field system is

$$\begin{eqnarray}\begin{array}{*{35}{l}} H & = & -n{\int}^{}\text{d}z\left[\Delta {{{\hat{\sigma}}}_{bb}}+\delta {{{\hat{\sigma}}}_{cc}}+g\sqrt{2\pi}\left({{{\hat{\sigma}}}_{ba}}\hat{E}+\text{H}.\text{c}.\right) \right. \\ {} & {} & \left. +\,\left(\Omega {{{\hat{\sigma}}}_{bc}}+\text{H}.\text{c}.\right)\right], \end{array}\end{eqnarray} \tag{ 28 }$$

where n is the linear density of atoms, g is the single photon-atom coupling strength, and ${{\hat{\sigma}}_{ij}}\equiv {{\hat{\sigma}}_{ij}}(z,t)$ are continuous average spin-flip operators in a small region around the position z. The propagation of the electric field is governed by the equation

$$\begin{eqnarray}&&\left({{\partial}_{t}}\,+\,v{{\partial}_{z}}\right)\hat{E}(z,t)=\text{i}\sqrt{2\pi}gn{{\hat{\sigma}}_{ab}}(z,t),\end{eqnarray} \tag{ 29 }$$

where v is the speed of light in the absence of atoms.

Zoom In Zoom Out Reset image size

Assuming that atoms are mostly in the ground state $|a\rangle$ and adiabatically eliminating the fast-varying atomic coherence ${{\hat{\sigma}}_{ab}}$ in the Heisenberg–Langevin equation, the propagation equation reduces to (for a derivation, see for example [193])

$$\begin{eqnarray}&&\left({{\partial}_{t}}\,\,+\,\,{{v}_{g}}{{\partial}_{z}}\right)\rm{\hat \Psi }(z,t)=\frac{\text{i}2\pi {{g}^{2}}n{{v}_{g}}\delta}{\left({{ \Omega }^{2}}-\text{i}\delta \Gamma \right)v}\rm{\hat \Psi }(z,t).\end{eqnarray} \tag{ 30 }$$

The field operator $\rm{\hat \Psi }$ describes a collective excitation of the quantum field $\hat{E}$ and the atomic excitation ${{\hat{\sigma}}_{ab}}$ which is called the dark state polariton. $\Gamma =\gamma /2-\text{i} \Delta$ where γ is the decay rate from the excited state $|b\rangle$ and ${{v}_{g}}\approx v{{ \Omega }^{2}}/\left({{ \Omega }^{2}}+2\pi {{g}^{2}}n\right)$ in the limit ${{ \Omega }^{2}}\gg |\delta \Gamma |$, with δ the two-photon detuning as depicted in figure 22. The dark state polariton operator reduces to $\rm{\hat \Psi }\approx g\sqrt{2\pi n}\hat{E}(z,t)/ \Omega$ if the energy is mostly stored in the form of atomic excitations, which happens when $2\pi {{g}^{2}}n\gg {{ \Omega }^{2}}$. In such a regime the group velocity v_g is much smaller than the speed of light in the empty waveguide, v.

Stationary light: Stationary light is formed when there are two counter-propagating control fields that form a Bragg grating for the quantum fields propagating in the opposite direction. The dark state polaritons are then trapped in the system while maintaining a finite photonic component, following the intensity profile of the control fields. Being stationary, these polaritons have longer time to interact with each other, which effectively enhances the nonlinearity of the system. To derive the equation of motion for these stationary polaritons, we need to first introduce some notations. Counter propagating fields will be denoted as ${{\hat{E}}_{\pm }}$ and ${{ \Omega }_{\pm }}$, where $\hat{E}={{\hat{E}}_{+}}{{\text{e}}^{\text{i}{{k}_{\text{c}}}z}}+{{\hat{E}}_{-}}{{\text{e}}^{-\text{i}{{k}_{\text{c}}}z}}$ and $\Omega ={{ \Omega }_{+}}{{\text{e}}^{\text{i}{{k}_{\text{c}}}z}}+{{ \Omega }_{-}}{{\text{e}}^{-\text{i}{{k}_{\text{c}}}z}}$. k_c is the wavenumber of the control field, which we assume to be similar to the wavenumber of the quantum field. In turn, this induces a slowly varying optical coherence ${{\hat{\sigma}}_{ab,\pm }}$ defined through ${{\hat{\sigma}}_{ab}}={{\hat{\sigma}}_{ab,+}}{{\text{e}}^{\text{i}{{k}_{0}}z}}+{{\hat{\sigma}}_{ab,-}}{{\text{e}}^{-\text{i}{{k}_{0}}z}}$. Because of the wavenumber mismatch, there is an additional term $\text{i} \Delta \omega$ in the rhs of the propagation equation (29), where $\Delta \omega =v\left({{k}_{0}}-{{k}_{\text{c}}}\right)$ (assuming the same speed of light v for both the quantum and classical light) [194].

If the atoms are hot enough, we can use the secular approximation, i.e. neglect the fast oscillating contributions proportional to ${{\text{e}}^{\pm \text{i}2{{k}_{0}}z}}$ to obtain propagation equations similar to (30) for counter propagating dark state polaritons ${{\rm{\hat \Psi }}_{\pm }}=g\sqrt{2\pi n}{{\hat{E}}_{\pm }}/{{ \Omega }_{\pm }}$, but with extra terms that couple the counter propagating fields [191, 193]. When ${{ \Omega }_{+}}={{ \Omega }_{-}}\equiv {{ \Omega }_{0}}$, these equations can be grouped into the equations for the sum, ${{\rm{\hat \Psi }}_{\text{S}}}=\left({{\rm{\hat \Psi }}_{+}}+{{\rm{\hat \Psi }}_{-}}\right)/\sqrt{2}$, and the difference, ${{\rm{\hat \Psi }}_{\text{D}}}=\left({{\rm{\hat \Psi }}_{+}}-{{\rm{\hat \Psi }}_{-}}\right)/\sqrt{2}$ modes:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\partial}_{t}}{{{\rm{\hat \Psi }}}_{\text{S}}}+2{{v}_{g}}{{\partial}_{z}}{{{\rm{\hat \Psi }}}_{\text{D}}} & = & \text{i}\left(\delta -2 \Delta \omega \frac{{{v}_{g}}}{v}\right){{{\rm{\hat \Psi }}}_{\text{S}}} \end{array}\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\partial}_{t}}{{{\rm{\hat \Psi }}}_{\text{D}}}+v{{\partial}_{z}}{{{\rm{\hat \Psi }}}_{\text{S}}} & = & -\frac{2\pi {{g}^{2}}n}{ \Gamma }{{{\rm{\hat \Psi }}}_{\text{D}}}. \end{array}\end{eqnarray} \tag{ 32 }$$

In the limit of interest, $2\pi {{g}^{2}}n\gg \Gamma$, the difference mode ${{\rm{\hat \Psi }}_{\text{D}}}$ can be adiabatically eliminated to give

$$\begin{eqnarray}&&\text{i}{{\partial}_{t}}{{\rm{\hat \Psi }}_{\text{S}}}=-\frac{1}{2{{m}_{\text{eff}}}}\partial _{z}^{2}{{\rm{\hat \Psi }}_{\text{S}}}+{{V}_{\text{eff}}}{{\rm{\hat \Psi }}_{\text{S}}},\end{eqnarray} \tag{ 33 }$$

where ${{m}_{\text{eff}}}=-\pi {{g}^{2}}n/\left[v{{v}_{g}}\left(\Delta +\text{i}\gamma /2\right)\right]$ and ${{V}_{\text{eff}}}=2 \Delta \omega {{v}_{g}}/v-\delta$ can be thought of as a (complex) effective mass and an effective potential in analogy to the Schrödinger equation. Note that the group velocity of the sum mode is zero (it becomes non-zero for non-identical control fields ${{ \Omega }_{\pm }}$), indicating that we have a stationary light.

Assuming $\Delta =-| \Delta |\gg \gamma$, the effective mass becomes real and can be written as ${{\gamma}_{1\text{D}}}n/4{{v}_{g}}| \Delta |$, where ${{\gamma}_{1\text{D}}}=4\pi {{g}^{2}}/v$ can be identified as the spontaneous emission rate of a single atom into the waveguide. $\eta ={{\gamma}_{1\text{D}}}/\gamma$ then defines a 'single-atom cooperativity'—the efficiency of coupling into the waveguide. With these definitions the optical depth of the waveguide is defined as $OD=\eta nl$, where l is the length of the waveguide. OD quantifies the strength of the interaction between photons and the medium and a large value of it ($\approx 1000$) is typically required to get into the strongly interacting regime considered in this section.

Preparation and detection of stationary pulses of light: Quantum simulations using stationary light largely consists of three steps: (1) loading a traveling pulse, (2) converting it to a stationary pulse that undergoes a prescribed evolution, and (3) detecting the resulting stationary light. Here let us illustrate this with a typical example. A quantum pulse of light ${{\hat{E}}_{+}}$ is sent in from the left with only the control field ${{ \Omega }_{+}}$ initially turned on. In this typical EIT scheme, the quantum pulse travels slowly with the group velocity ${{v}_{g}}\propto {{ \Omega }^{2}}$. When the pulse has completely entered the medium $\Omega (t)$ is adiabatically turned to zero, trapping the entire pulse in the system [192]. Following this loading process, stationary light is created by adiabatically increasing both ${{ \Omega }_{+}}(t)$ and ${{ \Omega }_{-}}(t)$ simultaneously, creating the long-lived sum mode ${{\rm{\hat \Psi }}_{\text{S}}}$ which evolves under an engineered Hamiltonian of choice. After a dynamical evolution or an adiabatic change of the initial state, the pulse is released by turning off ${{ \Omega }_{-}}$ which converts the stationary pulse to a propagating one. The out-going pulse can then be detected using standard quantum optical measurements; for example, a photon correlation measurement can detect spatial correlations of the stationary pulse that have been converted into temporal correlations.

4.1.2. EIT nonlinearities with four-level atoms.

The last essential ingredient in simulating strongly correlated effects with photons is a strong nonlinearity (between photons or more accurately polaritons). Normally, to get photons to interact with each other strongly, one requires them to be resonant with the atoms, which unfortunately entails strong dissipation. As we have already hinted, an effective way to overcome this problem is to use EIT to suppress the linear absorption while enhancing the Kerr-nonlinearity through a quantum interference effect [190]. However, the problem is that the magnitude of the nonlinearity is limited by dissipation in the usual three-lv scheme.

A scheme that uses four level atoms to generate a giant resonantly enhanced nonlinearities was first introduced in [195], where the atomic structure depicted in figure 23 was shown to yield a cross-Kerr nonlinearity between the fields ${{\hat{E}}_{1}}$ and $\widehat{{{E}_{2}}}$. The imaginary part of the third-order susceptibility ${{\chi}^{(3)}}$ can be made small by increasing the detuning ${{ \Delta }_{2}}$, while leaving the real part finite and, in fact, orders of magnitude larger than those obtained from the usual 3-lv scheme. Self-Kerr nonlinearity can also be obtained by replacing ${{\hat{E}}_{2}}$ with ${{\hat{E}}_{1}}$ [70]. These types of nonlinearity are exactly what one needs to simulate strongly interacting fields with photons as we will explain in the remainder of this section. We will continue to adopt the Heisenberg–Langevin equation approach to study the dynamics in the low-loss regime in this review, but a master equation description can also be employed as shown by Kiffner and Hartmann [196].

Zoom In Zoom Out Reset image size

Combining the techniques of stationary light and giant Kerr nonlinearity, Chang et al [52] showed that photons can crystallise and mimic the behaviour of strongly interacting 1D bosons called the Tonks–Girardeau (TG) gas [127, 197]. The system consists of an ensemble of four-level atoms interfaced with a tight-confining waveguide such as a tapered optical fibre or a hollow-core optical fibre. There has also been a proposal to observe a TG gas induced by a purely-dissipative interaction [198] in the same setting [199], but we will not discuss it here. Let us start with a gentle introduction to the physics of 1D Bose gas before delving into the proposal of Chang and coworkers.

Studies have shown that In 1 (spatial) dimension, effects of quantum fluctuations are enhanced and interacting particles show behaviour dramatically different from the higher dimensional analogues [200]. A famous example is the TG gas, where strongly interacting bosons are effectively fermionised. Dynamics of the simplest interacting bosons in 1D are governed by the Lieb–Liniger model [201], described by the Hamiltonian

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} -{\int}^{}\text{d}z{{ \Psi }^{\dagger}}(z)\frac{1}{2m}\partial _{z}^{2} \Psi (z)+U{{ \Psi }^{\dagger}}(z){{ \Psi }^{\dagger}}(z) \Psi (z) \Psi (z), \end{array}\end{eqnarray} \tag{ 34 }$$

which for a fixed linear density n is characterised by a single dimensionless parameter ${{\gamma}_{\text{LL}}}=nU/\left({{n}^{2}}/2m\right)$. When ${{\gamma}_{\text{LL}}}\gg 1$, the bosons are impenetrable, which means that there is an effective exclusion principle and the system behaves in many ways like a non-interacting Fermi gas. This fermionisation can be observed in the second-order density–density correlation function ${{g}^{(2)}}\left(z,{{z}^{\prime}}\right)=\langle I(z)I\left({{z}^{\prime}}\right)\rangle /\left(\langle I(z)\rangle \langle I\left({{z}^{\prime}}\right)\rangle \right)$, for example, as pointed out by Chang and coworkers. The TG gas has been realised with cold atoms in optical lattices [202, 203], but a realisation with photons brings in a new set of techniques such as the standard photon-correlation measurements.

A schematic diagram of the proposed setup to realise a TG gas of photons is depicted in figure 24. The two-photon detuning δ is assumed to be zero and each atomic coherence is coupled to two counter-propagating fields as in the stationary light setting. Furthermore, the atomic transition $c\leftrightarrow d$ is coupled to the same quantum fields as the transition $a\leftrightarrow b$ so that the self-Kerr nonlinearity is induced. In the absence of the fourth lv $|d\rangle$, we have already seen that the sum mode obeys the Schrödinger equation (33). With the fourth level, this equation becomes a nonlinear Schrödinger equation (assuming ${{V}_{\text{eff}}}=0$) [52]

$$\begin{eqnarray}&&\text{i}{{\partial}_{t}}\rm{\hat \Psi }(z,t)=-\frac{1}{2{{m}_{\text{eff}}}}\partial _{z}^{2}\rm{\hat \Psi }(z,t)+2\tilde{g}{{\rm{\hat \Psi }}^{\dagger}}(z,t){{\rm{\hat \Psi }}^{2}}(z,t),\end{eqnarray} \tag{ 35 }$$

where ${{m}_{\text{eff}}}=-\pi {{g}^{2}}n/\left[v{{v}_{g}}\left({{ \Delta }_{1}}+\text{i}\gamma /2\right)\right]$ as before and $2\tilde{g}=2\pi {{g}^{2}}{{v}_{g}}/\left[\left({{ \Delta }_{2}}+\text{i}\gamma /2\right)v\right]$. We have assumed for convenience that the two excited states have the same decay rate γ. Because the dark state polariton operators obey the bosonic commutation relations, the above equation describes the evolution under the Lieb–Liniger Hamiltonian (34). In particular, the dimensionless parameter ${{\gamma}_{\text{LL}}}$ is written in terms of bare optical parameters as

$$\begin{eqnarray}&&{{\gamma}_{\text{LL}}}=\frac{{{m}_{\text{eff}}}\tilde{g}}{{{n}_{\text{ph}}}}=-\frac{{{\pi}^{2}}{{g}^{4}}}{\left({{ \Delta }_{1}}+\text{i}\gamma /2\right)\left({{ \Delta }_{2}}+\text{i}\gamma /2\right)}\frac{n}{{{n}_{\text{ph}}}},\end{eqnarray} \tag{ 36 }$$

where n_ph is the linear density of photons at the centre of the pulse. Notice that there is a non-zero imaginary part, which will be neglected on the basis that $|{{ \Delta }_{1,2}}|\gg \gamma$.

Zoom In Zoom Out Reset image size

Among the characteristic features of a TG gas is an oscillating density–density correlation function, similar to the Friedel oscillations in electronic systems with an impurity [204]. In the TG limit, the second-order intensity correlation function ${{g}^{(2)}}\left(z,{{z}^{\prime}}\right)=\langle I(z)I\left({{z}^{\prime}}\right)\rangle /\left(\langle I(z)\rangle \langle I\left({{z}^{\prime}}\right)\rangle \right)$ of the ground state takes the form [205]

$$\begin{eqnarray}&&g_{\text{TG}}^{(2)}\left(z,{{z}^{\prime}}\right)=1-{{\left(\frac{\sin {{k}_{\text{F}}}\left(z-{{z}^{\prime}}\right)}{{{k}_{\text{F}}}\left(z-{{z}^{\prime}}\right)}\right)}^{2}},\end{eqnarray} \tag{ 37 }$$

where $I(z)={{\rm{\hat \Psi }}^{\dagger}}(z)\rm{\hat \Psi }(z)$ and ${{k}_{\text{F}}}=\pi {{n}_{\text{ph}}}$ is the 'Fermi momentum' of the photonic TG gas.

The 2k_F oscillations can be observed experimentally by following the preparation and detection procedure outlined in section 4.1.1: an initial pulse of coherent state is loaded into the waveguide and turned into a stationary pulse. Then the dimensionless interaction parameter ${{\gamma}_{\text{LL}}}$ is increased from zero to a large value by changing, for example, the detuning. The initial coherent state, which is an eigenstate of the free Hamiltonian, remains in an eigenstate of the Hamiltonian (with non-zero interaction strength) and therefore exhibits the Friedel oscillations. The stationary pulse can then be converted to a traveling pulse, whose temporal correlation function reveals the density–density correlation function of the simulated quantum state.

An example is shown in figure 25, taken from [52], for an initial density profile ${{n}_{\text{ph}}}(z)={{n}_{0}}{{\left(1-{{z}^{2}}/z_{0}^{2}\right)}^{1/2}}$ under an exponentially increasing interaction parameter ${{\gamma}_{\text{LL}}}$. The distance is in units of the inverse effective Fermi wave vector k_F and the final time is $t=10\omega _{f}^{-1}$, where ${{\omega}_{f}}$ is the Fermi energy on the order of $n_{\text{ph}}^{2}/{{m}_{\text{eff}}}$. For an optical depth of $\approx 2000$ and the cooperativity $\eta \approx 0.2$, the achievable maximum ${{\gamma}_{\text{LL}}}$ is around 10.

Zoom In Zoom Out Reset image size

Looking back to equation (33), it is not too difficult to see that a spatially dependent two-photon detuning would induce an effective potential for the stationary light. Indeed, this fact has been noticed in a proposal to simulate the Klein tunnelling with stationary light [22]. Upon introducing another transition, one can induce a small spatially-dependent two-photon detuning and as a consequence the resulting nonlinear Schrödinger equation would have an effective potential term. Another method is to induce a small change in the spatial profile of the atomic density [206]. The latter requires cold atomic gas so that the density profile does not change appreciably during the experiment. However, evolution of a stationary light in a cold atomic ensemble is different to that in a hot vapour because of the secular approximation discussed earlier in section 4.1.1 [207, 208] and one has to resort to a different atomic level scheme of a double-lambda type. Here, we discuss the effects of introducing a lattice potential to strongly interacting photons assuming that the lattice has been created via one of the two methods and defer a discussion on stationary light in cold atomic gases to the next section.

We follow Huo and Angelakis [206] and consider the Lieb–Liniger model with an additional potential term of the form $V(z)={{V}_{0}}{{\cos}^{2}}\left(\pi {{n}_{\text{ph}}}z\right)$. The periodicity of the lattice potential is set such that the number of sites equals the number of photons in the waveguide. Depending on the lattice depth and the interaction strength, the system shows either the Mott insulating phase or the superfluid phase [209] as depicted schematically in figure 26.

Zoom In Zoom Out Reset image size

There are two interesting limits in the phase diagram: (i) weakly interacting particles (${{\gamma}_{\text{LL}}}<1$) in a deep potential (${{V}_{0}}\gg {{E}_{r}}$) and (ii) strongly interacting particles (${{\gamma}_{\text{LL}}}>1$) in a shallow potential (${{V}_{0}}\sim {{E}_{r}}$), where ${{E}_{r}}={{\pi}^{2}}n_{\text{ph}}^{2}/2{{m}_{\text{eff}}}$ is the recoil energy. For case (i), the physics is well described by the well-known Bose–Hubbard model [75, 210], where the lowest localised Wannier mode in each lattice site is coupled to the modes in its nearest neighbour sites with depth-dependent on-site interaction strength. As we have seen earlier, the Bose–Hubbard hamiltonian reads

$$\begin{eqnarray}&&{{H}_{\text{BH}}}=-J\underset{\langle i,j\rangle}{\sum}\,a_{i}^{\dagger}{{a}_{j}}+\frac{U}{2}\underset{i}{\sum}\,a_{i}^{\dagger}a_{i}^{\dagger}{{a}_{i}}{{a}_{i}},\end{eqnarray} \tag{ 38 }$$

where a_i is now interpreted as the annihilation operator for the Wannier mode at site i, $\sqrt{\pi}J/{{E}_{r}}=4{{\left({{V}_{0}}/{{E}_{r}}\right)}^{3/4}}{{\text{e}}^{-2\sqrt{{{V}_{0}}/{{E}_{r}}}}}$, and $U/{{E}_{r}}=\sqrt{2/{{\pi}^{3}}}{{\left({{V}_{0}}/{{E}_{r}}\right)}^{1/4}}{{\gamma}_{\text{LL}}}$. Note that the effective parameters J and U only depend on V₀, E_r, and ${{\gamma}_{\text{LL}}}$, all of which can be tuned by changing quantum optical parameters. By tuning the value of U/E_r, through changing the value of the detuning ${{ \Delta }_{2}}$ for example, one can observe the Mott-to-superfluid transition.

There is also an interesting quantum phase transition called a 'pinning' transition in case (ii), when the interaction is strong. In the weak lattice limit, the particle density is weakly modulated and one is able to separate out the density perturbation and phase as canonically conjugate observables. Haldane was the first to consider this case and show that the quantum sine-Gordon model correctly describes the physics in this limit [211]. Following the work of Büchler, Blatter, and Zwerger [209], we briefly reproduce the derivation of the sine-Gordon Hamiltonian here.

First consider a hydrodynamic approximation where one describes the bosonic field operator in the low energy limit as $\Psi \approx \sqrt{{{\rho}_{0}}+{{\partial}_{z}}\theta /\pi}\exp \left(\text{i}\phi \right)$, where ${{\rho}_{0}}$ is the mean density of the atoms, ${{\partial}_{z}}\theta /\pi$ is the long-wavelength density fluctuation, and ϕ is the phase field. The latter two obeys the boson commutation relation $\left[{{\partial}_{z}}\theta (z),\phi \left({{z}^{\prime}}\right)\right]=\text{i}\pi \delta \left(z-{{z}^{\prime}}\right)$. Neglecting the lattice potential for now, the low-energy physics is described by the Hamiltonian

$$\begin{eqnarray}&&{{H}_{0}}=\frac{1}{2\pi}{\int}^{}\text{d}z\left[{{v}_{J}}{{\left({{\partial}_{z}}\phi \right)}^{2}}+{{v}_{N}}{{\left({{\partial}_{z}}\theta \right)}^{2}}\right].\end{eqnarray} \tag{ 39 }$$

The first part comes from the kinetic energy where ${{v}_{J}}=\pi {{\rho}_{0}}/m$ and the second term derives from the interaction energy with the coefficient ${{v}_{N}}={{\partial}_{n}}\mu /\pi$ determined by the inverse compressibility. To obtain the Hamiltonian due to the lattice potential, one needs to modify the density operator to account for its discrete character:

$$\begin{eqnarray}&&\rho (z)=\left[{{\rho}_{0}}+\frac{1}{\pi}{{\partial}_{z}}\theta \right]\left[1+2\underset{l=1}{\overset{\infty}{\sum}}\,\cos \left(l\theta +l\pi {{\rho}_{0}}z\right)\right].\end{eqnarray} \tag{ 40 }$$

In the commensurate limit (one particle per lattice cell), the dominant part of the potential arising from the lowest harmonic in the above equation takes the sine-Gordon form

$$\begin{eqnarray}&&{{H}_{V}}=\frac{{{n}_{\text{ph}}}{{V}_{0}}}{2}{\int}^{}\text{d}z\cos (2\theta ).\end{eqnarray} \tag{ 41 }$$

The Hamiltonian ${{H}_{0}}+{{H}_{V}}$ is the 1D quantum sine-Gordon equation which is known to exhibit the so-called 'pinning' transition [212]. This phase transition can best be described by the dimensionless parameter $K=\sqrt{{{v}_{J}}/{{v}_{N}}}$ that is related, in the absence of the periodic potential, to the Lieb–Liniger parameter ${{\gamma}_{\text{LL}}}$; $K(\gamma \leqslant 10)\approx \pi /\sqrt{{{\gamma}_{\text{LL}}}-(1/2\pi )\gamma _{\text{LL}}^{3/2}}$ and $K(\gamma \gg 10)\approx {{\left(1+2/{{\gamma}_{\text{LL}}}\right)}^{2}}$. When K  >  K_c  =  2, the weak potential is irrelevant and the bosons are in the superfluid phase whereas for K  <  K_c an arbitrarily weak potential can pin the bosons to lattice sites, achieving the Mott insulating state.

The procedure to observe these phases in the nonlinear waveguide system with a lattice potential is as follows. First a pulse containing a certain number of photons is loaded into the waveguide. Then one turns up the interaction strength and then the lattice potential adiabatically. This step can be realised experimentally by changing the detunings ${{ \Delta }_{2}}$ and δ. By the end of this step, the photonic state would have adiabatically followed the Hamiltonian and consequently would exhibit either the superfluid or the Mott-insulating behaviour which can be observed by releasing the stationary pulse and performing correlation measurement on them as in the observation of a TG gas of photons.

The realisation of strongly interacting photons in the previous section relied on having a hot vapour of atoms so that the secular approximation holds. However, it is also possible to create stationary light with cold atoms [213, 214] as demonstrated in [215]. In particular, Zimmer et al [214] showed how the linear Schrödinger equation can be obtained with cold atoms using a slightly different atomic level scheme shown in figure 27(a). In this scheme, ${{\hat{E}}_{+}}$ and ${{\hat{E}}_{-}}$ couple to different atomic transitions unlike in the usual $\Lambda$ scheme where the two couple to the same transition, so the secular approximation no longer needs to be invoked. For equal control field strengths, ${{ \Omega }_{+}}={{ \Omega }_{-}}$, the stationary dark state polariton is defined as

$$\begin{eqnarray}&&\rm{\hat \Psi }\left(\mathbf{r},t\right)=\frac{\cos \theta}{\sqrt{2}}\left({{\hat{E}}_{+}}\left(\mathbf{r},t\right)+{{\hat{E}}_{-}}\left(\mathbf{r},t\right)\right)-\sin \theta \hat{S}\left(\mathbf{r},t\right),\end{eqnarray} \tag{ 42 }$$

where ${{\tan}^{2}}\theta =2\pi {{g}^{2}}n/{{ \Omega }^{2}}$, ${{ \Omega }^{2}}= \Omega _{+}^{2}+ \Omega _{-}^{2}$, and $\hat{S}=\sqrt{n}{{\hat{\sigma}}_{ac}}\left(\mathbf{r},t\right)$. The Schrödinger equation obeyed by this field operator reads

$$\begin{eqnarray}&&\text{i}{{\partial}_{t}}\rm{\hat \Psi }\left(\mathbf{r},t\right)=-\left[\frac{1}{2{{m}_{\parallel}}}\partial _{z}^{2}\,+\,\frac{1}{2{{m}_{\bot}}}{{ \Delta }_{\bot}}\right]\rm{\hat \Psi }\left(\mathbf{r},t\right),\end{eqnarray} \tag{ 43 }$$

where ${{m}_{\bot}}=m{{v}_{\text{rec}}}/{{v}_{g}}$ and $m_{\parallel}^{-1}=m_{\bot}^{-1}{{k}_{p}}\left(\gamma v/2\pi {{g}^{2}}n\right)\left(2 \Delta /\gamma +\text{i}\right)$. Here k_p is the wavenumber of the probe field and ${{v}_{\text{rec}}}={{k}_{p}}/m$ is the recoil velocity of an atom of mass m.

Zoom In Zoom Out Reset image size

Introduction of additional levels as in figure 27(b) adds the Kerr-type interaction between the dark state polaritons whose elastic collision rate can be much greater than the loss rates and thereby allow thermalisation of polaritons without significant losses [216]. The optical depth required to reach this regime is only modestly high at about OD  >  30. For realistic parameters, the critical temperature to create a BEC out of dark state polaritons (as compared to the critical temperature for ideal bosonic atoms) is ${{T}_{\text{c}}}/T_{\text{c}}^{\text{atoms}}\approx 6\times {{10}^{5}}$ which is in the mK regime. The latter is in the same range as the temperature limit of about 1 mK imposed by the finite frequency window of EIT. Therefore, it is possible to observe the BEC of stationary light polaritons at temperatures several orders of magnitude higher than the critical temperature of atomic BECs.

In the stationary light setup, a BEC will be created not by reducing the temperature of the polariton gas, but by increasing the critical temperature of condensation through changing the effective mass. The steps to observe Bose–Einstein condensation of polaritons is then as follows. (i) A pulse of coherent light is trapped as stationary polaritons using the same method described in the previous section. Initially, the critical temperature is below the ambient temperature so there is no condensation. (ii) As the critical temperature is increased, by tuning the strength of the control laser, the polaritons thermalise through elastic collisions and reach BEC before significant losses occur. This causes the macroscopic occupation of $\mathbf{k}=0$ momentum state, which can be (iii) observed by converting the stationary polariton into a traveling one, whose transverse profile shows the onset of BEC. The latter is illustrated in figure 28, where the transverse emission profiles are shown at temperatures above (left) and below (right) the critical temperature. A similar proposal with Rydberg atoms—to induce dipolar interaction between dark state polaritons—have been made by Nikoghosyan et al [217].

Zoom In Zoom Out Reset image size

So far in this section, we have reviewed works on quantum simulations of strongly interacting bosons when there is only a single species of bosons. Soon after the pioneering work on a TG gas of stationary light, Angelakis et al realised that the physics of two interacting bosons can also be simulated in a similar waveguide setup [53, 218]. Here, we describe how this leads to a realisation of the two-component Lieb–Liniger model, which could allow a clear experimental observation of spin-charge separation. Furthermore, by introducing a lattice potential, one can also simulate a two-component Bose–Hubbard model, which exhibits interesting behaviours such as (a 1D analogue of) the famous BEC-BCS crossover [219, 220].

A schematic diagram of the atomic level scheme that yields the two-component Lieb–Liniger model is shown in figure 29. Two different quantum fields ${{\hat{E}}_{1}}$ and ${{\hat{E}}_{2}}$ with frequencies $\omega _{qu}^{(1)}$ and $\omega _{qu}^{(2)}$ interact with two different species of 4-lv atoms labeled a and b and forms two polariton operators via the usual relation ${{\rm{\hat \Psi }}_{j,\pm }}=g\sqrt{2\pi {{n}^{{{x}_{j}}}}}{{\hat{E}}_{j,\pm }}/{{ \Omega }_{{{x}_{j}}}}$. The figure depicts the level scheme for type-a atoms; atoms of type b have a similar coupling scheme except that the roles of ${{\hat{E}}_{1}}$ and ${{\hat{E}}_{2}}$ are reversed and the detunings are accordingly labeled by replacing a with b. We use $j=1,2$ to denote the polariton species, ${{x}_{1}}=a,{{x}_{2}}=b$ to label the atomic species, and we have assumed that the atom-photon coupling strength is independent of the atomic species and transition. n^a (n^b) denotes the density of atomic species a (b) and ${{\delta}_{ij}}=\omega _{qu}^{(i)}-\omega _{qu}^{(\,j)}$. By setting the detunings such that ${{\hat{E}}_{1}}$ (${{\hat{E}}_{2}}$) is only coupled to atoms of type-b (a) through the transition $|3\rangle \leftrightarrow |4\rangle$, one obtains a cross-species interaction term on top of the self-interaction terms. The dynamics of the polaritons are thus governed by the two-component Lieb–Liniger model

$$\begin{eqnarray}&&{{H}_{\text{tcll}}}\,=\,{\int}^{}\text{d}z\underset{j}{\sum}\,\left[-\frac{1}{2{{m}_{j}}}\rm{\hat \Psi }_{j}^{\dagger}(z)\partial _{z}^{2}{{\rm{\hat \Psi }}_{j}}(z)\right. \\ &&\left. +\frac{{{U}_{j}}}{2}\rm{\hat \Psi }_{j}^{\dagger}(z)\rm{\hat \Psi }_{j}^{\dagger}(z){{\rm{\hat \Psi }}_{j}}(z){{\rm{\hat \Psi }}_{j}}(z)\right]\,\,+\,\,\frac{V}{2}{{\rho}_{1}}(z){{\rho}_{2}}(z),\end{eqnarray} \tag{ 44 }$$

where ${{\rho}_{j}}=\rm{\hat \Psi }_{j}^{\dagger}(z){{\rm{\hat \Psi }}_{j}}(z)$. m_j and U_j are defined as m_eff and $4\tilde{g}$ introduced earlier, with species-dependent optical parameters, whereas the inter-species interaction parameter V is new. For simplicity, we will assume that all species-dependence are suppressed. Then the inter-species interaction strength is compactly written as

$$\begin{eqnarray}&&V=\frac{16\pi {{g}^{2}}n{{v}_{g}}}{{{ \Delta }_{2}}}.\end{eqnarray} \tag{ 45 }$$

Zoom In Zoom Out Reset image size

An alternative scheme where the two species correspond to different polarisations has also been proposed in [218]. We will come back to this setup in section 4.6 where a quantum simulation of an interacting relativistic model will be reviewed.

4.5.1. Luttinger liquid and spin-charge separation.

We have already seen how Lieb–Liniger bosons can be described, in the low-energy sector, in terms of the long-wave density function ${{\partial}_{z}}\theta$ and the phase field ϕ. It turns out that the same decomposition applied to the two-component Lieb–Liniger model separates the Hamiltonian into two independent components, describing spin and charge components that travel in different speeds. Let us assume that the two species have identical parameters, i.e. m_j  =  m, U_j  =  U, and $\langle {{\rho}_{j}}\rangle =\rho$. Following the notations in [221] (see also [222]), the Hamiltonian (44) becomes $H={{H}_{c}}+{{H}_{s}}$, where

$$\begin{eqnarray}&&{{H}_{c}}=\frac{1}{2\pi}{\int}^{}\text{d}z{{v}_{c}}{{K}_{c}}{{\left({{\partial}_{z}}{{\theta}_{c}}\right)}^{2}}+\frac{{{v}_{c}}}{{{K}_{c}}}{{\left({{\partial}_{z}}{{\phi}_{c}}\right)}^{2}}\end{eqnarray} \tag{ 46 }$$

and

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{s}}=\frac{1}{2\pi}{\int}^{}\text{d}z{{v}_{s}}{{K}_{s}}{{\left({{\partial}_{z}}{{\theta}_{s}}\right)}^{2}}+\frac{{{v}_{s}}}{{{K}_{s}}}{{\left({{\partial}_{z}}{{\phi}_{s}}\right)}^{2}} \\ +\frac{V}{{{\left(\pi {{\rho}^{-1}}\right)}^{2}}}\cos \left(\sqrt{8{{\phi}_{s}}}\right). \end{array}\end{eqnarray} \tag{ 47 }$$

Here the subscripts c and s refer to the charge and spin components and ${{v}_{c,s}}={{v}_{0}}\sqrt{1\pm VK/\pi {{v}_{0}}}$ and ${{K}_{c,s}}=K/\sqrt{1\pm VK/\left(\pi {{v}_{0}}\right)}$, where analytic expressions ${{v}_{0}}=\sqrt{U\rho /m}$ and $K=\sqrt{{{\pi}^{2}}\rho /(mU)}$ hold in the weakly interacting regime.

The physics of this system is fully determined by the velocities v_c,s and the so-called Luttinger parameters K_c,s. Note that the speed of the spin and charge degrees of freedom can be different and could even have different signs. In terms of polaritons, the spin and charge densities correspond to the difference and the sum of the densities of two stationary polariton components. Therefore one could create an excitation in the original polariton field and observe as it splits into two components. One would go through the preparation stage discussed before, then let the two components evolve under different speeds, followed by the conversion into traveling slow-light which can be detected via standard optical measurements and post-processed to determine the velocities of the spin and charge parts, for example. In [53], it was shown that the optical depth of $\approx 3000$ is sufficient to create the difference ${{v}_{c}}/{{v}_{s}}=2$ and observe spin-charge separation by measuring the single particle spectral function, which can be determined from the first order correlation function of one of the fields, say ${{\hat{E}}_{1}}$ at a certain quasi-momentum.

4.5.2. BCS-BEC like crossover.

Introducing an effective lattice potential to the two-component setup, one can also experimentally realise a 1D two-component Bose–Hubbard model [220], which exhibits some interesting properties. For example, upon inducing strong repulsive interaction between the same species, one can go into an effective Fermi–Hubbard regime where a 1D analogue of the famous BCS-BEC crossover can be observed. This requires an attractive inter-species interaction, which is easily achieved in the stationary light setup by changing the sign of a single-photon detuning. In the strongly attractive case, the two species tend to pair up to form a molecule although no more than a pair is allowed at a single site because of the effective Pauli exclusion due to the strong intra-species repulsion. This is reminiscent of the short-range-correlated bosonic molecules forming in the BEC regime. In the weakly attractive case, two species tend to stay away from each other and form a long-range correlated pair analogous to the Cooper pair. Thus, by continuously tuning the inter-species interaction strength, it is possible to observe a 1D analogue of the BEC-BCS cross-over. A schematic illustration of these phases is provided in figure 30.

Zoom In Zoom Out Reset image size

An optical depth of approximately 1000 was shown to be sufficient for simulating the crossover, where tuning the frequency difference between the quantum fields changes the ratio between the inter- and intra-species interaction strengths. Going through the standard preparation and release processes, one could measure the cross-species intensity correlation function to clearly distinguish all three (BCS-like, BEC-like, and the big-Boson) phases. In figure 30 we plot two types of correlation functions: $g_{\uparrow \downarrow}^{2}(l)={\sum}_{i}\langle {{n}_{i\uparrow}}{{n}_{i+l\downarrow}}\rangle$ and $g_{-}^{2}(l)={\sum}_{i}\langle \left({{n}_{i\uparrow}}-{{n}_{i\downarrow}}\right)\left({{n}_{i+l\uparrow}}-{{n}_{i+l\downarrow}}\right)\rangle$. The same-site cross-species density–density correlation function $g_{\uparrow \downarrow}^{2}(0)$ clearly shows the cross-over behaviour upon changing the inter-species interaction strength, whereas a sharp transition is observed as the inter-species attraction becomes stronger than the intra-species repulsion. Correlations in the population difference, $g_{-}^{2}$, is sensitive to the BCS-BEC cross over, but not the transition to the big-Boson phase.

Simulations of quantum field theories using nonlinear waveguides can go beyond non-relativistic models. In particular, it was shown that the Thirring model (TM) [223], a relativistic quantum field theory describing the self-interaction of a Dirac field in 2 spacetime dimensions, can be implemented [54]. The atomic level scheme required is shown in figure 31(a). It consists of two sets of the standard four-level arrangements $\left\{|a\rangle,|b,s\rangle,|c,s\rangle,|d,s,{{s}^{\prime}}\rangle \right\}$, with $s,{{s}^{\prime}}=\left\{\uparrow,\downarrow \right\}$ denoting the two oppositely-polarised quantum pulses ${{\hat{E}}_{\uparrow}}$ and ${{\hat{E}}_{\downarrow}}$. Two self-interacting stationary polaritons ${{\rm{\hat \Psi }}_{s}}$ appear in exactly the same manner as discussed earlier and interact through the extra level $|d,s,\bar{s}\rangle$, $\bar{s}\ne s$, with strength dependent on ${{ \Delta }_{s,{{s}^{\prime}}}}$. Making the usual assumptions, these polaritons can be shown to obey the two-component Lieb–Liniger model [218]. To allow for the linear momentum term required in relativistic models, a small but finite difference between the forward and backward propagating control fields is necessary. Accommodating this change, the (quasi-) stationary lights are now defined as ${{\rm{\hat \Psi }}_{s}}={{\alpha}_{s,+}}{{\rm{\hat \Psi }}_{s,+}}+{{\alpha}_{s,-}}{{\rm{\hat \Psi }}_{s,-}}$ with ${{\alpha}_{s,\pm }}= \Omega _{s,\pm }^{2}/\left(\Omega _{s,+}^{2}+ \Omega _{s,-}^{2}\right)$ and ${{\rm{\hat \Psi }}_{s,\pm }}=\sqrt{2\pi n}g{{\hat{E}}_{s,\pm }}/{{ \Omega }_{s,\pm }}$. The dynamics of these fields are governed by

$$\begin{eqnarray}\begin{array}{*{20}{l}} {{\rm{i}}{\partial _t}{{{\rm{\hat \Psi }}}_s}}& = &{ - \frac{1}{{2{m_{nr,s}}}}\partial _z^2{{{\rm{\hat \Psi }}}_s} + {\rm{i}}{\eta _s}{\partial _z}{{{\rm{\hat \Psi }}}_s} + {{\rm{\Omega }}_0}{{{\rm{\hat \Psi }}}_{\bar s}}}\\ {}&{}&{ + \,{\chi _{ss}}{\rm{\hat \Psi }}_s^\dagger {\rm{\hat \Psi }}_s^2 + {\chi _{s\bar s}}{\rm{\hat \Psi }}_{\bar s}^\dagger {{{\rm{\hat \Psi }}}_{\bar s}}{{{\rm{\hat \Psi }}}_s},}\end{array}\end{eqnarray} \tag{ 48 }$$

where ${{m}_{nr,s}}=-\rm{\bar \Omega }_{s}^{2}/\left(4v_{s}^{2}{{ \Delta }_{s}}{{\sin}^{2}}2{{\varphi}_{s}}\right)$, ${{\eta}_{s}}=-2{{v}_{s}}\cos 2{{\varphi}_{s}}$, ${{\chi}_{ss}}=4\rm{\bar \Omega }_{s}^{2}/\left({{ \Delta }_{ss}}n\right)$, and ${{\chi}_{s\bar{s}}}=2\rm{\bar \Omega }_{s}^{2}\left[2+\cos \left({{\varphi}_{{\bar{s}}}}-{{\varphi}_{s}}\right)\right]/\left({{ \Delta }_{s\bar{s}}}n\right)$. v_s is the reduced group velocity for the dark state polaritons, ${{\varphi}_{s}}$ is determined from the relation ${{\tan}^{2}}{{\varphi}_{s}}= \Omega _{s,-}^{2}/ \Omega _{s,+}^{2}$, and $\rm{\bar \Omega }_{s}^{2}=\left(\Omega _{s,+}^{2}+ \Omega _{s,-}^{2}\right)/2$, and ${{ \Omega }_{0}}$ is the coupling rate between $|c,\uparrow \rangle$ and $|c,\downarrow \rangle$ as shown in figure 31(a).

Zoom In Zoom Out Reset image size

To transform the above equation into the bosonic TM, we need to throw away the quadratic dispersion term and set ${{\eta}_{\uparrow}}=-{{\eta}_{\downarrow}}=\eta$. The latter can be achieved by setting ${{v}_{\uparrow}}={{v}_{\downarrow}}$ and $\cos 2{{\varphi}_{\uparrow}}=-\cos 2{{\varphi}_{\downarrow}}$, whereas the former requires the non-relativistic kinetic energy to be much smaller than the relativistic kinetic energy, which can be achieved by an appropriate tuning of the single photon detunings ${{ \Delta }_{s}}$. Defining a spinor field $\boldsymbol{\Psi}={{\left({{\rm{\hat \Psi }}_{\uparrow}},{{\rm{\hat \Psi }}_{\downarrow}}\right)}^{T}}$, the Hamiltonian generating the equation can be written as

$$\begin{eqnarray}\begin{array}{*{35}{l}} H = & {\int}^{}\text{d}z\overline{\boldsymbol{\Psi}}\left(-\text{i}|\eta |{{\gamma}_{1}}{{\partial}_{z}}\,+\,{{m}_{0}}{{\eta}^{2}}\right)\boldsymbol{\Psi} \\ {} & +\,\underset{s}{\sum}\,\frac{{{\chi}_{ss}}}{2}\rm{\hat \Psi }_{s}^{\dagger}\rm{\hat \Psi }_{s}^{\dagger}{{{\rm{\hat \Psi }}}_{s}}{{{\rm{\hat \Psi }}}_{s}}+\frac{\chi}{2}\overline{\boldsymbol{\Psi}}{{\gamma}^{\mu}}\overline{\boldsymbol{\Psi}}{{\boldsymbol{\Psi}}}\,{{\gamma}_{\mu}}\boldsymbol{\Psi}, \end{array}\end{eqnarray} \tag{ 49 }$$

where $\overline{\boldsymbol{\Psi}}=\boldsymbol{}{{\mathbf{\Psi}}^{\dagger}}{{\gamma}_{0}}$ with the gamma matrices defined as ${{\gamma}_{0}}={{\gamma}^{0}}={{\sigma}_{x}}$, ${{\gamma}_{1}}=-{{\gamma}^{1}}=\text{i}{{\sigma}_{y}}$ and ${{m}_{0}}=-{{ \Omega }_{0}}/{{\eta}^{2}}$. Thus, neglecting the small quadratic dispersion term, the dynamics of the stationary polaritons is governed by the bosonic TM. Furthermore, by entering the hardcore regime in which the self-interaction energy dominates over all other terms, one can get into the fermionic regime where, as discussed earlier, the 1D bosons behave in many ways as fermions. Figure 32 illustrates that this regime (small kinetic energy  +  hardcore) can be reached in a possible experiment, assuming achievable values for the optical parameters: $|\cos 2{{\varphi}_{s}}|=0.004$, $n={{10}^{7}}~{{\text{m}}^{-1}}$, ${{\rm{\bar \Omega }}_{2}}\approx 1.5\gamma$, $0<\chi /|\eta |<\pi$, and the single atom cooperativity ${{\gamma}_{1\text{D}}}/\gamma =0.2$.

Zoom In Zoom Out Reset image size

One of the methods to extract interesting physical information in field theories is to calculate correlation functions. In the massless fermionic TM, the n-point correlation function defined as $\langle 0|{\prod}_{i=1}^{n}\rm{\hat \Psi }_{\uparrow}^{\dagger}\left({{z}_{i}}\right){{\rm{\hat \Psi }}_{\downarrow}}\left({{z}_{i}}\right)\rm{\hat \Psi }_{\downarrow}^{\dagger}\left(z_{i}^{\prime}\right){{\rm{\hat \Psi }}_{\uparrow}}\left(z_{1}^{\prime}\right)|0\rangle$ has been explicitly calculated [224], and serves as an excellent check for the quantum simulator. After a satisfactory agreement has been obtained in the massless case, one can increase the value of the mass and explore the correlation functions in the relatively unexplored massive regime. To see how the correlation function can be measured, it is helpful to define the fields ${{\rm{\hat \Psi }}_{x,\pm }}=\left({{\rm{\hat \Psi }}_{\uparrow}}\pm {{\rm{\hat \Psi }}_{\downarrow}}\right)/\sqrt{2}$ and ${{\rm{\hat \Psi }}_{y,\pm }}=\left({{\rm{\hat \Psi }}_{\uparrow}}\mp \text{i}{{\rm{\hat \Psi }}_{\downarrow}}\right)/\sqrt{2}$. When written in terms of these fields, the original correlation function becomes a combination of various density–density correlation functions involving the densities of x and y, which can be straightforwardly measured. The detection scheme is depicted in figure 31(b). The waveplates mix the $\uparrow$ and  ↓  components into x and y fields, one of which is then detected by using a polarisation beam splitter in each arm. In this way, one could probe the scaling behaviour of the correlation functions and moreover infer the renormalisation of mass due to the interaction.

There are largely two promising experimental techniques to produce the waveguide-atom interface considered in this section. One is loading atoms in a hollow-core photonic band gap fibres (HCPBFs) and the other is coupling atoms to a tapered-fibre or more recently a 2D photonic crystal. The purpose of this section is to discuss briefly the experimental progress in these fields.

The HCPBF is a waveguide which has a hollow-core surrounded by a photonic crystal structure [225, 226] and can therefore accommodate atoms in its core, where the guided mode resides. Strong confinement of light at the core then facilitates the sought-after strong atom-light interaction. The first experimental effort used pressure difference to pump hydrogen atoms into a HCPBF and thereby demonstrate stimulated Raman scattering [227]. It was soon followed by the generation of optical solitons in a xenon-filled fibre [228]. A step towards strong light-matter interaction has been made by demonstrating EIT using acetylene molecules [229] and rubidium molecules [230]. Optical guiding of thermal [231] and cold atoms [232] through a HCPBF and trapping of ultracold atoms in a HCPBF [233] have also been achieved. Most importantly for the subject of this review, guiding and trapping of laser cooled atoms was achieved using a magnetic guide [234]. In the same experiment, authors demonstrated all-optical switching using the four-level giant cross-Kerr effect introduced in section 4.1.1 with a moderate value of optical depth $\approx 30$. More recently, the use of hollow-optical-beam has allowed more efficient guiding of atoms, leading to an improved optical depth of 180 [235].

The second promising system consists of atoms coupled to the guided mode of a tapered fibre. When an optical fibre is drawn out such that its diameter is smaller than the wavelength of the guided light [236], light mostly resides outside the fibre in the form of an evanescent wave, which enables easy coupling between the guided light and surrounding atoms. Coupling trapped cold atoms to the guided mode of a tapered fibre has been demonstrated [237, 238] and a hot vapour of atoms has been used to demonstrate EIT with a very low pump power [239]. Trapping and interrogation of cold atoms on the surface of a tapered-fibre by creating an optical lattice above the nanofibre surface has been achieved with an optical depth as high as $\approx 30$ [240, 241]. Most recently, state-insensitive trapping of atoms has been achieved with an optical depth of 66.

More recently, people have started looking at 2D nanophotonic structures involving photonic-crystal patterns [242, 243]. These types of systems allow strong light-atom coupling and may lead to an alternative platform to realise the models reviewed in this section. Moreover, other types of models can be readily simulated in this setup, such as the conventional cavity-QED physics [244] and tunable many-body spin models [245]. As such, the platform is also a good candidate for a quantum information processor.

Before we close the section, we would like to briefly mention another interesting paradigm to realise strong photon–photon interaction that has attracted much interest recently: Rydberg EIT. In Rydberg EIT, the strong long-range Rydberg interaction between atoms are converted to an effective photon—photon interaction through the appropriate use of EIT—one chooses the third state ($|c\rangle$) of the atom to be a Rydberg state (typically requiring the level scheme to be of the ladder type). Due to an energy shift induced by the interaction between Rydberg states, the presence of a single photon in the system creates a range (called blockade radius) around the excitation where light propagation is forbidden. This means that two Rydberg excitations cannot travel within the blockade radius, or in other words they repel each other. The nonlinear interaction not only allows one to build photonic quantum information processors [246–252], but also provides an opportunity to investigate many-body physics. One can for example create a BEC of stationary-light dark state Rydberg polaritons [217], or observe Wigner Crystallisation of photons in cold Rydberg ensembles [253]. Experimental progress in the field has been equally impressive; there have been demonstrations of EIT in a Rydberg medium [254], a strong Rydberg blockade [255–257], and even a microwave-controlled polariton–polariton interaction [258]. It is also possible to make Rydberg polaritons in a cavity, which has been demonstrated in [259, 260]. It will be interesting to further investigate how the unique advantages of the Rydberg polaritons can be harnessed for the purpose of quantum simulations and studying nonequilibrium many-body physics.