Arrays of waveguide-coupled optical cavities that interact strongly with atoms

The interaction of light and matter is of central importance to research in quantum optics [1]. The strength of this interaction depends on the intensity of the light at the location of the atom. To enhance the interaction, photons can be trapped in a high-finesse cavity of small volume, leading to high-intensity fields inside the cavity even with small photon numbers. This arrangement forms the central paradigm of cavity quantum electrodynamics (CQED) [2–4]. Early experiments in CQED used atoms that were dropped through the cavity and hence passed briefly through the region of strong interaction [5]. With the advent of laser cooling and trapping techniques, it became possible to position atoms or ions at a desired location within the cavity and to achieve 'strong' atom–photon coupling, i.e. a coupling rate higher than both the atomic decay rate and the cavity decay rate. Phenomena such as vacuum Rabi splitting [6], photon blockade [7, 8] and single-photon sources [9] have now been observed.

So far, experiments in CQED have almost exclusively been carried out with single cavities. In recent years, however, it has become possible to fabricate microcavities in regular arrays in various settings [10–17]. Due to their small size, these cavities all have small mode volume and correspondingly a strong interaction between light and matter. Indeed the 'strong coupling' regime has already been demonstrated for many of these devices [13, 18, 19].

An important next step is to find a way of coupling these cavities so that they form an array in which photons can tunnel from one cavity to another. Recently, various approaches were put forward for using such arrays of coupled cavities as a platform for quantum emulators [20–24]. These include a scheme for emulating the Bose–Hubbard Hamiltonian [20, 25–27], models of interacting Jaynes–Cummings Hamiltonians [22, 23] and an effective spin Hamiltonian [28]. The phase diagrams of these models have been studied [29–32] and the existence of a glassy phase has been predicted [29].

In this paper, we describe the design for a practical device in which atoms interact strongly with high-finesse microcavities and, at the same time, photons can tunnel with low loss from one cavity to the next through an interconnecting waveguide chip. The main components of the device are shown schematically in figure 1.

Zoom In Zoom Out Reset image size

This paper is organized as follows. In section 2, we describe the design of the device in detail and specify some realistic parameters. In section 3, we calculate the spectrum of one composite cavity. We then show (section 4) that the dynamics of the device is well approximated by a Jaynes–Cummings Hamiltonian. We describe experiments on a small scale to demonstrate two types of quantum emulator using coupled cavity arrays. As a first example (section 5) we show that a lattice of interacting Jaynes–Cummings Hamiltonians could be implemented in our device and analyse its steady states in a driven dissipative regime for a two-cavity setup [33–36]. In a second example (section 6), we show that our device is also suitable for implementing effective spin Hamiltonians and we study the dynamics in a two-cavity setup. Finally, in section 7, we outline future applications that can be based on these capabilities but go beyond the two-cavity setting.

We envisage a device with quantum emitters in several separate microcavities, coupled by waveguides on a coupling chip as illustrated in figure 1. The microcavities are based on the design described in [19]. Hemispherical micro-mirrors, wet-etched into a silicon chip [15, 19] and reflection-coated, form one side of a set of Fabry–Perot microcavities. These cavities are closed by the plane reflection-coated ends of waveguides integrated on a photonic chip. The other end of each waveguide is also reflection-coated to form a second set of optical resonators. Each waveguide/microcavity pair, coupled together through the shared mirror, forms a composite cavity. Henceforth, we will sometimes simply call this a cavity, and we will speak of a microcavity or a waveguide cavity if we wish to distinguish the component parts. The microcavities are open in the transverse direction, giving access to lasers that trap and manipulate atoms at the position of maximum interaction with cavity mode. This design benefits from the intrinsic scalability of microfabrication to achieve controlled nearest-neighbour coupling of many optical cavities.

Table 1. Definitions of relevant symbols. All frequencies are angular: e.g. 2π ×  1 GHz means 6.28 × 10⁹ s⁻¹.

Definitions and symbols
Symbol	Definition	Typical values
LW, LW	Length of the cavities	LC = (10–100) μm
		LW > (10–20) mm
ri	Amplitude reflection
	coefficients of the mirrors
Ri = |ri|2	Intensity reflection	99.95% > RC = 99.9% > 99.0%
	coefficients of the mirrors	99.95% > RW = 99.0%
		RCW ≈ 98%
FSRi = πc/Li	Free spectral range	FSRC ≈ 2π × (0.2–2) THz
		FSRW ≈ 2π × (1–2) GHz
$F_{{\rm C}}=\frac {\pi }{1-\sqrt {R_{{\rm C}} R_{{\rm CW}}}}$	Cavity finesse
$F_{{\rm W}}=\frac {\pi }{1-\sqrt {R_{{\rm W}} R_{{\rm CW}}}}$	Waveguide finesse
κi = FSRi/2Fi	Resonance linewidth
wC	Microcavity mode waist	wC = 3–5 μm
wW	Waveguide mode size	Ideally matched to wC
$g_{{\rm AC}} =\sqrt {\frac {3c \lambda ^2\gamma }{\pi ^2 w_{{\rm C}}^2 L_{{\rm C}}}}$	Atom–cavity coupling	2π × (0.1–1) GHz
γ	Atom amplitude decay rate	Rb: γ = 2π × 3 MHz
λ	Wavelength	Rb: λ = 780 nm
JWW	Tunnelling rate between	0 < JWW < 2π × 2 GHz
	adjacent waveguide resonators

In the following discussion, we assume a wavelength of λ = 780 nm, which is the D2 line of rubidium-87 atoms whose amplitude decay rate is γ = 2π × 3 MHz. In a microcavity of length L_C, the maximum coupling rate g_AC between a two-level atom and the standing-wave field of one photon is

$$\begin{equation} g_{\mathrm{AC}} = \sqrt{\frac{3 c \lambda^{2} \gamma}{\pi^{2} w_{0}^{2} L_{\mathrm{C}}}}, \end{equation} \tag{ 1 }$$

where w₀ is the 1/e² radius of the Gaussian intensity profile at the waist and c is the speed of light. Microcavities of the type envisaged can have mode waists down to 2 μm, whereas the length of the microcavity L_C is in the range of 10–100 μm. Therefore, the coupling rate is of the order of 2π × 0.1–1 GHz, and hence g_AC ≫ γ. The reflectivity of each mirror is given by R_i = 1 − (T_i + A_i). Assuming that the power transmission and absorption coefficients fulfil T_i,A_i ≪ 1, the cavity field amplitude decays at a rate

$$\begin{equation} \kappa_{\mathrm{C}} = \frac{c \xi_{\mathrm{C}}}{2 L_{\mathrm{C}}}, \quad \mathrm{with} \quad \xi_{\mathrm{C}} \approx \frac{T_{\mathrm{C}} + A_{\mathrm{C}} + T_{\mathrm{CW}} + A_{\mathrm{CW}}}{2}. \end{equation} \tag{ 2 }$$

