Quantum statistics of light transmitted through an intracavity Rydberg medium

In an optically nonlinear atomic medium, the dispersion and absorption of a classical light beam depends on the powers of its amplitude [1]. At the quantum level, dispersive optical nonlinearities translate into effective interactions between photons. The ability to achieve such strong quantum optical nonlinearities is of prominent importance in quantum communication and computation, as it allows us to implement photonic conditional two-qubit gates. The standard Kerr dispersive non linearities obtained in noninteracting atomic ensembles, either in off-resonant two-level or resonant three-level configurations involving electromagnetically induced transparency (EIT), are too small to allow for quantum nonlinear optical manipulations. To further enhance such nonlinearities, EIT protocols were put forward in which the upper level of the ladder is a Rydberg level. In such schemes, the strong van der Waals interactions between Rydberg atoms result in a cooperative Rybderg blockade phenomenon [2–4], where each Rydberg atom prevents the excitation of its neighbours inside a 'blockade sphere'. This Rydberg blockade deeply changes the EIT profile and leads to magnified nonlinear susceptibilities [5–8]. In particular, giant dispersive nonlinear effects were experimentally obtained in an off-resonant Rydberg–EIT scheme using cold rubidium atoms placed in an optical cavity [9, 10]. In this paper, we theoretically investigate the quantum statistical properties of the light generated in the latter protocol. Note that, contrary to other theoretical works, e.g. [11, 12], here, we are interested in the dispersive region. Moreover, since we place the atoms in a cavity rather than in free-space, the theoretical framework and calculations we perform also differ from [11, 12]. In particular, a technical benefit of our approach is that we are not restricted to considering only photon pairs but could, in principle, investigate higher-order correlations.

We first write the dynamical equations for the system of interacting three-level atoms coupled to the strong control field and the nonresonant cavity mode fed by the probe beam. We show that, under some assumptions, the system effectively behaves as a large spin coupled to the cavity mode [13]. We then compute the steady-state second-order correlation function to characterize the emission of photons out of the cavity. Our numerical simulations suggest that the bunched or antibunched nature of the outgoing light, as well as its coherence time, may be controlled through adjusting the detuning between the cavity mode and probe field frequencies.

The paper is structured as follows. In section 2, we present our setup and the assumptions we make to compute its dynamics. We also explain the analytical and numerical methods we employ to calculate the second-order ${{g}^{\left( 2 \right)}}$ correlation function of the outgoing light beam. In section 3, we present and interpret the results of the simulations we obtained for ${{g}^{\left( 2 \right)}}\left( 0 \right)$ and ${{g}^{\left( 2 \right)}}\left( \tau >0 \right)$ on the specific experimental case considered in [9]. Finally, we conclude in section 4 by evoking open questions and perspectives of our work. The appendices address the supplementary technical details that are omitted in the text for readability.

The system we consider comprises N atoms that present a three-level ladder structure with a ground $\left| g \right\rangle$, intermediate $\left| e \right\rangle$, and Rydberg states $\left| r \right\rangle$ (see figure 1). The energy of the atomic level $\left| k=g,e,r \right\rangle$ is denoted by ℏω_k (by convention ${{\omega }_{g}}=0$), and the dipole decay rates from the intermediate and Rydberg states are denoted by γ_e and γ_r, respectively. The transitions $\left| g \right\rangle \leftrightarrow \left| e \right\rangle$ and $\left| e \right\rangle \leftrightarrow \left| r \right\rangle$ are, respectively, driven by a weak probe field of frequency ω_p and a strong control field of frequency ω_cf. To limit absorption, both fields are off-resonant; the respective detunings are given by ${{\Delta }_{e}}\equiv \left( {{\omega }_{p}}-{{\omega }_{e}} \right)$ and ${{\Delta }_{r}}\equiv \left( {{\omega }_{p}}+{{\omega }_{cf}}-{{\omega }_{r}} \right)$. Moreover, to enhance dispersive effects while keeping a high input–output coupling efficiency, the atoms are placed in an optical low-finesse cavity. The transition $\left| g \right\rangle \leftrightarrow \left| e \right\rangle$ is supposed to be in the neighbourhood of a cavity resonance. The frequency and annihilation operator of the corresponding mode are denoted by ω_c and a, respectively; the detuning of this mode with the probe laser is defined by ${{\Delta }_{c}}\equiv \left( {{\omega }_{p}}-{{\omega }_{c}} \right)$, and α denotes the feeding rate of the cavity mode with the probe field, which is assumed real for simplicity. Finally, we introduce g and Ω_cf, which are the single-atom coupling constant of the transition $\left| g \right\rangle \leftrightarrow \left| e \right\rangle$ with the cavity mode and the Rabi frequency of the control field on the transition $\left| e \right\rangle \leftrightarrow \left| r \right\rangle$, respectively. In the following paragraphs, we study the dynamics of the system, which, under some assumptions, is equivalent to a damped harmonic oscillator, i.e. the cavity mode, coupled to an assembly of spins $\frac{1}{2}$, and the Rydberg bubbles corresponding to the 'super-atoms' delimited by the Rydberg blockade spheres.

Zoom In Zoom Out Reset image size

Starting from the full Hamiltonian, we perform the rotating wave approximation and adiabatically eliminate the intermediate state $\left| e \right\rangle$ as described in appendix A. Note that the result we obtain coincides with the lowest order of the EIT model—the nonlinearity of the three-level atoms is neglected, and the leading nonlinear effect comes from the Rydberg–Rydberg collisional effects. The system therefore consists of N effective two-level atoms $\left\{ \left| g \right\rangle ,\left| r \right\rangle \right\}$, with an effective power-broadened dipole decay rate from the Rydberg level

$$\begin{eqnarray*}&&{{\tilde{\gamma }\,}_{r}}=\left( {{\gamma }_{r}}+\frac{\Omega _{cf}^{2}{{\gamma }_{e}}}{4\left( \Delta _{e}^{2}+\gamma _{e}^{2} \right)} \right),\end{eqnarray*}$$

coupled to the cavity mode of the effective decay rate

$$\begin{eqnarray*}&&{{\tilde{\gamma }\,}_{c}}=\left( {{\gamma }_{c}}+\frac{{{g}^{2}}N{{\gamma }_{e}}}{\Delta _{e}^{2}+\gamma _{e}^{2}} \right)\end{eqnarray*}$$

increased by the coupling to the atomic ensemble. The Hamiltonian reads