The subscript C denotes the concave microcavity mirror, while the subscript CW denotes the plane coupling mirror between the cavity and the waveguide. As we intend to fabricate the concave mirror by isotropic etching of silicon, we expect the losses due to surface roughness to be of the order of A_C ≈ 10⁻⁴ without additional polishing. The coupling mirror, on the other hand, can have losses of the order of A_CW ≈ 5 × 10⁻⁶, assuming the waveguide facet is super-polished and a dielectric mirror formed by ion-assisted deposition is used. The decay rate of the microcavity could then be made as small as κ_C ≈ 2π × 0.01 GHz for a 100 μm long cavity. This length will give enough space to permit external optical access to the atoms. However, in order to couple effectively to the waveguide resonator, it is desirable to increase the transmission of the coupling mirror so that T_CW ≫ T_C + A_C. Therefore in practice, κ_C ⩾ 2π × 0.1 GHz. This does not necessarily imply higher loss for the composite cavity since photons that go through this mirror enter the waveguide and are not necessarily lost.

Similarly, the field in the waveguide resonator decays at a rate

$$\begin{equation} \kappa_{\mathrm{W}} = \frac{c \xi_{\mathrm{W}}}{2 L_{\mathrm{W}}} \quad \mathrm{with} \quad \xi_{\mathrm{W}} \approx \frac{T_{\mathrm{W}} + A_{\mathrm{W}} + T_{\mathrm{CW}} + A_{\mathrm{CW}} + K_{\mathrm{W}}}{2}, \end{equation} \tag{ 3 }$$

where L_W is the optical length of the waveguide cavity (refractive index times physical length) and K_W accounts for the waveguide propagation loss over one round-trip. While there is little use of integrated waveguides at 780 nm, we find in the literature [37] that propagation losses can be less than 0.02 dB cm⁻¹ for wavelengths around 800 nm in polymer waveguides. This platform offers all the necessary technological components for our device. We estimate that to achieve the necessary coupling lengths and separations between the waveguides without incurring additional bend losses, the waveguide resonator will have a length of L_W > 1 cm. We will therefore conservatively consider L_W = 2 cm, which gives a fractional round-trip loss of 1.8% and a corresponding decay rate of κ_W > 2π × 10 MHz.

The coupling from the microcavity into the waveguide resonator depends on the transmission T_CW, but also on the spatial overlap of the waveguide and microcavity modes, which will need to be optimized experimentally. In principle, this can reach unity, and we estimate that at least 90% should be achievable in practice. The photon tunnelling rate between the microcavity and the waveguide can easily exceed the free spectral range of the waveguide cavity. In this case, it is appropriate to consider the eigenmodes of the composite cavity, rather than viewing the microcavity and waveguide cavity as individual devices. We calculate the spectrum of one such coupled cavity in section 3.

Each waveguide is coupled to its two nearest neighbours by short regions of evanescent field overlap to produce a tunnelling rate J_WW. The maximum value $J^{(\max )}_{{\rm WW}}$, determined by the maximum field overlap and the length of the coupling region, can be as high as 10¹⁰ s⁻¹. This overlap can be tuned using thermal phase shifters, one near each end, to move the standing wave without changing the cavity length. In this way, J_WW can be tuned in principle over the range 0 to $J^{(\max )}_{{\rm WW}}$. In practice, a small amount of crosstalk will be unavoidable because of scattering, tuning noise and linewidth effects, but we estimate that $J^{(\min )}_{{\rm WW}} < 10^{-3} \times J^{(\max )}_{{\rm WW}}$. The maximum coupling strength must therefore be chosen judiciously to ensure access to the desired range of coupling. The thermo-optic phase shifters also allow the optical length of each waveguide cavity to be tuned at a rate of up to a few wavelengths per ms. Figure 2 illustrates how two coupled cavities exhibit a normal mode splitting and shows how the evanescent-wave coupling can be adjusted by shifting the phase of one standing wave relative to the other. In particular, the coupling is switched off for a relative phase of π/2.

Zoom In Zoom Out Reset image size

As our starting point, we consider the spectrum of a single composite cavity. Several such cavities are to be concatenated in our setup to form the coupled-cavity array.

Suppose that we couple light into one of the cavities, through the top waveguide mirror as drawn in figure 1. Without the curved microcavity mirror at the bottom, the field reflected by the waveguide cavity is given by

$$\begin{equation} E_{\mathrm{r,W}} = \frac{r_{\mathrm{W}}-r_{\mathrm{CW}} \,\mathrm{e}^{2\mathrm{i}\phi_{\mathrm{W}}}}{1-r_{\mathrm{W}}r_{\mathrm{CW}} \,\mathrm{e}^{2\mathrm{i}\phi_{\mathrm{W}}}} \, E_{\mathrm{in}}, \end{equation} \tag{ 4 }$$

where E_in is the incident field, r_W and r_CW are the amplitude reflection coefficients of the two mirrors, and 2ϕ_W = 4πL_W/λ is the round-trip propagation phase. Here we have assumed lossless mirrors for simplicity. Resonance occurs when ϕ_W = mπ with m being an integer. Similarly, a travelling field E_in,C that is incident on the microcavity produces a reflected field in the waveguide given by

$$\begin{equation} \frac{E_{\mathrm{r,C}}}{E_{\mathrm{in},\mathrm{C}}} = \frac{r_{\mathrm{CW}}-r_{\mathrm{C}} \,\mathrm{e}^{2\mathrm{i}\phi_{\mathrm{C}}}}{1-r_{\mathrm{CW}}r_{\mathrm{C}} \,\mathrm{e}^{2\mathrm{i}\phi_{\mathrm{C}}}} \equiv \tilde{r}_{\mathrm{C}} \,\mathrm{e}^{\mathrm{i}\theta}, \end{equation} \tag{ 5 }$$

where 2ϕ_C = 4πL_C/λ. This ratio defines the reflection coefficient $\tilde {r}_{{\rm C}}$ and phase shift θ at the bottom end of the waveguide cavity when the microcavity mirror is in place:

$$\begin{eqnarray}&&\begin{array}{l} \displaystyle \tilde{r}_{\mathrm{C}}^2 = \frac{r_{\mathrm{CW}}^2+r_{\mathrm{C}}^2-2r_{\mathrm{CW}}r_{\mathrm{C}}\cos(2\phi_{\mathrm{C}})}{1+r_{\mathrm{CW}}^2r_{\mathrm{C}}^2-2r_{\mathrm{CW}}r_{\mathrm{C}}\cos(2\phi_{\mathrm{C}})}\\ \displaystyle \theta = \arctan\kern-2pt\left(\frac{(r_{\mathrm{CW}}^2-1)r_{\mathrm{C}}\sin(2\phi_{\rm C})}{r_{\mathrm{CW}}(r_{\mathrm{C}}^2+1)-r_{\mathrm{C}}(r_{\mathrm{CW}}^2+1)\cos(2\phi_{\mathrm{C}})}\right).\end{array} \end{eqnarray} \tag{ 6 }$$

On replacing r_CW in equation (4) by $\tilde {r}_{{\rm C}} \,{\rm e}^{{\rm i}\theta }$, we obtain the reflected field from a composite cavity:

$$\begin{equation} E_{\mathrm{r,CW}} = \frac{r_{\mathrm{W}}-\tilde{r}_{\mathrm{C}}\,\mathrm{e}^{\mathrm{i}(\theta+2\phi_{\mathrm{W}})}}{1-r_{\mathrm{W}}\tilde{r}_{\mathrm{C}}\,\mathrm{e}^{\mathrm{i}(\theta+2\phi_{\mathrm{W}})}} \, E_{\mathrm{in}}, \end{equation} \tag{ 7 }$$

which is in resonance when

$$\begin{equation} 2 \phi_{\mathrm{W}} + \theta = 2 \pi n, \end{equation} \tag{ 8 }$$

for any integer n. This set of equations describes the system fully. It cannot generally be solved analytically, but we can determine when the reflected intensity becomes zero. This is of interest because the spectrum is generally sensitive to the presence of atoms under this condition. Zero reflected field is achieved when the microcavity can be tuned to make $r_{w}^{2} = \tilde {r}_{{\rm C}}^{2}$. Such a tuning is possible when the following inequality is satisfied:

$$\begin{equation} |E_{\mathrm{r,C}}^{\min}|^2 = \frac{(r_{\mathrm{C}}-r_{\mathrm{CW}})^2}{(1-r_{\mathrm{C}}r_{\mathrm{CW}})^2}\leqslant r_{\mathrm{W}}^2 \leqslant |E_{\mathrm{r,C}}^{\max}|^2=\frac{r_{\mathrm{C}}^2+r_{\mathrm{CW}}^2}{1+(r_{\mathrm{C}}r_{\mathrm{CW}})^2}, \end{equation} \tag{ 9 }$$

and occurs when the cavity round-trip phase satisfies

$$\begin{equation} 2\phi_{\mathrm{C}}({\rm opt}) = \arccos\kern-2pt{\left( \frac{r_{\mathrm{W}}^2(1+(r_{\mathrm{C}}r_{\mathrm{CW}})^2)-(r_{\mathrm{C}}^2+r_{\mathrm{CW}}^2)}{2r_{\mathrm{CW}}r_{\mathrm{C}}(r_{\mathrm{W}}^2-1)}\right)}. \end{equation} \tag{ 10 }$$

This determines the phase shift θ, and the resonance of the whole cavity is then ensured by adjusting the length of the waveguide to satisfy equation (8).

The semi-logarithmic graph on the top left of figure 3 shows how the reflected intensity varies with a change Δω in the angular frequency of the laser for a cavity with L_W = 15.6 mm and L_C = 156 μm and with mirror reflectivities of R_C = 0.999, R_CW = 0.98 and R_W = 0.99. Without any coupling, both constituent cavities are resonant at Δω = 0. In the reflection spectrum, we see a series of dips, symmetrically disposed around Δω = 0. The two central dips represent normal modes of the composite cavity. These would be degenerate in the absence of coupling but are split apart here. Further away from this doublet, the resonances are not so strongly coupled to the microcavity and their spacing approaches the free spectral range of the bare waveguide cavity. The horizontal line at 0.99 shows the reflection coefficient r²_W of the input mirror, while the dashed line indicates the value of $\tilde {r}_{{\rm C}}^2$. These are most nearly equal at the second reflection dip on each side, making those dips the strongest ones. The width of each resonance can be understood in a simple way because the waveguide resonator is very much longer than the microcavity. Consequently, $\tilde {r}_{{\rm C}}$ and θ are essentially constant over any given resonance line of the coupled system. By analogy with equation (2) this gives the cavity damping rate as

$$\begin{equation} \tilde{\kappa} \approx \frac{c}{2 L_{\mathrm{W}}}\xi, \quad \mathrm{with} \; \xi \approx 1 - r_{\mathrm{W}} \tilde{r}_{\mathrm{C}}. \end{equation} \tag{ 11 }$$

Zoom In Zoom Out Reset image size

The top right graph in figure 3 explores how the reflected intensity near one of the central resonances depends on the detuning of the two constituent cavities. We see that the reflection goes through a minimum whose position, indicated by the small circle, is as expected from equation (10). The two lower graphs show how the circulating power inside the waveguide resonator (left) and the microcavity (right) varies with the same detuning of the constituent cavities. In particular, we see that the condition of minimum reflected power (again indicated by a circle) corresponds closely to having maximum power inside the microcavity, as required for good atom–cavity coupling. At this point, the energy density in the microcavity is ten times that in the waveguide (although the total energy is ten times lower because of the 100-fold disparity in lengths).

In this section, we describe a series of approximations that allow us to obtain a description of our system in terms of atom–photon interactions of Jaynes–Cummings form in each cavity together with the tunnelling of photons between adjacent cavities.

4.1. Composite cavity Hamiltonian

The Hamiltonian for the field in one composite cavity is

$$\begin{equation} H_{\mathrm{cav}} = \sum_{\alpha} \omega_{\alpha} a_{j,\alpha}^{\dagger} a_{j,\alpha}, \end{equation} \tag{ 12 }$$

where the index j labels the particular cavity and index α labels its eigenmodes. The frequencies ω_α are those of the eigenmodes discussed in section 3.

4.2. Photon tunnelling between adjacent waveguides

The rate J_WW for photons to tunnel between two adjacent waveguides can be tuned over a wide range. Let us choose to make it small compared with the waveguide free spectral range, which also ensures that |J_WW| ≪ |ω_α − ω_α'|. In this regime, we can write the coupling as a tunnelling term between resonant modes of the two composite cavities,

$$\begin{equation*} H_{\mathrm{WW}} = - \sum_{\alpha, \beta} J_{\alpha, \beta} (a_{j,\alpha}^{\dagger} a_{j+1,\beta} + \mathrm{H.c.} ) \approx - \sum_{\alpha} J_{\alpha, \alpha}(a_{j,\alpha}^{\dagger} a_{j+1,\alpha} + \mathrm{H.c.}). \end{equation*} \noindent$$

When the cavity network only contains photons that are near-resonant with the normal modes α₀ and when the Rabi frequencies of atom–photon coupling are also ≪|ω_α − ω_α'|, only these modes are populated and H_WW can be simplified further to read

$$\begin{equation} H_{\mathrm{WW}} \approx - J_{\alpha_{0}, \alpha_{0}} (a_{j,\alpha_{0}}^{\dagger} a_{j+1,\alpha_{0}} + \mathrm{H.c.}). \end{equation} \noindent \tag{ 13 }$$

The tunnelling rate J_α0,α0 is related to J_WW. Yet, since J_WW is harder to determine experimentally than J_α0,α0, we do not specify the conversion here.

4.3. Atom–photon coupling

The circulating fields in the microcavity and waveguide, E_C and E_W, respectively, are given by

$$\begin{equation} E_{\mathrm{W}} = E_{\mathrm{in}} \frac{t_{\mathrm{W}}}{1-r_{\mathrm{W}} \tilde{r}_{\mathrm{C}} \,\mathrm{e}^{\mathrm{i} \left(2 \phi_{\mathrm{W}}+\theta \right)}}, \end{equation} \tag{ 14 }$$

$$\begin{equation} E_{\mathrm{C}} = E_{\mathrm{in}} \frac{t_{\mathrm{CW}} t_{\mathrm{W}} \,\mathrm{e}^{\mathrm{i} \phi_{\mathrm{W}}}}{1-r_{\mathrm{C}} r_{\mathrm{CW}} \,\mathrm{e}^{2\mathrm{i} \phi_C}+r_{\mathrm{C}} r_{\mathrm{W}} \,\mathrm{e}^{2\mathrm{i} \left(\phi_{\mathrm{C}}+\phi_{\mathrm{W}}\right)}-r_{\mathrm{CW}} r_{\mathrm{W}} \,\mathrm{e}^{2\mathrm{i} \phi_{\mathrm{W}}}}, \end{equation} \noindent \tag{ 15 }$$