$$\begin{eqnarray*} \tilde{H}\, & = & -\hbar {{\tilde{\Delta }\,}_{r}}\left( \underset{n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\sigma _{rr}^{\left( n \right)} \right)+\underset{m<n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\hbar {{\kappa }_{mn}}\sigma _{rr}^{\left( m \right)}\sigma _{rr}^{\left( n \right)} \\ {} & {} & -\hbar {{\tilde{\Delta }\,}_{c}}{{a}^{\dagger }}a+\hbar \alpha \left( a+{{a}^{\dagger }} \right)+\hbar {{g}_{\text{eff}}}\left\{ a\left( \underset{n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\sigma _{rg}^{\left( n \right)} \right)+h.c. \right\} \\ \end{eqnarray*}$$

In this expression, we have introduced the atomic operators $\sigma _{kl}^{\left( n \right)}\equiv {{\mathbb{I}}^{\left( 1 \right)}}\otimes \ldots \otimes {{\mathbb{I}}^{\left( n-1 \right)}}\otimes \left| k \right\rangle \left\langle l \right|\otimes {{\mathbb{I}}^{\left( n+1 \right)}}\otimes \ldots \otimes {{\mathbb{I}}^{\left( N \right)}}$ for $\left( k,l \right)=g,e,r$ as well as the effective detunings

$$\begin{eqnarray*}&&{{\tilde{\Delta }\,}_{r}}\equiv {{\Delta }_{r}}-\frac{\Omega _{cf}^{2}{{\Delta }_{e}}}{4\left( \Delta _{e}^{2}+\gamma _{e}^{2} \right)}\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{{\tilde{\Delta }\,}_{c}}\equiv {{\Delta }_{c}}-\frac{{{g}^{2}}N{{\Delta }_{e}}}{\Delta _{e}^{2}+\gamma _{e}^{2}},\end{eqnarray*}$$

respectively shifted from Δ_r and Δ_c by the light-induced Stark shift of the control beam and by the linear atomic susceptibility. The quantity ${{\kappa }_{mn}}\equiv {{C}_{6}}/\|{{{\vec{r}}_{m}}-{{\vec{r}}_{n}}}\|^{6}$ is the van der Waals interaction between atoms $\left( m,n \right)$ in their Rydberg level—when atoms are in the ground or intermediate states, their interactions are neglected—while

$$\begin{eqnarray*}&&{{g}_{\text{eff}}}=\frac{g{{\Omega }_{cf}}}{2{{\Delta }_{e}}}\end{eqnarray*}$$

is the effective coupling strength of the two-photon transition $\left| g \right\rangle \to \left| r \right\rangle$ driven by the cavity mode and the control laser.

At this point, following [13], we introduce the Rydberg bubble approximation. In this approach, the strong Rydberg interactions are assumed to effectively split the sample into ${{\mathcal{N}}_{b}}$ bubbles $\left\{ {{\mathcal{B}}_{\alpha =1,\ldots ,{{\mathcal{N}}_{b}}}} \right\}$, each of which contains ${{n}_{b}}=\left( \frac{N}{{{\mathcal{N}}_{b}}} \right)$ atoms but can only accomodate a single Rydberg excitation, delocalized over the bubble. Note that, within this approximation, all bubbles have the same radius, which is fixed by the coefficient C₆ and the detunings Δ_r and Δ_e, as well as the control field Rabi frequency Ω_cf; assuming that the atomic ensemble is homogeneous, the number of atoms per bubble n_b is approximately given by [9]

$$\begin{eqnarray*}&&{{n}_{b}}=\frac{2{{\pi }^{2}}{{\rho }_{\text{at}}}}{3}\sqrt{\frac{\left| {{C}_{6}} \right|}{{{\Delta }_{r}}-\Omega _{cf}^{2}/(4{{\Delta }_{e}})}}\end{eqnarray*}$$

where ${{\rho }_{\text{at}}}$ is the atomic density. Each bubble can, therefore, be viewed as an effective spin $\frac{1}{2}$ whose Hilbert space is spanned by

$$\begin{eqnarray*} \left| {{-}_{\alpha }} \right\rangle & = & \left| {{G}_{\alpha }} \right\rangle \equiv \underset{{{i}_{\alpha }}\in {{\mathcal{B}}_{\alpha }}}{\mathop{\otimes }}\,\left| {{g}_{{{i}_{\alpha }}}} \right\rangle \\ \left| {{+}_{\alpha }} \right\rangle & = & \left| {{R}_{\alpha }} \right\rangle \equiv \frac{1}{\sqrt{{{n}_{b}}}}\left\{ \left| rg\ldots g \right\rangle +\ldots +\left| g\ldots gr \right\rangle \right\}, \\ \end{eqnarray*}$$

the ground state of the bubble ${{\mathcal{B}}_{\alpha }}$ and its symmetric singly Rydberg excited state, respectively. Introducing the bubble spin-$\frac{1}{2}$ operators $\text{s}_{-}^{\left( \alpha \right)}=\hbar \left| {{-}_{\alpha }} \right\rangle \left\langle {{+}_{\alpha }} \right|$, where the operator $\text{s}_{-}^{\left( \alpha \right)}$ corresponds to the lowering operator of the spin and the annihilation of a Rydberg excitation, one can write the Hamiltonian under the approximate form (see appendix A)

$$\begin{eqnarray*}\begin{matrix} \tilde{H}\, & \approx & -\hbar {{\tilde{\Delta }\,}_{c}}{{a}^{\dagger }}a+\hbar \alpha \left( a+{{a}^{\dagger }} \right) \\ {} & {} & -\hbar {{\tilde{\Delta }\,}_{r}}\left( \frac{{{\mathcal{N}}_{b}}}{2}+\frac{{{\text{J}}_{z}}}{\hbar } \right) \\ {} & {} & +{{g}_{\text{eff}}}\sqrt{{{n}_{b}}}\left( a{{\text{J}}_{+}}+{{a}^{\dagger }}{{\text{J}}_{-}} \right) \\ &&\end{matrix}\end{eqnarray*}$$

where we introduced the collective angular momentum ${{\text{J}}_{-}}\equiv \mathop{\sum }^{}_{\alpha =1}^{{{\mathcal{N}}_{b}}}\text{s}_{-}^{\left( \alpha \right)}$. The system is, therefore, equivalent to a large spin, i.e. the assembly of spin-$\frac{1}{2}$ Rydberg bubbles, coupled to a harmonic oscillator. Its density matrix satisfies the master equation

$$\begin{eqnarray*} {{\partial }_{t}}\tilde{\rho }\, & = & \mathcal{L}\tilde{\rho }\, \\ {} & = & \frac{1}{\text{i}\hbar }\left[ \tilde{H}\,,\tilde{\rho }\, \right]+{{\tilde{\gamma }\,}_{c}}\left\{ 2a\tilde{\rho }\,{{a}^{\dagger }}-{{a}^{\dagger }}a\tilde{\rho }\,-\tilde{\rho }\,{{a}^{\dagger }}a \right\} \\ {} & {} & +{{\tilde{\gamma }\,}_{r}}\underset{\alpha =1}{\overset{{{\mathcal{N}}_{b}}}{\mathop{\mathop{\sum }^{}}}}\,\left\{ 2\text{s}_{-}^{\left( \alpha \right)}\tilde{\rho }\,\text{s}_{+}^{\left( \alpha \right)}-\text{s}_{+}^{\left( \alpha \right)}\text{s}_{-}^{\left( \alpha \right)}\tilde{\rho }\,-\tilde{\rho }\,\text{s}_{+}^{\left( \alpha \right)}\text{s}_{-}^{\left( \alpha \right)} \right\} \\ \end{eqnarray*} \tag{ 1 }$$

One can also write the Heisenberg–Langevin equations for the time-dependent operators $a\left( t \right),{{\text{J}}_{-}}\left( t \right)$

$$\begin{eqnarray}&&{{\partial }_{t}}a=\left( \text{i}{{\tilde{\Delta }\,}_{c}}-{{\tilde{\gamma }\,}_{c}} \right)a-\text{i}\alpha +\text{i}{{g}_{eff}}\sqrt{{{n}_{b}}}\frac{{{\text{J}}_{-}}}{\hbar }+{{\tilde{a}\,}_{in}}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{\partial }_{t}}{{\text{J}}_{-}}=\left( \text{i}{{\tilde{\Delta }\,}_{r}}-{{\tilde{\gamma }\,}_{r}} \right){{\text{J}}_{-}}+\text{i}\hbar {{g}_{eff}}\sqrt{N{{\mathcal{N}}_{b}}}a+{{\tilde{\text{J}}\,}_{in}}\end{eqnarray} \tag{ 3 }$$

where ${{\tilde{a}\,}_{in}},{{\tilde{\text{J}}\,}_{in}}\equiv \mathop{\sum }^{}_{n=1}^{N}\tilde{F}\,_{gr}^{\left( n \right)}$ are the Langevin forces associated with a and ${{\text{J}}_{-}}$, respectively. Note that we neglected the effect of extra dephasing due to, for example, collisions or laser fluctuations.

To study the quantum properties of the light transmitted through the cavity, we shall compute the function $g_{\text{out}}^{\left( 2 \right)}$, which characterizes the two-photon correlations. In the input–output formalism [14], one shows that this function simply equals the function ${{g}^{\left( 2 \right)}}$ for the intracavity field (see appendix B for details) given by

$$\begin{eqnarray*}&&{{g}^{\left( 2 \right)}}\left( \tau \right)=\frac{\text{Tr}\left\{ {{a}^{\dagger }}a{{e}^{\mathcal{L}\tau }}\left[ a{{\rho }_{ss}}{{a}^{\dagger }} \right] \right\}}{\text{Tr}{{[{{a}^{\dagger }}a{{\rho }_{ss}}]}^{2}}}\end{eqnarray*}$$

where ρ_ss denotes the steady state of the system defined by $\mathcal{L}{{\rho }_{ss}}=0$; see equation (1).

In the region of the small feeding parameter α, one can compute ρ_ss numerically by propagating in time the initial state ${{\rho }_{0}}\equiv \left| {{N}_{r}}=0 \right\rangle \left\langle {{N}_{r}}=0 \right|\otimes \left| {{n}_{c}}=0 \right\rangle \left\langle {{n}_{c}}=0 \right|$, where $\left| {{N}_{r}}=0,1,\ldots ,{{\mathcal{N}}_{b}} \right\rangle$ represents the symmetric state in which ${{N}_{r}}\equiv \left( \frac{{{\mathcal{N}}_{b}}}{2}+\frac{{{\text{J}}_{z}}}{\hbar } \right)$ bubbles are excited, and $\left| {{n}_{c}}=0,1,\ldots \right\rangle$ are the Fock states of the cavity mode. To this end, one applies the Liouvillian evolution operator ${{e}^{\mathcal{L}t}}$ in a truncated basis, restricted to states of low numbers of excitations (typically with ${{n}_{c}}+{{N}_{r}}\leqslant 6$). The steady state is reached in the limit of large times,ideally when $t\to \infty$. The denominator of the ratio in equation (2) is directly obtained from ρ_ss. To compute its numerator, one first propagates in time $a{{\rho }_{ss}}{{a}^{\dagger }}$ from t = 0 to τ, using the same procedure as earlier, then applies the operator ${{a}^{\dagger }}a$ and takes the trace.

In the region of weak feeding, it is also possible to get a perturbative expression for ${{g}^{\left( 2 \right)}}\left( 0 \right)$ by computing the expansion of $\left\langle {{a}^{\dagger }}{{a}^{\dagger }}a{{a}}\right\rangle _{ss}$ and $\left\langle {{a}^{\dagger }}{{a}}\right\rangle _{ss}$ in powers of α. To this end, one uses the Heisenberg equations of the system in equations (2) and (3) to derive the hierarchy of equations relating the different mean values and correlations $\left\langle {\ldots }\right\rangle _{ss}$ up to the fourth order in α. After straightforward, though lengthy, algebra, one gets an expression for ${{g}^{\left( 2 \right)}}\left( 0 \right)$ that is too cumbersome to be reproduced here, but allows for faster calculations than the numerical approach. Such a fully analytical treatment, however, cannot, to our knowledge, be extended to ${{g}^{\left( 2 \right)}}\left( \tau >0 \right)$; for $\tau >0$, therefore, we rely entirely on numerical simulations.

To conclude this section, we consider the region of a large number of bubbles and low number of excitations, i.e. ${{\mathcal{N}}_{b}}\gg 1$ and $\frac{{{\text{J}}_{z}}}{\hbar }\ll {{\mathcal{N}}_{b}}$. As shown in appendix A, the operator $b\equiv \frac{{{J}_{-}}}{\hbar \sqrt{{{\mathcal{N}}_{b}}}}$ is then approximately bosonic, and the term $\left( \frac{{{\mathcal{N}}_{b}}}{2}+\frac{{{\text{J}}_{z}}}{\hbar } \right)$ can be put under the form

$$\begin{eqnarray*}\begin{matrix} \left( \frac{{{\mathcal{N}}_{b}}}{2}+\frac{{{\text{J}}_{z}}}{\hbar } \right) & \approx & \frac{{{\text{J}}_{+}}{{\text{J}}_{-}}}{{{\hbar }^{2}}\left( {{\mathcal{N}}_{b}}+1 \right)}+\frac{{{\left( {{\text{J}}_{+}}{{\text{J}}_{-}} \right)}^{2}}}{{{\hbar }^{4}}{{\left( {{\mathcal{N}}_{b}}+1 \right)}^{3}}} \\ {} & \approx & \frac{{{\mathcal{N}}_{b}}}{\left( {{\mathcal{N}}_{b}}+1 \right)}{{b}^{\dagger }}b+\frac{\mathcal{N}_{b}^{2}}{{{\left( {{\mathcal{N}}_{b}}+1 \right)}^{3}}}{{b}^{\dagger }}b{{b}^{\dagger }}b \\ {} & \approx & {{b}^{\dagger }}b+\frac{1}{{{\mathcal{N}}_{b}}}{{b}^{\dagger }}{{b}^{\dagger }}bb \\ &&\end{matrix}\end{eqnarray*}$$

Finally, we get the following approximate expression for the effective Hamiltonian

$$\begin{eqnarray*}\begin{matrix} \tilde{H}\, & \approx & -\hbar {{\tilde{\Delta }\,}_{c}}{{a}^{\dagger }}a+\hbar \alpha \left( a+{{a}^{\dagger }} \right)-\hbar {{\tilde{\Delta }\,}_{r}}{{b}^{\dagger }}b-\frac{\hbar \bar{\kappa }}{2}{{b}^{\dagger }}{{b}^{\dagger }}bb+\hbar {{g}_{\text{eff}}}\sqrt{N}\left( a{{b}^{\dagger }}+{{a}^{\dagger }}b \right) \\ &&\end{matrix}\end{eqnarray*}$$

where $\bar{\kappa }\equiv 2{{\tilde{\Delta }\,}_{r}}/{{\mathcal{N}}_{b}}$. In this region, the system thus behaves as two coupled oscillators: one is harmonic, the cavity field; the other is anharmonic, the Rydberg bubble field. The cavity resonance is therefore shifted depending on the number of excitations in the system.

In the following section, we present and discuss the results we obtained with the specific system used in [9]. It appears that one can choose the bunched or antibunched behaviour of the light transmitted through the cavity by adjusting the detuning Δ_c. We also show that the time behaviour of the function ${{g}^{\left( 2 \right)}}\left( \tau \right)$ depends on the region considered, and can be roughly understood as resulting from the damped exchange of a single excitation between the atoms and field.

We now consider the physical setup presented in [9], i.e. an ensemble of $^{87}\text{Rb}$ atoms whose state space is restricted to the levels $\left| g \right\rangle =\left| 5{{s}_{\frac{1}{2}}};F=2 \right\rangle$, $\left| e \right\rangle =\left| 5{{p}_{\frac{3}{2}}};F=3 \right\rangle$, and $\left| r \right\rangle =\left| 95{{d}_{\frac{5}{2}}};F=4 \right\rangle$ with the decay rates ${{\gamma }_{e}}=2\pi \times 3$ MHz and ${{\gamma }_{r}}=2\pi \times 0.03$ MHz. The other physical parameters must be designed so that strong nonlinearities may be observed at the single-photon level. In the specific system considered here, we find that this is achieved for a cavity decay rate ${{\gamma }_{c}}=2\pi \times 1$ MHz, a volume of the sample $V=40\pi \times 15\times 15\;\mu {{\text{m}}^{3}}$, a sample density ${{n}_{at}}=0.4\;\mu {{\text{m}}^{-3}}$, a control laser Rabi frequency ${{\Omega }_{cf}}=10{{\gamma }_{e}}$, a cooperativity C = 1000, a detuning of the intermediate level ${{\Delta }_{e}}=-35{{\gamma }_{e}}$, a detuning of the Rydberg level ${{\Delta }_{r}}=0.4{{\gamma }_{e}}$, and a cavity feeding rate $\alpha =0.01{{\gamma }_{e}}$. For these parameters, the cavity detuning $\Delta _{c}^{\left( 0 \right)}=-6.1{{\gamma }_{e}}$ corresponds to the maximal average number of photons in the cavity. Note that these physical parameters are experimentally realistic and feasible.

Let us first focus on the second-order correlation function at zero time ${{g}^{\left( 2 \right)}}\left( 0 \right)$, represented in figure 2(a) as a function of the reduced detuning $\theta \equiv \left( {{\Delta }_{c}}-\Delta _{c}^{\left( 0 \right)} \right)/{{\gamma }_{e}}$. The numerical and analytical results are in such good agreement for the region considered that the corresponding curves cannot be distinguished. One notes a strong bunching peak (B) ${{\theta }_{B}}=-4.9$ and a deep antibunching area centered on (A) ${{\theta }_{A}}=0$. This suggests that around (A), photons are preferably emitted one by one, while around (B) they are preferably emitted by pairs. Note, however, that as a ratio, ${{g}^{\left( 2 \right)}}\left( 0 \right)$ only gives information on the relative importance of pair and single-photon emissions. Its peaks, therefore, do not correspond to the maxima of photon pair emission, but to the best possible compromises between $\left\langle {a}^{\dagger }{a}^{\dagger }a{a}\right\rangle _{ss}$ and $\left\langle {{a}^{\dagger }}a\right\rangle _{ss}^{2}$, which can be checked by comparing figures 2(a) and (b). Hence, pair emission might dominate in a region where the number of photons coming out from the cavity is actually very small.

Zoom In Zoom Out Reset image size

We now investigate the behaviour of ${{g}^{\left( 2 \right)}}\left( \tau >0 \right)$ for two different values of the detuning, i.e. ${{\theta }_{B}}=-4.9$ and ${{\theta }_{A}}=0$, which respectively correspond to the peak (B) and minimum (A) of ${{g}^{\left( 2 \right)}}\left( 0 \right)$. The numerical simulations we obtained are given in figure 3. The plot relative to (B) exhibits damped oscillations, alternatively showing a bunched $\left( {{g}^{\left( 2 \right)}}\left( \tau \right)>1 \right)$ or antibunched $\left( {{g}^{\left( 2 \right)}}\left( \tau \right)<1 \right)$ behaviour. The plot corresponding to (A) always remains on the antibunched side, though asymptotically tending to 1.

Zoom In Zoom Out Reset image size

The dynamical features observed can be understood and satisfactorily accounted for by a simple three-level model. Indeed, due to the weakness of α, the system in its steady state is expected to contain, at most, two excitations (either photonic or atomic). After a photon detection at t = 0, it contains, at most, one excitation, which can be exchanged between the cavity field and atoms, as it has been known for a long time [15, 16]. In other words, the operator $a{{\rho }_{ss}}{{a}^{\dagger }}$ can be expanded in the space restricted to the three states $\left\{ \left| 00 \right\rangle \equiv \left| {{N}_{r}}=0,{{n}_{c}}=0 \right\rangle ,\left| 01 \right\rangle \equiv \left| {{N}_{r}}=0,{{n}_{c}}=1 \right\rangle ,\left| 10 \right\rangle \equiv \left| {{N}_{r}}=1,{{n}_{c}}=0 \right\rangle \right\}$ and the effective non-Hermitian Hamiltonian for the system, in this subspace, takes the following form:

$$\begin{eqnarray*}{{\text{H}}_{3}}=\hbar \left[ \begin{matrix} 0 & \alpha & 0 \\ \alpha & -{{\tilde{\Delta }\,}_{c}}-\text{i}{{\tilde{\gamma }\,}_{c}} & {{g}_{\text{eff}}}\sqrt{N} \\ 0 & {{g}_{\text{eff}}}\sqrt{N} & -{{\tilde{\Delta }\,}_{r}}-\text{i}{{\tilde{\gamma }\,}_{r}} \\ &&\end{matrix} \right]\end{eqnarray*}$$

The order of magnitude of the frequencies and decay rates of the oscillations observed for ${{g}^{\left( 2 \right)}}(t)$ in the specific cases (A,B) are satisfactorily recovered by this Hamiltonian, which validates the schematic model we used and suggests it comprises the main physical processes at work.

To conclude this section, it is worth mentioning that the two-boson approximation, though strictly speaking not applicable here—the parameters considered in this section indeed correspond to a number of bubbles ${{\mathcal{N}}_{b}}\simeq 2$—yields, however, the qualitative behaviour for ${{g}^{\left( 2 \right)}}\left( 0 \right)$. The minimum is correctly located, though slightly higher than in the spin model; the antibunching peak is slightly shifted towards positive detunings and is weaker than in the previous treatment. These discrepancies result from too low a value of the nonlinearity parameter $\bar{\kappa }$; they can be corrected by replacing $\bar{\kappa }=2\tilde{\Delta }\,/{{\mathcal{N}}_{b}}$ with $\bar{\kappa }\prime =2\tilde{\Delta }\,/\left( {{\mathcal{N}}_{b}}-1 \right)$ in the two-boson Hamiltonian. We first note that $\bar{\kappa }$ and $\bar{\kappa }\prime$ coincide in the region of a large number of bubbles. Moreover, $\bar{\kappa }\prime$ makes sense in the region of a low number of bubbles: in particular, when ${{\mathcal{N}}_{b}}\to 1$, i.e. when only one bubble is available, the nonlinearity proportional to $\bar{\kappa }\prime$ diverges accordingly, thus forbidding the boson field to contain more than one excitation. Finally, let us mention that $\bar{\kappa }\prime$ can also be recovered via a perturbative treatment of the full model, which will be presented in a future paper.

In this work, we studied how the strong Rydberg–Rydberg van der Waals interactions in an atomic medium may affect the quantum statistical properties of an incoming light beam. In our model, atoms are located in a low-finesse cavity and subject to a weak signal beam and a strong control field. These two fields nonresonantly drive the transition from the ground to a Rydberg level. The system was shown to effectively behave as a large spin coupled to a damped harmonic oscillator, i.e. the assembly of Rydberg bubbles and the cavity mode, respectively. The strong anharmonicity of the atomic spin affects the quantum statistics of the outgoing light beam. To demonstrate this effect, we performed analytical and numerical calculations of the second-order correlation function ${{g}^{\left( 2 \right)}}\left( \tau \geqslant 0 \right)$. The results we obtained on a specific physical example with rubidium atoms indeed show that the transmitted light presents either bunched or antibunched characters, depending on the detuning between the cavity mode and the probe field. This suggests that in such a setup, one could design light of arbitrary quantum statistics through appropriately adjusting the physical parameters.

In this work, we performed the Rydberg bubble approximation, which allowed us to derive a tractable effective Hamiltonian. This scheme is, however, questionable: interactions between bubbles are indeed neglected, and the different spatial arrangements of the bubbles in the sample are not considered. Though challenging, it would be interesting to run full simulations of the system, rejecting those states that are too far off-resonant due to Rydberg–Rydberg interactions. Besides validating the assumption of the present work, this would indeed enable us to consider other regions, such as, for instance, the case of resonant transition towards the Rydberg level. We also implicitly made the assumption that the cavity mode and control beam were homogeneous. Spatial variations should be included in the model, and their potential influence studied in future work. Finally, due to the very weak probe field region considered in this paper, we only presented results on the function ${{g}^{\left( 2 \right)}}\left( \tau \right)$: the production of $n=3,4,\ldots$ correlated photons is indeed very unlikely. In principle, we can, however, numerically compute ${{g}^{\left( n \right)}}\left( \tau \right)$ for any $n>2$, which might be relevant in a future work, if addressing stronger probe fields.

This work was supported by the EU through the ERC Advanced Grant 246669 'DELPHI' and the Collaborative Project 600645 'SIQS'.

A.1. Rotating wave approximation

The full Hamiltonian of the system can be written under the form

$$\begin{eqnarray*} H & = & {{H}_{a}}+{{H}_{c}}+{{V}_{a-c}} \\ {{H}_{a}} & = & \hbar {{\omega }_{e}}\underset{n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\sigma _{ee}^{\left( n \right)}+\hbar {{\omega }_{r}}\underset{n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\sigma _{rr}^{\left( n \right)} \\ {} & {} & +\hbar {{\Omega }_{cf}}\;\cos \left( {{\omega }_{cf}}t \right)\underset{n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\left( \sigma _{re}^{\left( n \right)}+\sigma _{er}^{\left( n \right)} \right) \\ {} & {} & +\underset{m<n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\hbar {{\kappa }_{mn}}\sigma _{rr}^{\left( m \right)}\sigma _{rr}^{\left( n \right)} \\ {{H}_{c}} & = & \hbar \left[ {{\omega }_{c}}{{a}^{\dagger }}a+2\alpha \;\cos \left( {{\omega }_{p}}t \right)\left( a+{{a}^{\dagger }} \right) \right] \\ {{V}_{a-c}} & = & \underset{n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\hbar g\left( a+{{a}^{\dagger }} \right)\left( \sigma _{eg}^{\left( n \right)}+\sigma _{ge}^{\left( n \right)} \right) \\ \end{eqnarray*}$$

where $\sigma _{\alpha \beta }^{\left( n \right)}\equiv {{\mathbb{I}}^{\left( 1 \right)}}\otimes \ldots \otimes {{\mathbb{I}}^{\left( n-1 \right)}}\otimes \left| \alpha \right\rangle \left\langle \beta \right|\otimes {{\mathbb{I}}^{\left( n+1 \right)}}\otimes \ldots \otimes {{\mathbb{I}}^{\left( N \right)}}$, $\hbar {{\omega }_{\alpha }}$ is the energy of the atomic level $\left| \alpha \right\rangle$ for $\alpha =e,r$ (with the convention ${{\omega }_{g}}=0$), and ${{\kappa }_{mn}}\equiv \frac{{{C}_{6}}}{\|{{{\vec{r}}}_{m}}-{{{\vec{r}}}}\|}_{n}^{6}$ denotes the van der Waals interaction between atoms in the Rydberg level. When atoms are in the ground or intermediate states, their interactions are neglected.

We switch to the rotating frame defined by $\left| \psi \right\rangle \to \left| \tilde{\psi }\, \right\rangle =\exp \left( -\frac{\text{i}t}{\hbar }{{H}_{0}} \right)$ where

$$\begin{eqnarray*}&&{{H}_{0}}\equiv \hbar {{\omega }_{p}}{{a}^{\dagger }}a+\hbar {{\omega }_{p}}\underset{n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\sigma _{ee}^{\left( n \right)}+\hbar \left( {{\omega }_{p}}+{{\omega }_{cf}} \right)\sigma _{rr}^{\left( n \right)}\end{eqnarray*}$$

and perform the rotating wave approximation to get the new Hamiltonian $\tilde{H}\,={{\tilde{H}\,}_{a}}+{{\tilde{H}\,}_{c}}+{{\tilde{V}\,}_{a-c}}$, where

$$\begin{eqnarray*} {{\tilde{H}\,}_{a}} & = & -\hbar {{\Delta }_{e}}\underset{n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\sigma _{ee}^{\left( n \right)}-\hbar {{\Delta }_{r}}\underset{n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\sigma _{rr}^{\left( n \right)} \\ {} & {} & +\frac{\hbar {{\Omega }_{cf}}}{2}\underset{N}{\overset{n=1}{\mathop{\mathop{\sum }^{}}}}\,\left( \sigma _{re}^{\left( n \right)}+\sigma _{er}^{\left( n \right)} \right)+\underset{N}{\overset{m<n=1}{\mathop{\mathop{\sum }^{}}}}\,\hbar {{\kappa }_{mn}}\sigma _{rr}^{\left( m \right)}\sigma _{rr}^{\left( n \right)} \\ {{\tilde{H}\,}_{c}} & = & -\hbar {{\Delta }_{c}}{{a}^{\dagger }}a+\hbar \alpha \left( a+{{a}^{\dagger }} \right) \\ {{\tilde{V}\,}_{a-c}} & = & \underset{n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\hbar g\left( a\sigma _{eg}^{\left( n \right)}+{{a}^{\dagger }}\sigma _{ge}^{\left( n \right)} \right) \\ \end{eqnarray*}$$

with the detunings ${{\Delta }_{c}}\equiv \left( {{\omega }_{p}}-{{\omega }_{c}} \right)$, ${{\Delta }_{e}}\equiv \left( {{\omega }_{p}}-{{\omega }_{e}} \right)$, and ${{\Delta }_{r}}\equiv \left( {{\omega }_{p}}+{{\omega }_{cf}}-{{\omega }_{r}} \right)$.

The corresponding Heisenberg–Langevin equations are:

$$\begin{eqnarray}&&\frac{d}{dt}a=\left( \text{i}{{\Delta }_{c}}-{{\gamma }_{c}} \right)a-\text{i}\alpha -\text{i}g\underset{i}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\sigma _{ge}^{\left( i \right)}+{{a}_{in}}\end{eqnarray} \tag{ A.1 }$$

$$\begin{eqnarray}&&\frac{d}{dt}\sigma _{ge}^{\left( i \right)}=\left( \text{i}{{\Delta }_{e}}-{{\gamma }_{e}} \right)\sigma _{ge}^{\left( i \right)}-\text{i}\frac{{{\Omega }_{cf}}}{2}\sigma _{gr}^{\left( i \right)}+\text{i}ga\left( \sigma _{ee}^{\left( i \right)}-\sigma _{gg}^{\left( i \right)} \right)+F_{ge}^{\left( i \right)}\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray} \frac{d}{dt}\sigma _{gr}^{\left( i \right)} & = & \left( \text{i}{{\Delta }_{r}}-{{\gamma }_{r}} \right)\sigma _{gr}^{\left( i \right)}-\text{i}\frac{{{\Omega }_{cf}}}{2}\sigma _{ge}^{\left( i \right)}+\text{i}ga\sigma _{er}^{\left( i \right)} \\ {} & {} & -\text{i}\sigma _{gr}^{\left( i \right)}\underset{j\ne i}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,{{\kappa }_{ij}}\sigma _{rr}^{\left( j \right)}+F_{gr}^{\left( i \right)} \\ \end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray} \frac{d}{dt}\sigma _{er}^{\left( i \right)} & = & \left\{ \text{i}\left( {{\Delta }_{r}}-{{\Delta }_{e}} \right)-{{\gamma }_{er}} \right\}\sigma _{er}^{\left( i \right)}+\text{i}\frac{{{\Omega }_{cf}}}{2}\left( \sigma _{rr}^{\left( i \right)}-\sigma _{ee}^{\left( i \right)} \right) \\ {} & {} & +\text{i}g{{a}^{\dagger }}\sigma _{gr}^{\left( i \right)}-\;\text{i}\sigma _{er}^{\left( i \right)}\underset{j\ne i}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,{{\kappa }_{ij}}\sigma _{rr}^{\left( j \right)}+F_{er}^{\left( i \right)} \\ \end{eqnarray} \tag{ A.4 }$$

where a_in and $F_{\alpha \beta }^{\left( i \right)}$ denote Langevin forces.

A.2. Elimination of the intermediate state

Let us now simplify the system. First, one deduces from equation (A.4) that σ_er is of second order in the small feeding constant α. The term $a\sigma _{er}^{\left( i \right)}$ can, therefore, be neglected in equation (A.3). Moreover, since the ground-state population remains dominant during the evolution of the system, we can write $\sigma _{ee}^{\left( i \right)}-\sigma _{gg}^{\left( i \right)}\simeq -\mathbb{I}$ ; from equation (A.2), the steady-state solution for $\sigma _{ge}^{\left( i \right)}$ in the far detuned region is thus

$$\begin{eqnarray*}&&\sigma _{ge}^{\left( i \right)}\simeq \frac{{{\Omega }_{cf}}}{2\left( {{\Delta }_{e}}+\text{i}{{\gamma }_{e}} \right)}\sigma _{gr}^{\left( i \right)}+\frac{g}{\left( {{\Delta }_{e}}+\text{i}{{\gamma }_{e}} \right)}a+\frac{\text{i}}{\left( {{\Delta }_{e}}+\text{i}{{\gamma }_{e}} \right)}F_{ge}^{\left( i \right)}\end{eqnarray*}$$

Finally, substituting this relation into equations (A.1) and (A.3) one gets

$$\begin{eqnarray}&&\frac{d}{dt}a=\left( \text{i}{{\tilde{\Delta }\,}_{c}}-{{\tilde{\gamma }\,}_{c}} \right)a-\text{i}\alpha +\text{i}{{g}_{eff}}\left( \underset{i}{\mathop{\mathop{\sum }^{}}}\,\sigma _{gr}^{\left( i \right)} \right)+{{\tilde{a}\,}_{in}}\end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray}&&\frac{d}{dt}\sigma _{gr}^{\left( i \right)}=\left( \text{i}{{\tilde{\Delta }\,}_{r}}-{{\tilde{\gamma }\,}_{r}} \right)\sigma _{gr}^{\left( i \right)}+\text{i}{{g}_{eff}}a-\text{i}\sigma _{gr}^{\left( i \right)}\left( \underset{j\ne i}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,{{\kappa }_{ij}}\sigma _{rr}^{\left( j \right)} \right)+\tilde{F}\,_{gr}^{\left( i \right)}\end{eqnarray} \tag{ A.6 }$$

where

$$\begin{eqnarray*}\begin{matrix} {{\tilde{\Delta }\,}_{c}} & = & {{\Delta }_{c}}-{{\Delta }_{e}}\frac{{{g}^{2}}N}{\left( \Delta _{e}^{2}+\gamma _{e}^{2} \right)} \\ {{\tilde{\gamma }\,}_{c}} & = & {{\gamma }_{c}}+{{\gamma }_{e}}\frac{{{g}^{2}}N}{\left( \Delta _{e}^{2}+\gamma _{e}^{2} \right)} \\ {{\tilde{\Delta }\,}_{r}} & = & {{\Delta }_{r}}-{{\Delta }_{e}}\frac{\Omega _{cf}^{2}}{4\left( \Delta _{e}^{2}+\gamma _{e}^{2} \right)} \\ {{\tilde{\gamma }\,}_{r}} & = & {{\gamma }_{r}}+{{\gamma }_{e}}\frac{\Omega _{cf}^{2}}{4\left( \Delta _{e}^{2}+\gamma _{e}^{2} \right)} \\ {{g}_{eff}} & = & \frac{g{{\Omega }_{cf}}}{2\left( {{\Delta }_{e}}+\text{i}{{\gamma }_{e}} \right)}\approx \frac{g{{\Omega }_{cf}}}{2{{\Delta }_{e}}} \\ &&\end{matrix}\end{eqnarray*}$$

are the parameters for the effective two-level model and ${{\tilde{a}\,}_{in}},\tilde{F}\,_{gr}^{\left( i \right)}$ are the modified Langevin noise operators

$$\begin{eqnarray*}\begin{matrix} {{\tilde{a}\,}_{in}} & = & {{a}_{in}}+\frac{g}{\left( {{\Delta }_{e}}+\text{i}{{\gamma }_{e}} \right)}\underset{i}{\mathop{\mathop{\sum }^{}}}\,F_{ge}^{\left( i \right)}\approx {{a}_{in}}+\frac{g}{{{\Delta }_{e}}}\underset{i}{\mathop{\mathop{\sum }^{}}}\,F_{ge}^{\left( i \right)} \\ \tilde{F}\,_{gr}^{\left( i \right)} & = & F_{gr}^{\left( i \right)}+\frac{{{\Omega }_{cf}}}{2\left( {{\Delta }_{e}}+\text{i}{{\gamma }_{e}} \right)}F_{ge}^{\left( i \right)}\approx F_{gr}^{\left( i \right)}+\frac{{{\Omega }_{cf}}}{2{{\Delta }_{e}}}F_{ge}^{\left( i \right)} \\ &&\end{matrix}\end{eqnarray*}$$

Note that, in the absence of collisional terms, one simply recovers the standard three-level EIT susceptibility in the far-detuned region

$$\begin{eqnarray*}&&\frac{da}{dt}=\left( \text{i}{{\Delta }_{c}}-{{\gamma }_{c}}-\frac{{{g}^{2}}N}{\frac{\Omega _{cf}^{2}}{4({{\gamma }_{r}}-\text{i}{{\Delta }_{r}})}-\text{i}{{\Delta }_{e}}} \right)a-\text{i}\alpha +{{\tilde{a}\,}_{in}}\end{eqnarray*}$$

Finally, we get the effective Hamiltonian

$$\begin{eqnarray} \tilde{H}\, & = & -\hbar {{\tilde{\Delta }\,}_{r}}\left( \underset{n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\sigma _{rr}^{\left( n \right)} \right)+\underset{m<n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\hbar {{\kappa }_{mn}}\sigma _{rr}^{\left( m \right)}\sigma _{rr}^{\left( n \right)} \\ {} & {} & -\hbar {{\tilde{\Delta }\,}_{c}}{{a}^{\dagger }}a+\hbar \alpha \left( a+{{a}^{\dagger }} \right)+\hbar {{g}_{\text{eff}}}\left\{ a\left( \underset{n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\sigma _{rg}^{\left( n \right)} \right)+h.c. \right\} \\ \end{eqnarray}$$

A.3. Rybderg bubble approximation

As described in the main text, we introduce the Rydberg bubble approximation. In this approach, the strong Rydberg interactions are assumed to effectively split the sample into ${{\mathcal{N}}_{b}}$ bubbles $\left\{ {{\mathcal{B}}_{\alpha =1,\ldots ,{{\mathcal{N}}_{b}}}} \right\}$, each of which contains ${{n}_{b}}=\left( \frac{N}{{{\mathcal{N}}_{b}}} \right)$ atoms but can only accomodate a single Rydberg excitation, delocalized over the bubble. Note that the number of atoms per bubble n_b is approximately given by [9]

$$\begin{eqnarray*}&&{{n}_{b}}=\frac{2{{\pi }^{2}}{{\rho }_{\text{at}}}}{3}\sqrt{\frac{\left| {{C}_{6}} \right|}{{{\Delta }_{r}}-\Omega _{cf}^{2}/(4{{\Delta }_{e}})}}\end{eqnarray*}$$

where ${{\rho }_{\text{at}}}$ is the atomic density. Each bubble can, therefore, be viewed as an effective spin $\frac{1}{2}$ whose Hilbert space is spanned by

$$\begin{eqnarray} \left| {{-}_{\alpha }} \right\rangle & = & \left| {{G}_{\alpha }} \right\rangle \equiv \underset{{{i}_{\alpha }}\in {{\mathcal{B}}_{\alpha }}}{\mathop{\otimes }}\,\left| {{g}_{{{i}_{\alpha }}}} \right\rangle \\ \left| {{+}_{\alpha }} \right\rangle & = & \left| {{R}_{\alpha }} \right\rangle \equiv \frac{1}{\sqrt{{{n}_{b}}}}\left\{ \left| rg\ldots g \right\rangle +\ldots +\left| g\ldots gr \right\rangle \right\} \\ \end{eqnarray}$$

the ground state of the bubble ${{\mathcal{B}}_{\alpha }}$, and its symmetric singly Rydberg excited state, respectively. Introducing the bubble Pauli operators $\text{s}_{-}^{\left( \alpha \right)}=\hbar \left| {{-}_{\alpha }} \right\rangle \left\langle {{+}_{\alpha }} \right|$, the operator $\text{s}_{-}^{\left( \alpha \right)}$ corresponds to the lowering operator of the spin and the annihilation of a Rydberg excitation, one can write

$$\begin{eqnarray} \underset{n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\sigma _{gr}^{\left( n \right)} & = & \underset{\alpha =1}{\overset{{{\mathcal{N}}_{b}}}{\mathop{\mathop{\sum }^{}}}}\,\underset{{{i}_{\alpha }}\in {{\mathcal{B}}_{\alpha }}}{\mathop{\mathop{\sum }^{}}}\,\sigma _{gr}^{\left( {{i}_{\alpha }} \right)} \\ {} & \approx & \underset{\alpha =1}{\overset{{{\mathcal{N}}_{b}}}{\mathop{\mathop{\sum }^{}}}}\,\frac{\text{s}_{-}^{\left( \alpha \right)}}{\hbar }\left\langle {{-}_{\alpha }} \right|\underset{{{i}_{\alpha }}\in {{\mathcal{B}}_{\alpha }}}{\mathop{\mathop{\sum }^{}}}\,\sigma _{gr}^{\left( {{i}_{\alpha }} \right)}\left| {{+}_{\alpha }} \right\rangle \\ {} & \approx & \sqrt{{{n}_{b}}}\underset{\alpha =1}{\overset{{{\mathcal{N}}_{b}}}{\mathop{\mathop{\sum }^{}}}}\,\frac{\text{s}_{-}^{\left( \alpha \right)}}{\hbar } \\ {} & = & \sqrt{{{n}_{b}}}\frac{{{\text{J}}_{-}}}{\hbar } \\ \end{eqnarray}$$

where we introduced the collective angular momentum ${{\text{J}}_{-}}\equiv \mathop{\sum }^{}_{\alpha =1}^{{{\mathcal{N}}_{b}}}\text{s}_{-}^{\left( \alpha \right)}$. In the same way,

$$\begin{eqnarray} \underset{n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,\sigma _{rr}^{\left( n \right)} & = & \underset{\alpha =1}{\overset{{{\mathcal{N}}_{b}}}{\mathop{\mathop{\sum }^{}}}}\,\underset{{{i}_{\alpha }}\in {{\mathcal{B}}_{\alpha }}}{\mathop{\mathop{\sum }^{}}}\,\sigma _{rr}^{\left( {{i}_{\alpha }} \right)} \\ {} & \approx & \underset{\alpha =1}{\overset{{{\mathcal{N}}_{b}}}{\mathop{\mathop{\sum }^{}}}}\,\left| {{+}_{\alpha }} \right\rangle \left\langle {{+}_{\alpha }} \right|\left\langle {{+}_{\alpha }} \right|\underset{{{i}_{\alpha }}\in {{\mathcal{B}}_{\alpha }}}{\mathop{\mathop{\sum }^{}}}\,\sigma _{rr}^{\left( {{i}_{\alpha }} \right)}\left| {{+}_{\alpha }} \right\rangle \\ {} & \approx & \underset{\alpha =1}{\overset{{{\mathcal{N}}_{b}}}{\mathop{\mathop{\sum }^{}}}}\,\left( \frac{1}{2}+\frac{\text{s}_{z}^{\left( \alpha \right)}}{\hbar } \right) \\ {} & \approx & \left( \frac{{{\mathcal{N}}_{b}}}{2}+\frac{{{\text{J}}_{z}}}{\hbar } \right) \\ \end{eqnarray}$$

where we used $\left| {{+}_{\alpha }} \right\rangle \left\langle {{+}_{\alpha }} \right| \equiv \left( \frac{1}{2}+\frac{\text{s}_{z}^{\left( \alpha \right)}}{\hbar } \right)$. Finally, the Hamiltonian of the system takes the approximate form

$$\begin{eqnarray*} \tilde{H}\, & \approx & -\hbar {{\tilde{\Delta }\,}_{c}}{{a}^{\dagger }}a+\hbar \alpha \left( a+{{a}^{\dagger }} \right) \\ {} & {} & -\hbar {{\tilde{\Delta }\,}_{r}}\left( \frac{{{\mathcal{N}}_{b}}}{2}+\frac{{{\text{J}}_{z}}}{\hbar } \right) \\ {} & {} & +{{g}_{\text{eff}}}\sqrt{{{\mathcal{N}}_{b}}}\left( a{{\text{J}}_{+}}+{{a}^{\dagger }}{{\text{J}}_{-}} \right) \\ \end{eqnarray*}$$

which represents the interaction of the large spin ${{\text{J}}_{-}}$ with the cavity mode a.

A.4. Region of a large number of bubbles and a low number of excitations

From the well-known relation ${{\text{J}}_{+}}{{\text{J}}_{-}}={{\vec{J}}^{2}}-\text{J}_{z}^{2}+\hbar {{\text{J}}_{z}}$, we deduce the second-order operator equation

$$\begin{eqnarray*}&&\text{J}_{z}^{2}-\hbar {{\text{J}}_{z}}-{{\hbar }^{2}}\frac{{{\mathcal{N}}_{b}}}{2}\left( \frac{{{\mathcal{N}}_{b}}+2}{2} \right)+{{\text{J}}_{+}}{{\text{J}}_{-}}=0\end{eqnarray*}$$

In the region of a large number of bubbles ${{\mathcal{N}}_{b}}\gg 1$ and for low excitation numbers, i.e. eigenstates of the total angular momentum $\left| j=\frac{{{\mathcal{N}}_{b}}}{2};m=-\frac{{{\mathcal{N}}_{b}}}{2}+k \right\rangle$ with $k\ll {{\mathcal{N}}_{b}}$, the solution of this equation is approximately given by

$$\begin{eqnarray*}&&{{\text{J}}_{z}}\approx \hbar \left\{ -\frac{{{\mathcal{N}}_{b}}}{2}+\frac{{{\text{J}}_{+}}{{\text{J}}_{-}}}{{{\hbar }^{2}}\left( {{\mathcal{N}}_{b}}+1 \right)}+\frac{{{\left( {{\text{J}}_{+}}{{\text{J}}_{-}} \right)}^{2}}}{{{\hbar }^{4}}{{\left( {{\mathcal{N}}_{b}}+1 \right)}^{3}}} \right\}\end{eqnarray*}$$

whence, at the lowest order in the excitation number,

$$\begin{eqnarray}&&\left( \frac{{{\mathcal{N}}_{b}}}{2}+\frac{{{\text{J}}_{z}}}{\hbar } \right)\approx \frac{{{\text{J}}_{+}}{{\text{J}}_{-}}}{{{\hbar }^{2}}\left( {{\mathcal{N}}_{b}}+1 \right)}+\frac{{{\left( {{\text{J}}_{+}}{{\text{J}}_{-}} \right)}^{2}}}{{{\hbar }^{4}}{{\left( {{\mathcal{N}}_{b}}+1 \right)}^{3}}}\end{eqnarray} \tag{ A.7 }$$

$$\begin{eqnarray}&&\left[ {{\text{J}}_{+}},{{\text{J}}_{-}} \right]\approx -{{\hbar }^{2}}{{\mathcal{N}}_{b}}\end{eqnarray} \tag{ A.8 }$$

Injecting equation (A.7) into the previous form of the Hamiltonian, we get

$$\begin{eqnarray*} \tilde{H}\, & \approx & -\hbar {{\tilde{\Delta }\,}_{c}}{{a}^{\dagger }}a+\hbar \alpha \left( a+{{a}^{\dagger }} \right) \\ {} & {} & -\hbar {{\tilde{\Delta }\,}_{r}}\left( \frac{{{\text{J}}_{+}}{{\text{J}}_{-}}}{{{\hbar }^{2}}\left( {{\mathcal{N}}_{b}}+1 \right)}+\frac{{{\left( {{\text{J}}_{+}}{{\text{J}}_{-}} \right)}^{2}}}{{{\hbar }^{4}}{{\left( {{\mathcal{N}}_{b}}+1 \right)}^{3}}} \right) \\ {} & {} & +\hbar {{g}_{\text{eff}}}\sqrt{N}\left( a\frac{{{\text{J}}_{+}}}{\hbar \sqrt{{{\mathcal{N}}_{b}}}}+{{a}^{\dagger }}\frac{{{\text{J}}_{-}}}{\hbar \sqrt{{{\mathcal{N}}_{b}}}} \right) \\ \end{eqnarray*}$$

Moreover, from equation (A.8), we deduce that the operator $b\equiv \frac{{{J}_{-}}}{\hbar \sqrt{{{\mathcal{N}}_{b}}}}$ is approximately bosonic and, therefore, the Hamiltonian can finally be put under the form

$$\begin{eqnarray}&&\tilde{H}\,\approx -\hbar {{\tilde{\Delta }\,}_{c}}{{a}^{\dagger }}a+\hbar \alpha \left( a+{{a}^{\dagger }} \right)-\hbar {{\tilde{\Delta }\,}_{r}}{{b}^{\dagger }}b-\frac{\hbar \bar{\kappa }}{2}{{b}^{\dagger }}{{b}^{\dagger }}bb+\hbar {{g}_{\text{eff}}}\sqrt{N}\left( a{{b}^{\dagger }}+{{a}^{\dagger }}b \right)\end{eqnarray} \tag{ A.9 }$$

where $\bar{\kappa }\equiv 2{{\tilde{\Delta }\,}_{r}}/{{\mathcal{N}}_{b}}$.

By definition, the second-order correlation function for the outgoing field is

$$\begin{eqnarray*}&&g_{out}^{\left( 2 \right)}({{t}_{1}},{{t}_{2}})=\frac{\left\langle a_{out}^{\dagger }\left( {{t}_{1}} \right)a_{out}^{\dagger }\left( {{t}_{2}} \right){{a}_{out}}\left( {{t}_{2}} \right){{a}_{out}}\left( {{t}_{1}} \right) \right\rangle }{\left\langle a_{out}^{\dagger }\left( {{t}_{2}} \right){{a}_{out}}\left( {{t}_{2}} \right) \right\rangle \left\langle a_{out}^{\dagger }\left( {{t}_{1}} \right){{a}_{out}}\left( {{t}_{1}} \right) \right\rangle }\end{eqnarray*}$$

Using the relations [14]

$$\begin{eqnarray*}\begin{matrix} \left\langle a_{out}^{\dagger }\left( t \right){{a}_{out}}\left( t \right) \right\rangle & = & {{\gamma }_{c}}\left\langle {{a}^{\dagger }}\left( t \right)a\left( t \right)\right\rangle \\ {{a}_{out}}\left( t \right) & = & \sqrt{2{{\gamma }_{c}}}a\left( t \right)-{{a}_{in}}\left( t \right) \\ &&\end{matrix}\end{eqnarray*}$$

and keeping only nonzero terms (all terms like $\left\langle a_{in}^{\dagger }....\right\rangle$ and $\left\langle ....{{a}_{in}}\right\rangle$ equal zero), one obtains in the numerator four nonzero terms

$$\begin{matrix} \left\langle {{a}^{\dagger }}({{t}_{1}}){{a}^{\dagger }}({{t}_{2}}){{a}_{in}}({{t}_{2}})a({{t}_{1}})\right\rangle \\ \left\langle {{a}^{\dagger }}({{t}_{1}}){{a}^{\dagger }}({{t}_{2}})a({{t}_{2}})a({{t}_{1}})\right\rangle \\ \left\langle {{a}^{\dagger }}({{t}_{1}})a_{in}^{\dagger }({{t}_{2}})a({{t}_{2}})a({{t}_{1}})\right\rangle \\ \left\langle {{a}^{\dagger }}({{t}_{1}})a_{in}^{\dagger }({{t}_{2}}){{a}_{in}}({{t}_{2}})a({{t}_{1}})\right\rangle \\ \end{matrix}$$

Let us consider the first term. Using the standard commutation relations between a and a_in operators, we have:

$$\begin{eqnarray*} \left\langle {{a}^{\dagger }}({{t}_{1}}){{a}^{\dagger }}({{t}_{2}}){{a}_{in}}({{t}_{2}})a({{t}_{1}})\right\rangle & = & \left\langle {{a}^{\dagger }}({{t}_{1}}){{a}^{\dagger }}({{t}_{2}})a({{t}_{1}}){{a}_{in}}({{t}_{2}})\right\rangle \\ {} & {} & \left\langle +{{a}^{\dagger }}({{t}_{1}}){{a}^{\dagger }}({{t}_{2}})\left[ {{a}_{in}}({{t}_{2}}),a({{t}_{1}}) \right]\right\rangle \\ {} & = & \sqrt{2{{\gamma }_{c}}}\theta ({{t}_{1}}-{{t}_{2}})\left\langle {{a}^{\dagger }}({{t}_{1}}){{a}^{\dagger }}({{t}_{2}})\left[ a({{t}_{2}}),a({{t}_{1}}) \right]\right\rangle \\ \end{eqnarray*}$$

Here we used the relation

$$\begin{eqnarray*}&&\left[ X\left( {{t}_{1}} \right),{{a}_{in}}\left( {{t}_{2}} \right) \right]=\sqrt{2{{\gamma }_{c}}}\theta \left( {{t}_{1}}-{{t}_{2}} \right)\left[ X,a \right]\end{eqnarray*}$$

where X is any system operator [14] and where $\theta \left( \tau \right)$ is the Heaviside step-function (with $\theta \left( 0 \right)=\frac{1}{2}$). Evaluating the other terms in the same way, one finally obtains

$$\begin{eqnarray*}&&g_{out}^{(2)}\left( {{t}_{1}},{{t}_{2}} \right)=\frac{\left\langle {{a}^{\dagger }}\left( {{t}_{m}} \right){{a}^{\dagger }}\left( {{t}_{M}} \right)a\left( {{t}_{M}} \right)a\left( {{t}_{m}} \right) \right\rangle }{\left\langle {{a}^{\dagger }}\left( {{t}_{1}} \right)a\left( {{t}_{1}} \right) \right\rangle \left\langle {{a}^{\dagger }}\left( {{t}_{2}} \right)a\left( {{t}_{2}} \right) \right\rangle }\end{eqnarray*}$$

where ${{t}_{m}}\equiv \min \left( {{t}_{1}},{{t}_{2}} \right)$ and ${{t}_{M}}\equiv \max \left( {{t}_{1}},{{t}_{2}} \right)$.