where E_in is the input field from the top end of the waveguide; see figure 1. Since one photon in the composite cavity has a total energy content of ℏω, the absolute values of the circulating field amplitudes per photon in the microcavity and waveguide are

$$\begin{equation} \mathcal{E}_{\mathrm{C}} = \sqrt{\frac{\hbar \omega}{2 \varepsilon_{0} \pi w_{0}^{2}}} \sqrt{\frac{\left|E_{\mathrm{C}}\right|^{2}}{L_{\mathrm{C}} \left|E_{\mathrm{C}}\right|^{2} + L_{\mathrm{W}} \left|E_{\mathrm{W}}\right|^{2}}}, \end{equation} \tag{ 16 }$$

$$\begin{equation} \mathcal{E}_{\mathrm{W}} = \sqrt{\frac{\hbar \omega}{2 \varepsilon_{0} \pi w_{0}^{2}}} \sqrt{\frac{\left|E_{\mathrm{W}}\right|^{2}}{L_{\mathrm{C}} \left|E_{\mathrm{C}}\right|^{2} + L_{\mathrm{W}} \left|E_{\mathrm{W}}\right|^{2}}}. \end{equation} \tag{ 17 }$$

Hence the atom–photon coupling at an antinode of the microcavity is

$$\begin{equation} g = g_{\mathrm{AC}} \, \sqrt{\frac{\left|E_{\mathrm{C}}\right|^{2}}{\left|E_{\mathrm{C}}\right|^{2} + \frac{L_{\mathrm{W}}}{L_{\mathrm{C}}} \left|E_{\mathrm{W}}\right|^{2}}}, \end{equation} \tag{ 18 }$$

where g_AC is given in equation (1) for a two-level atom. This formula for a two-level atom applies because we are considering the closed cycling transition $5^{2}S_{1/2}(F=2,m_F=\pm 2) \leftrightarrow 5^{2}P_{3/2}(F'=3,m'_F=\pm 3)$ on the D₂ line of ⁸⁷Rb.

4.4. Jaynes–Cummings lattice

With the above approximations, the dynamics of atoms coupled to the composite cavity normal modes a_j,α0 can be described by a Jaynes–Cummings lattice model. Dropping the index α₀ on the field operators, the Hamiltonian for N composite cavities reads

$$\begin{eqnarray} H_{\mathrm{JCarray}} & = & \omega_{\mathrm{A}} \sum_{j=1}^{N} \sigma_{j}^{+} \sigma_{j}^{-} + \omega_{\mathrm{C}} \sum_{j=1}^{N} a_{j}^{\dagger} a_{j} \nonumber \\ && + g \sum_{j=1}^{N} (\sigma_{j}^{+} a_{j} + \sigma_{j}^{-} a_{j}^{\dagger} ) \\ && - J \sum_{j=1}^{N-1} (a_{j}^{\dagger} a_{j+1} + a_{j} a_{j+1}^{\dagger}) \nonumber. \end{eqnarray} \tag{ 19 }$$

Here, ω_C = ω_α0, σ⁻_j = |g_j〉〈e_j| is the transition operator between the excited state |e_j〉 and the ground state |g_j〉 of the atom in cavity j and a_j is the annihilation operator for photons in mode α₀ of that cavity. The losses for the system described by this Hamiltonian arise through spontaneous emission from the excited states of the atoms at rate γ and photon loss from the normal modes a_j at rate $\tilde {\kappa }$ given by equation (11).

We now consider an array of composite cavities connected to each other by nearest-neighbour coupling and excited at the end of one waveguide by a laser tuned to the normal mode a₁. The resonant pumping is described by an additional term $H_{{\rm D}} = \frac {\eta }{2} a_{1}^{\dagger } + \frac {\eta ^{*}}{2} a_{1}$ in the Hamiltonian, whereas photon leakage and spontaneous emission losses are taken into account by Markovian damping terms. The dynamics of the resulting driven dissipative system is then described by the master equation,

$$\begin{eqnarray} \frac{\mathrm{d}\rho}{\mathrm{d}t} &=& -\mathrm{i} \left[H,\rho\right] + \gamma \sum_{j} (2 \sigma_{j}^{-} \rho \sigma_{j}^{+} - \sigma_{j}^{+}\sigma_{j}^{-} \rho - \rho \sigma_{j}^{+}\sigma_{j}^{-}) \nonumber \\ && + \tilde{\kappa} \sum_{j} (2 a_{j} \rho a_{j}^{\dagger} - a_{j}^{\dagger} a_{j} \rho - \rho a_{j}^{\dagger} a_{j}), \end{eqnarray} \tag{ 20 }$$

where H = H_D + H_JCarray. For one cavity, a driving amplitude η gives a steady-state energy in the cavity of $\hbar \omega \langle a^{\dagger }a\rangle =\hbar \omega |\eta |^2/(2 \tilde {\kappa }^2)$. Also, the ratio of energy in the cavity to input power P is $2\left (|E_{{\rm C}}|^{2}L_{{\rm C}}+|E_{{\rm W}}|^{2}L_{{\rm W}}\right )/(|E_{{\rm in}}|^2c)$, where E_C and E_W are the circulating fields in the cavity. Hence, the power P of the driving laser is related to the driving amplitude η by

$$\begin{equation} P = \hbar \omega_{L} \, \frac{|\eta|^2}{\tilde{\kappa}^{2}} \, \frac{E_{\mathrm{in}}^2c}{2(|E_{\mathrm{C}}|^{2}L_{\mathrm{C}}+|E_{\mathrm{W}}|^{2}L_{\mathrm{W}})}. \end{equation} \tag{ 21 }$$

Below, we calculate the spectrum and photon statistics that can be observed at the output ports of the coupled waveguides [33–36]. These properties correspond to the most straightforward experiments that might be conducted using such a cavity array. As we shall see, they nonetheless reveal interesting physics including significantly entangled photon output states.

5.1. Composite cavity modes with strong coupling and high single-atom cooperativity

Since we are interested in having large atom–photon coupling and small photon losses, both from the cavities and from spontaneous emission, we wish to work with a mode that has high single-atom cooperativity, $C_1 = g^{2}/(2 \tilde {\kappa } \gamma )$, and is resonant with the atomic transition. The required combinations of microcavity and waveguide lengths are found by numerical optimization. In figure 4, the cavity resonances are indicated by vertical lines. The resonance near δω = 0, with the highest cooperativity (solid line) of C₁ = 17.3, coincides with the atomic transition frequency. The cavity parameters (L_C = 156.05 μm, L_W = 20.000 mm, R_C = 99.9%, R_CW = 98.0% and R_W = 99.8%) are very close to the parameters used in figure 3. This resonance corresponds to the strong reflection dip in figure 3 on the high-frequency side of the central doublet. The atom–photon coupling g and photon loss rate $\tilde {\kappa }$ are g = 0.3304 × g_AC = 2π × 33.04 MHz and $\tilde {\kappa } = 2 \pi \times 10.6\,$MHz.

Zoom In Zoom Out Reset image size

5.2. Steady state of a two-site Jaynes–Cummings array

Figure 5 shows the steady-state behaviour when two of these cavities are coupled together at a rate J = 0.2 × g_AC = 2π × 20 MHz. Each contains a rubidium atom and cavity number 1 is driven by a laser with η = 0.1 × g_AC = 2π × 10 MHz. The atoms are resonant with the normal modes a_j. Figure 5(a) shows the number of photons in each cavity: n₁ = 〈a^†₁a₁〉 and n₂ = 〈a^†₂a₂〉. The four resonance lines in the spectrum correspond to the four singly excited eigenstates of the two-site Jaynes–Cummings Hamiltonian, i.e. equation (19) with N = 2. Here, the resonances at δω ≈ ±2π × 44 MHz correspond to states where the excitation is more likely to be found in one of the cavity modes, whereas for the resonances at δω ≈ ±2π × 24 MHz the excitation is more likely to be in the atom than in the cavity field.

Zoom In Zoom Out Reset image size

In figure 5(b), we study the photon density correlations g⁽²⁾(1,1) = 〈a^†₁a^†₁a₁a₁〉/n²₁, g⁽²⁾(1,2) = 〈a^†₁a^†₂a₂a₁〉/n₁n₂ and g⁽²⁾(2,2) = 〈a^†₂a^†₂a₂a₂〉/n²₂. On the two outer resonances (δω ≈ ±2π × 44 MHz) we see small dips below unity, indicating photon anti-bunching in both cavities, g⁽²⁾(1,1) < 1 and g⁽²⁾(2,2) < 1, as well as anti-correlations between photons in distinct cavities, $\left |g^{(2)} (1,2)\right | <1$. This indicates that the coupling in both cavities is strong enough to generate an optical nonlinearity that converts coherent classical input light into a manifestly non-classical state of the cavity photons. For the resonances at δω ≈ ±2π × 44 MHz the state of the two cavity modes is approximately a superposition of one photon in cavity 1 with cavity 2 empty and one photon in cavity 2 with cavity 1 empty, which gives rise to the anti-correlations shown by g⁽²⁾(1,1), g⁽²⁾(2,2) and g⁽²⁾(1,2). Since cavity 2 is not directly driven, it experiences a lower intensity of incoming photons than cavity 1 and hence exhibits slightly stronger anti-bunching for the same nonlinearity. For this state of the cavity modes, the photons in the two cavities are expected to become entangled. This is confirmed by figure 5(c), which shows the entanglement between modes a₁ and a₂ as quantified by the logarithmic negativity E_N [39],

$$\begin{equation} E_{\mathrm{N}} = \log_{2} (\mathrm{Tr}\sqrt{\rho_{\mathrm{pt}}^{\dagger} \rho_{\mathrm{pt}}}), \end{equation} \tag{ 22 }$$

where ρ_pt is the partial transpose of the reduced density matrix of the two modes a₁ and a₂. For comparison, we note that Bell states, the maximally entangled states of two qubits, have E_N = 1.

On the inner resonances (δω ≈ ±2π × 24 MHz), by contrast, only photons involving both cavities are anti-correlated, whereas photons within each cavity bunch. This indicates that photons prefer to stick together in either of the cavities and the associated state contains superpositions of two or more photons in cavity 1 with cavity 2 empty and the reverse configuration. Hence photons in the output of one cavity tend to come in bunches, whereas if one cavity emits a photon the other cavity is unlikely to emit at the same time. Consequently, we also find entanglement between the photon modes for these resonances although less than for the resonances at δω ≈ ±2π × 44 MHz, which feature higher photon densities.

Close to δω = 0, the photons in the cavities show pronounced bunching. This emerges due to the nonlinear spectrum of our device, which here causes the laser drive to be detuned with respect to single-photon transitions but resonant with multi-photon transitions. In practice this bunching will, however, be hard to observe since the photon densities generated in these multi-photon processes are vanishingly small, as apparent from figure 5(a).

The findings shown in figures 5(b) and (c) clearly demonstrate the non-classical nature of the light fields generated in our device.

In this section we describe a second experiment that could be performed using such a device. We show how effective spin–spin interactions, as proposed in [28], can be implemented. Here, each cavity interacts with one rubidium atom, whose energy levels a,b,e associated with the D₂ line form the lambda structure depicted in figure 6. These levels are coupled both by the cavity photons (with couplings g_a and g_b) and by external laser fields (with angular frequencies ν_a, ν_b and Rabi frequencies Ω_a, Ω_b).

Zoom In Zoom Out Reset image size

Importantly, the transition $|a\rangle \leftrightarrow |b\rangle$ ($5^{2}S_{1/2}(F=1,m_{F}=1) \leftrightarrow 5^{2}S_{1/2}(F=2,m_{F}=1)$) is dipole forbidden and level |b〉 is thus metastable. Under conditions spelled out below, the dynamics can be constrained to the subspace formed by levels |a_b〉 and |b_C〉 of each atom and we can identify |a_j〉 with spin down, |↓_j〉 ≡ |a_j〉, and |b_j〉 with spin up, |↑_j〉 ≡ |b_j〉.

6.1. Outline of the approach

Following the arguments presented in section 3, we assume once again that the atom in the jth composite cavity only couples to one normal mode a_j (again, we skip the index α_j), for which we maximize the single-atom cooperativity as in section 5. The photons tunnel between adjacent composite cavities at a rate J. The Hamiltonian for this system can thus be written as

$$\begin{equation} H = H_{\mathrm{A}} + H_{\mathrm{C}} + H_{\mathrm{AC}}, \end{equation} \tag{ 23 }$$

where

$$\begin{equation} H_{\mathrm{A}} = \sum_{j=1}^N \omega_e |e_j\rangle \langle e_j| + \omega_{b} |b_j\rangle \langle b_j|, \end{equation} \tag{ 24 }$$

$$\begin{equation} H_{\mathrm{C}} = \omega_{\mathrm{C}} \sum_{j=1}^N a_j^{\dagger} a_j - J \sum_{j=1}^N (a_j^{\dagger} a_{j+1} + a_j a_{j+1}^{\dagger}), \end{equation} \tag{ 25 }$$

$$\begin{equation} H_{\mathrm{AC}} = \sum_{j=1}^N \sum_{x = a,b} \left[ \left(\frac{\Omega_x}{2} \mathrm{e}^{-\mathrm{i }\nu_x t} + g_x a_j \right) |e_j\rangle \langle x_j| + \mathrm{H.c.} \right]. \end{equation} \tag{ 26 }$$

Here, H_j describes the atoms, H₀ the cavity field and H_A the interaction between the atoms, lasers and cavity field. Also, ω_C is the a–e transition frequency and ω_AC the a–b transition frequency.

The Hamiltonian H_e can be decomposed into non-interacting collective photon modes, $H_{{\rm C}} = \sum _{k} \omega _{k} a_{k}^{\dagger } a_{k}$, where ω_b = ω_C − 2Jcos k and $a_{k} = \sqrt {\frac {2}{N+1}} \sum _{j=1}^{N} \sin (k j) a_{j}$ with $k = \frac {\pi l}{N+1}$ and l = 1,2,...,N. It is helpful to move to an interaction picture,

$$\begin{equation} H(t) = \,\mathrm{e}^{\mathrm{i} H_{0} t} (H-H_{0})\mathrm{ e}^{-\mathrm{i} H_{0} t}, \end{equation} \noindent \tag{ 27 }$$

where H is given in equation (23) and H_k reads,

$$\begin{equation} H_0 = \sum_{j=1}^N \left( \omega_e |e_j\rangle \langle e_j| + (\omega_{b} - \delta) |b_j\rangle \langle b_j|\right) + \sum_{k} \omega_{k} a_k^{\dagger} a_k, \end{equation} \tag{ 28 }$$

with the value of δ to be chosen later. In this picture, the Hamiltonian (23) becomes

$$\begin{eqnarray} &&\fl H(t) = \delta \sum_{j=1}^N |b_j\rangle \langle b_j| + \sum_{j=1}^N \sum_{x = a,b} \sum_{k} \left[ \left(\frac{\Omega_x}{2} \,\mathrm{e}^{\mathrm{i} \Delta_x t} + g_{x,j,k} \,\mathrm{e}^{\mathrm{i} \delta_{k}^{x} t} a_{k} \right) |e_j\rangle \langle x_j| + \mathrm{H.c.} \right], \end{eqnarray} \noindent \tag{ 29 }$$

where Δ_C = ω₀ − ν_a, Δ_e = ω_a − (ω_b − δ) − ν_e, as illustrated in figure 6, δ^a_b = ω_b − ω_k and δ^b_e = ω_k − (ω_k − δ) − ω_e. The coupling constants g_b and g_k are related to the couplings g_a,j,k, respectively, g_b,j,k via $g_{x,j,k} = \sqrt {\frac {2}{N+1}} \sin (k j) g_{x}$ for x = a,b.

A judicious choice of rotating frame is $\delta =\omega _{b}-\frac {1}{2}(\nu _{a}-\nu _{b})$, which ensures that δ^a_a − Δ_b = δ^b_k − Δ_b for all k and allows a convenient separation of fast rotating terms from near-resonant terms in what follows. In this rotating frame, a second-order adiabatic elimination of the excited levels |e_k〉 and the photons a_a, see appendix A and [40], yields the effective spin- $\frac {1}{2}$ Hamiltonian

$$\begin{equation} H_{\mathrm{spin}} = B \sum_{j=1}^{N} \sigma_{j}^{z} + \sum_{j \neq l} (J_{j,l} \sigma_{j}^{+} \sigma_{l}^{+} + J_{j,l}^{\star} \sigma_{j}^{-} \sigma_{l}^{-} + K_{j,l} \sigma_{j}^{+} \sigma_{l}^{-}), \end{equation} \tag{ 30 }$$

where the effective transverse field B reads

$$\begin{eqnarray*} &&\fl B = \frac{\delta}{2} - \frac{|\Omega_{b}|^{2}}{8 \Delta_{b}^{2}} \left[ \Delta_{b} - \frac{|\Omega_{b}|^{2}}{2 \Delta_{b}} - \frac{|\Omega_{a}|^{2}}{4 (\Delta_{a} - \Delta_{b})} - \sum_{k} \left( \frac{|g_{b,j,k}|^{2}}{\delta_{k}^{b} - \Delta_{b}} + \frac{|g_{a,j,k}|^{2}}{\delta_{k}^{a} - \Delta_{b}} \right) \right] \\ &&+ \frac{|\Omega_{a}|^{2}}{8 \Delta_{a}^{2}} \left[ \Delta_{a} - \frac{|\Omega_{a}|^{2}}{2 \Delta_{a}} - \frac{|\Omega_{b}|^{2}}{4 (\Delta_{b} - \Delta_{a})} - \sum_{k} \left( \frac{|g_{a,j,k}|^{2}}{\delta_{k}^{a} - \Delta_{a}} + \frac{|g_{b,j,k}|^{2}}{\delta_{k}^{b} - \Delta_{a}} \right) \right], \end{eqnarray*}$$

and the coupling constants are

$$\begin{eqnarray*} &&\fl J_{j,l} = \frac{\Omega_{a} \Omega_{b}^{\star}}{4 \Delta_{a} \Delta_{b}} \sum_{k} \frac{g_{b,j,k}^{\star} g_{a,l,k}}{\delta_{k}^{b} - \Delta_{a}},\quad K_{j,l} = \frac{|\Omega_{a}|^{2}}{4 \Delta_{a}^{2}} \sum_{k} \frac{g_{b,l,k}^{\star} g_{b,j,k}}{\delta_{k}^{b} - \Delta_{a}} + \frac{|\Omega_{b}|^{2}}{4 \Delta_{b}^{2}} \sum_{k} \frac{g_{a,j,k}^{\star} g_{a,l,k}}{\delta_{k}^{a} - \Delta_{b}}. \end{eqnarray*}$$

The Hamiltonian (30) is an appropriate description provided that |Δ_j − δ_j| ≪ FSR_x for x,y = a,b, $|\frac {\Omega _{x}}{2 \Delta _{x}}| \ll 1$ for x = a,b and $|\frac {\Omega _{x} g_{y,j,k}}{2 \Delta _{x}(\Delta _{x} - \delta _{k}^{y})}| \ll 1$ for x,y = a,b, see conditions (A.5). Furthermore, all the eigenmodes of the composite cavity should be sufficiently detuned from the |e〉 → |a〉 and |e〉 → |b〉 transitions to avoid Purcell-enhanced atomic relaxation.

The sums over all photon modes, $\sum _{k}$, depend on the number of cavities in the array. As a starting point for experiments, we focus here on the N = 2 case. The explicit expressions for this case are given in appendix B.

6.2. Effective spin dynamics for two coupled cavities

As an example of dynamical evolution under this effective spin Hamiltonian, we consider once again the two coupled cavities described in section 5.2, each containing an ⁸⁷Rb atom. The atom in cavity number 1 is prepared in state |b_y〉, whereas the other atom is placed in state |a_W〉 to form an initial effective spin state |↑₁,↓₂〉. The subsequent evolution of this state under the Hamiltonian of equation (30) is illustrated by the dotted and dash-dotted lines in figure 7, which plot the probabilities of the spin-up and spin-down states in cavity 1: $P_{\uparrow ,1} = {\rm Tr}(\mid \uparrow _{1}\rangle \langle \uparrow _{1}\mid \tilde {\rho })$ and $P_{\downarrow ,1} = {\rm Tr}(\mid \downarrow _{1}\rangle \langle \downarrow _{1}\mid \tilde {\rho })$, $\tilde {\rho }$ being the state of the effective spin system. We see that the spin oscillates in this case with a period of approximately 70 μs, determined by the coupling constant K in equation (30). The symmetry of the problem ensures that P₁ = P₂ and P_↑,2 = P_↓,1. The values of g_↓,2 and g_↑,1 are derived from equation (18) (via equations (1), (14) and (15)), together with the appropriate weightings ($1/\sqrt {2}$ and $1/\sqrt {6}$) relative to the cycling transition. The values of these and the other relevant constants are listed in the caption of figure 7.

Zoom In Zoom Out Reset image size

In deriving equation (30) we adiabatically eliminated the excited state |e〉 and therefore ignored the spontaneous emission from the atoms. In reality, this emission dephases the effective spin and causes loss of probability from the three-level system {a,b,e}. In addition, equation (30) ignores the loss of photons from the cavity. If a real experiment is to emulate the dynamics of a spin chain, as given by the Hamiltonian in equation (30), these rates of dissipation must be small enough. To test whether this is the case, we have calculated the full dynamics of the relevant atomic levels, including spontaneous emission and photon leakage using the master equation given in appendix C. The solid and dashed lines in figure 7 show the probabilities for the atom in cavity 1 to be in states a and b: $P_{a_{1}} = {\rm Tr}(|a_{1}\rangle \langle a_{1}|\tilde {\rho })$ and $P_{b_{1}} = {\rm Tr}(|b_{1}\rangle \langle b_{1}|\tilde {\rho })$. These should correspond to the spin-up and -down states of the equivalent spin model, but we see that they do not because the coherence is damped, leaving an incoherent mixture of states |a〉 and |b〉.

The situation can be rectified, as illustrated in figure 8(a). Here, the reflectivity of the microcavity mirror has been increased from 99.9 to 99.99% (together with a slight shortening of the microcavity and a minor improvement in the waveguide mirror). This improvement in the microcavity mirror allows the real system to provide a reasonable approximation to the ideal spin evolution, albeit with some residual damping. Figure 8(b) shows the excited state population P₁ on site number 1, which determines the lifetime $(P_{e}\times \frac {1}{3}\times 2\gamma )^{-1}\simeq 260\,\mu$s for atoms to be lost from the lambda system by spontaneous decay to the state 5²S(F = 2,m₂ = 2). Also shown is the small population n_e,1 of cavity photons. Significant further improvement in the coherence is possible in principle if the dissipation in the waveguide cavity is also reduced. In practice, this will require an advance in integrated waveguide technology to achieve smaller propagation losses than the current state of the art.

Zoom In Zoom Out Reset image size

Following an initial demonstration of quantum coherent coupling between two neighbouring cavities in such a system, a variety of more sophisticated experiments to ascertain coherent interactions across longer chains become accessible. A stringent test of the coherence properties of transport in extended chains of cavities is provided by the transport of quantum entanglement through such a chain. Most promising in this context appears to be the generation of entanglement between an atomic qubit and the photonic degree of freedom of the cavity in which the atom resides. Once created, the coherent coupling between constituents of the cavity array will lead to transport of the state of the photonic degree of freedom as in a harmonic chain [42]. Alternatively, one may use the natural dynamics of a chain of interacting quantum systems following a sudden quench to generate highly entangled states. Significant two-particle entanglement will, for example, build up between sites whose distance is proportional to twice the propagation speed in the chain [46]. A basic demonstration of such an experiment is already possible with three sites but may be scaled to larger arrays as the strength of the entanglement decreases only slowly with distance. The presence of entanglement between distant sites may then be verified on the basis of the measurement of simple two-site correlations [41]. An alternative that avoids the verification of entanglement would exploit the quantum coherent effect of dynamical localization [43, 44] in a harmonic chain, induced by sinusoidal modulation of the cavity resonance frequencies. Depending on the frequency of this modulation the excitation transport through the chain will exhibit sharp resonances whose depth can be used to infer the level of quantum coherence in the system without the need for tomography or entanglement verification [45].

Our proposed cavity quantum emulator could also be used to explore quantum many-body physics and to exploit it for the generation of long-distance entanglement. Indeed, systems with finite correlation length, such as the one-dimensional (1D) Heisenberg and XX models, allow sizeable ground-state end-to-end entanglement, independent of the size of the system, provided that simple patterns of site-dependent couplings are selected [48, 49]. Already a chain of four cavities together with local control of interactions would suffice to demonstrate this effect.

In summary, we have presented a practical way to realize an array of coupled cavities that could be used for quantum emulation. The device consists of open Fabry–Perot microcavities, coupled together by a waveguide chip. We have demonstrated that, under suitable conditions, the dynamical evolution is well approximated by a Jaynes–Cummings lattice Hamiltonian and that it is suitable for implementing effective spin Hamiltonians. The quantum emulator discussed here provides a very promising experimental platform for achieving controlled dynamics using photons and quantum emitters. Moreover, it has the potential to go beyond simple demonstrations to achieve significant computational power based on long-range entanglement.

We gratefully acknowledge financial support from the European Union (HIP), the Austrian Nano-initiative (PLATON-NAP), the FWF, the UK EPSRC and the Royal Society. MJH acknowledges financial support from the German Research Foundation (DFG) via the Emmy Noether project HA 5593/1-1 and the SFB 631. MBP was supported by the Alexander von Humboldt Foundation and EAH by the Royal Society.

The Schrödinger equation containing the Hamitonian H as in equation (29) reads

$$\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t} |\Psi,t \rangle = -\mathrm{ i} H(t) |\Psi,t \rangle. \end{equation} \noindent \tag{ A.1 }$$

This equation can formally be integrated to yield

$$\begin{equation} |\Psi,t + T \rangle = |\Psi,t\rangle - \mathrm{i} \int_{t}^{t+T} \mathrm{d}s H(s) |\Psi,s \rangle. \end{equation} \noindent \tag{ A.2 }$$

Iterating the right-hand side results in

$$\begin{eqnarray} &&\fl |\Psi,t + T \rangle = |\Psi,t \rangle + \sum_{n=1}^{\infty} (-\mathrm{ i})^{n} \int_t^{t+T} \mathrm{d}t_{n} H(t_{n}) \int_{t}^{t_{n}} \mathrm{d}t_{n-1} H(t_{n-1}) \cdots \int_{t}^{t_{1}} \mathrm{d}t_{0} H(t_{0})|\Psi,t \rangle. \end{eqnarray} \tag{ A.3 }$$

An effective Hamiltonian that accurately describes processes which happen on a time scale T can now be found by performing the time integrations and identifying the dominant terms. For our derivation, we assume that at time t, the excited states of the atoms are not occupied and that no photons are present,

$$\begin{equation} \langle e | \Psi,t \rangle = 0 \quad \mathrm{and} \quad a_{k} | \Psi,t \rangle = 0. \end{equation} \noindent \tag{ A.4 }$$

Furthermore, we assume for the parameters in H(t)

$$\begin{equation} \left|\frac{\Omega_{x}}{2 \Delta_{x}}\right| \ll 1 \quad \mathrm{for} \; x = a,b, \end{equation} \tag{ A.5 }$$

$$\begin{equation} \left|\frac{\Omega_{x} g_{y,j,k}}{2 \Delta_{x}(\Delta_{x} - \delta_{k}^{y})}\right| \ll 1 \quad \mathrm{for} \; x,y = a,b. \end{equation} \tag{ A.6 }$$

We keep terms up to n = 3 on the right-hand side of equation (A.3) and, by virtue of equation (A.5), neglect all oscillating terms to arrive at

$$\begin{equation} |\Psi,t + T \rangle = \left(1 + \sum_{\mu=1}^{4}\frac{\left(-\mathrm{ i} H_{0} T\right)^{\mu}}{\mu!} + \sum_{\nu=1}^{2}\frac{\left(-\mathrm{ i} H_{1} T\right)^{\nu}}{\nu!} - \mathrm{i} H_{2} T \right) |\Psi,t\rangle , \end{equation} \noindent \tag{ A.7 }$$

where

$$\begin{eqnarray*} H_{0} & = & \delta \sum_{j=1}^{N} |b_{j}\rangle \langle b_{j}|, \\ H_{1} & = & - \frac{|\Omega_{b}|^{2}}{4 \Delta_{b}} \sum_{j=1}^{N} |b_{j}\rangle \langle b_{j}| - \frac{|\Omega_{a}|^{2}}{4 \Delta_{a}} \sum_{j=1}^{N} |a_{j}\rangle \langle a_{j}| , \\ H_{2} & = & \sum_{j=1}^{N} \frac{|\Omega_{b}|^{2}}{4 \Delta_{b}^{2}} \left[ \frac{|\Omega_{b}|^{2}}{2 \Delta_{b}} + \frac{|\Omega_{a}|^{2}}{4 (\Delta_{a} - \Delta_{b})} + \sum_{k} \left( \frac{|g_{b,j,k}|^{2}}{\delta_{k}^{b} - \Delta_{b}} + \frac{|g_{a,j,k}|^{2}}{\delta_{k}^{a} - \Delta_{b}} \right) \right] |b_{j}\rangle \langle b_{j}| \nonumber \\ && + \sum_{j=1}^{N} \frac{|\Omega_{a}|^{2}}{4 \Delta_{a}^{2}} \left[ \frac{|\Omega_{a}|^{2}}{2 \Delta_{a}} + \frac{|\Omega_{b}|^{2}}{4 (\Delta_{b} - \Delta_{a})} + \sum_{k} \left( \frac{|g_{a,j,k}|^{2}}{\delta_{k}^{a} - \Delta_{a}} + \frac{|g_{b,j,k}|^{2}}{\delta_{k}^{b} - \Delta_{a}} \right) \right] |a_{j}\rangle \langle a_{j}| \nonumber \\ && + \sum_{j \neq l} \left( \frac{\Omega_{a} \Omega_{b}^{\star}}{4 \Delta_{a} \Delta_{b}} \sum_{k} \frac{g_{b,j,k}^{\star} g_{a,l,k}}{\delta_{k}^{b} - \Delta_{a}} \, |b_{j}\rangle \langle a_{j}| \otimes |b_{l}\rangle \langle a_{l}| + \mathrm{H.c.} \right) \nonumber \\ && + \sum_{j \neq l} \left( \frac{|\Omega_{b}|^{2}}{4 \Delta_{b}^{2}} \sum_{k} \frac{g_{a,j,k}^{\star} g_{a,l,k}}{\delta_{k}^{a} - \Delta_{b}} \, |b_{j}\rangle \langle a_{j}| \otimes |a_{l}\rangle \langle b_{l}| + \mathrm{H.c.} \right). \nonumber \end{eqnarray*} \noindent$$

Provided we can choose the time scale T such that ||H_{FT|| ≪ 1, ||H}T|| ≪ 1 and ||H₁T|| ≪ 1 but Δ₀T ≫ 1 and δ^x₁T ≫ 1 (x = a,b), we can write an effective Schrödinger equation,

$$\begin{equation} \frac{|\Psi,t + T \rangle - |\Psi,t\rangle}{T} \approx -\mathrm{i} \left( H_{0} + H_{1} + H_{2} \right) |\Psi,t\rangle. \end{equation} \noindent \tag{ A.8 }$$

On time scales T, the dynamics of our system is thus accurately described by the effective Hamiltonian $\mathcal {H} = H_{0} + H_{1} + H_{2}$, which is identical to the effective spin Hamiltonian of equation (30) up to an irrelevant global constant.

$$\begin{eqnarray} &&\fl B|_{N=2} = \frac{\delta}{2} - \frac{|\Omega_{b}|^{2}}{8 \Delta_{b}^{2}} \left[ \Delta_{b} - \frac{|\Omega_{b}|^{2}}{2 \Delta_{b}} - \frac{|\Omega_{a}|^{2}}{4 (\Delta_{a} - \Delta_{b})} - \frac{|g_{b}|^{2} (\delta_{b} - \Delta_{b})}{(\delta_{b} - \Delta_{b})^{2} - J^{2}} - \frac{|g_{a}|^{2} (\delta_{a} - \Delta_{b})}{(\delta_{a} - \Delta_{b})^{2} - J^{2}}\right] \nonumber \\ &&+ \frac{|\Omega_{a}|^{2}}{8 \Delta_{a}^{2}} \left[ \Delta_{a} - \frac{|\Omega_{a}|^{2}}{2 \Delta_{a}} - \frac{|\Omega_{b}|^{2}}{4 (\Delta_{b} - \Delta_{a})} - \frac{|g_{a}|^{2} (\delta_{a} - \Delta_{a})}{(\delta_{a} - \Delta_{a})^{2} - J^{2}} - \frac{|g_{b}|^{2} (\delta_{b} - \Delta_{a})}{(\delta_{b} - \Delta_{a})^{2} - J^{2}} \right],\nonumber\\ \end{eqnarray} \tag{ B.1 }$$

$$\begin{equation} J_{j,l}|_{N=2} = \frac{\Omega_{a} \Omega_{b}^{\star}}{4 \Delta_{a} \Delta_{b}} \frac{g_{b}^{\star} g_{a} J}{(\delta_{b} - \Delta_{a})^{2} - J^{2}}, \end{equation} \tag{ B.2 }$$

$$\begin{equation} K_{j,l}|_{N=2} = \frac{|\Omega_{b}|^{2}}{4 \Delta_{b}^{2}} \frac{|g_{a}|^{2} J}{(\delta_{a} - \Delta_{b})^{2} - J^{2}} + \frac{|\Omega_{a}|^{2}}{4 \Delta_{a}^{2}} \frac{|g_{b}|^{2} J}{(\delta_{b} - \Delta_{a})^{2} - J^{2}}. \end{equation} \tag{ B.3 }$$

Here, δ₂ = ω_x − ω_k and δ_a = ω_e − (ω_C − δ) − ω_b.

The three levels of ⁸⁷Rb depicted in figure 6 are two of the 5s²S_e ground states, |a〉 = |F = 1, m_b = 1〉, |b〉 = |F = 2,m_C = 1〉, and the 5p²P_1/2 excited state |e〉 = |F' = 2,m_F = 2〉. In order to compute the dynamical evolution of this system in the presence of light (here we consider only the σ⁺ component of the fields), one needs to include spontaneous emission on the π transition to the 5s²S_F ground state |x〉 = |F = 2, m_3/2 = 2〉, which is outside the three-level system under consideration. The dynamics is thus described by the master equation

$$\begin{eqnarray} &&\fl \frac{\mathrm{d}\rho}{\mathrm{d}t} = -\mathrm{i} \left[H,\rho\right] + \frac{\gamma}{2} \sum_{j} ( 2 \sigma_{j}^{ae} \rho \sigma_{j}^{ea} - \sigma_{j}^{ee} \rho - \rho \sigma_{j}^{ee})+ \frac{\gamma}{6} \sum_{j} ( 2 \sigma_{j}^{be} \rho \sigma_{j}^{eb} - \sigma_{j}^{ee} \rho - \rho \sigma_{j}^{ee})\nonumber\\ &&+ \frac{\gamma}{3} \sum_{j} ( 2 \sigma_{j}^{xe} \rho \sigma_{j}^{ex} - \sigma_{j}^{ee} \rho - \rho \sigma_{j}^{ee}) + \tilde{\kappa} \sum_{j} ( 2 a_{j} \rho a_{j}^{\dagger} - a_{j}^{\dagger} a_{j} \rho - \rho a_{j}^{\dagger} a_{j}), \end{eqnarray} \tag{ C.1 }$$

where σ^αβ_F' = |α_1/2〉〈β_F| and the atomic transition rates are weighted according to the squares of the respective dipole matrix elements.