Creation and detection of photonic molecules in Rydberg gases

TL;DR

This paper explores how photons in Rydberg gases can form bound states called photonic molecules under electromagnetically induced transparency, analyzing their creation, dynamics, and detection through theoretical and numerical methods.

Contribution

It introduces a detailed Green's function and numerical wave-function approach to analyze photonic molecule formation and their interplay with scattering states in Rydberg gases.

Findings

Bound states of photons can be formed and detected.

Scattering states acquire a robust phase useful for separation.

Analytic solutions clarify conditions for molecule formation.

Abstract

We consider the propagation of photons in a gas of Rydberg atoms under conditions of electromagnetically induced transparency, where they form strongly interacting massive particles, termed Rydberg polaritons. Depending on the strength of the van der Waals-type interactions of the atoms either bunching or anti-bunching of photons can be observed when driving the atoms off-resonantly. The bunching is associated with the formation of bound states. We employ a Green's function approach and numerical wave-function simulations to analyze the conditions for the creation and the dynamics of these photonic molecules and their interplay with the scattering continuum which can also show photon bunching. Analytic solutions of the pair-propagation problem obtained from a pseudopotential approximation and verified numerically provide a detailed understanding of bound and scattering states. We find that the scattering contributions acquire asymptotically a robust relative phase which can be employed to separate bound-state and scattering contributions by a homodyne detection scheme.

Equations

$$\mathrm{i}\frac{\partial }{\partial t}\hat{\mathcal{E}}(\mathbf{r})$$

$$\mathrm{i}\frac{\partial }{\partial t}\hat{\mathcal{P}}(\mathbf{r})$$

$$\mathrm{i}\frac{\partial }{\partial t}\hat{\mathcal{S}}(\mathbf{r})$$

$$\mathrm{i}\frac{\partial }{\partial t}\hat{\mathcal{E}}(z)$$

$$\mathrm{i}\frac{\partial }{\partial t}\hat{\mathcal{S}}(z)$$

$$+\int \mathrm{d}{z}^{\prime }V(z-z^{\prime }){\hat{\mathcal{S}}}^{†}(z^{\prime })\hat{\mathcal{S}}(z^{\prime })\hat{\mathcal{S}}(z),$$

$$\mathrm{i}\frac{\partial }{\partial t}\left(\begin{array}{c}\hat{\mathcal{E}}\\ \hat{\mathcal{S}}\end{array}\right)={\hat{H}}_0\left(\begin{array}{c}\hat{\mathcal{E}}\\ \hat{\mathcal{S}}\end{array}\right),{\hat{H}}_0=-\mathrm{i}\left(\begin{array}{cc}c\frac{\partial }{\partial z}+\frac{g^2}{\Gamma }& \frac{g\Omega }{\Gamma }\\ \frac{g\Omega }{\Gamma }& \frac{{\Omega }^2}{\Gamma }\end{array}\right).$$

$$\mathrm{i}\frac{\partial }{\partial t}\mathbf{\Psi }=\left\{{\hat{\mathcal{H}}}_0(z_1,z_2)+V(z_1-z_2){\hat{\mathrm{P}}}_{\mathcal{SS}}\right\}\mathbf{\Psi }$$

$$\mathbf{\Psi }(K,r,0)=f(K,r)|{\phi }_1⟩,$$

$$\mathcal{EE}(K,r,t)=\frac{{\cos }^4\theta }{2\pi \mathrm{i}}\iint \mathrm{d}\omega \mathrm{d}{r}^{\prime }{\mathrm{e}}^{-\mathrm{i}\omega t}G(r,r^{\prime },\omega )f(K,r^{\prime }).$$

$$\begin{gathered} G(r,r^{\prime },\omega )=G_0(r,r^{\prime },\omega ) \\ -{\sin }^4\theta \int \mathrm{d}{r}^{\prime \prime }{G}_0(r,r^{\prime \prime },\omega )W(r^{\prime \prime })G(r^{\prime \prime },r^{\prime },\omega ), \end{gathered}$$

$$G_0(r,r^{\prime },\omega )=-\frac{\exp \{\mathrm{i}\sqrt{2m(\omega -v_{\mathrm{g}}K)}|r-r^{\prime }|\}}{2\mathrm{i}\sqrt{2m(\omega -v_{\mathrm{g}}K)}}.$$

$${\hat{\mathcal{H}}}_{\mathrm{eff}}=-\frac{1}{2m}\frac{{\mathrm{d}}^2}{\mathrm{d}{r}^2}+W(r){\sin }^4\theta .$$

$$m=\mathrm{i}\frac{g^2}{4c\Gamma v_{\mathrm{g}}}=\frac{\mathrm{sign}(\Delta )}{4v_{\mathrm{g}}{L}_{\text{abs}}}\left(1+i\frac{\gamma }{\Delta }\right)$$

$$W(r)\equiv \frac{V(r)}{1+\alpha V(r)},\alpha =\frac{i\gamma -\Delta }{2{\Omega }^2}.$$

$$\chi (r)={\chi }^{\prime }+\mathrm{i}{\chi }^{\prime \prime }=-\mathrm{i}\frac{g^2}{{\Omega }^2}W(r).$$

$$G(r,r^{\prime },\omega )=\sum _{n=1}^N\frac{{\psi }_n(r){\psi }_n^{\ast }(r^{\prime })}{\omega -E_n}+\int \mathrm{d}E\frac{{\psi }_E(r){\psi }_E^{\ast }(r^{\prime })}{\omega -E},$$

$${\int }_{-\infty }^{\infty }\mathrm{d}rmW(r)<0.$$

$$N\le 1+2\left|m\right|{\int }_{-\infty }^{\infty }\mathrm{d}r\left|r\right|W(r).$$

$$\xi \le \sqrt{3\sqrt{3}/\pi }\approx 1.2861$$

$$\left|E_n\right|\le \frac{|m|}{4}{\left({\int }_{-\infty }^{\infty }\mathrm{d}r\left|W(r)\right|\right)}^2.$$

$$|E_0|\le \frac{1}{2}{\xi }^2{\left(\frac{2\pi }{3}\right)}^2\frac{2{\Omega }^2}{\Delta }\lesssim \frac{2{\Omega }^2}{\Delta },$$

$$\Delta p={\int }_{-\infty }^{\infty }\mathrm{d}r{\left(\frac{d{\psi }_0(r)}{dr}\right)}^2\lesssim {\int }_{-\infty }^{\infty }\mathrm{d}r|2mW(r)|.$$

$$r_{\mathrm{b}}={\int }_{-\infty }^{\infty }\mathrm{d}r{r}^2{\psi }_0^2(r)\gtrsim \frac{1}{2}{\left({\int }_{-\infty }^{\infty }\mathrm{d}r|2mW(r)|\right)}^{-1},$$

$$r_{\mathrm{b}}\gtrsim \frac{3}{4\pi }\frac{L_{\text{abs}}}{\xi }>L_{\text{abs}}>R_{\mathrm{B}},$$

$$\frac{\mathcal{EE}(r,t)}{{\cos }^4\theta }=C_0{\mathrm{e}}^{-\mathrm{i}{E}_{0}t}{\psi }_0(r)+\int C(E){\mathrm{e}}^{-\mathrm{i}Et}{\psi }_E(r)\mathrm{d}E,$$

$${\mathcal{E}\mathcal{S}}_{+}(r,t)=\mathcal{E}\mathcal{S}+\mathcal{S}\mathcal{E}\approx -2C_0{\cos }^3\theta {\psi }_0(r){\mathrm{e}}^{-\mathrm{i}{E}_{0}t}.$$

$$\mathcal{SS}(r,t)\approx C_0\frac{{\cos }^2\theta }{1-\frac{\Delta }{2{\Omega }^2}V(r)}{\psi }_0(r){\mathrm{e}}^{-\mathrm{i}{E}_{0}t}.$$

$$\frac{\left|{\mathcal{E}\mathcal{S}}_{-}(r,t)\right|}{2C_0{\cos }^3\theta {\psi }_0(r)}\approx \left|\frac{L_{\text{abs}}}{{\psi }_0(r)}\frac{d{\psi }_0(r)}{dr}\right|\to 0$$

$$\mathbf{\Psi }(r,t)={\cos }^2\theta C_0\left(\begin{array}{c}{\cos }^2\theta \\ -\cos \theta \\ -\cos \theta \\ \frac{1}{1-\frac{\Delta }{2{\Omega }^2}V(r)}\end{array}\right){\psi }_0(r){\mathrm{e}}^{-\mathrm{i}{E}_{0}t},$$

$$\begin{array}{c}\mathrm{i}\frac{\partial }{\partial t}G(r,r^{\prime },t)=\left(\frac{1}{2m}\frac{{\partial }^2}{\partial r^2}+W(r)\right)G(r,r^{\prime },t),\\ G(r,r^{\prime },0)=\delta (r-r^{\prime }).\end{array}$$

$$W(r)\to W_{\text{eff}}(r)=\frac{2\pi R_{\mathrm{B}}}{3}\frac{2{\Omega }^2}{|\Delta |}\frac{1}{(1+\mathrm{i}\gamma /|\Delta |{)}^{5/6}}\delta (r).$$

Full text

I Introduction

Rydberg gases are of great interest in quantum optics as they enable to mediate strong and long-range nonlinearities between photons. The van der Waals-type interactions between Rydberg states Saffman et al. (2010) in a gas of three-level atoms can be used to create strong photon-photon interactions for light fields coupling to the atomic medium in an scheme of electromagnetically induced transparency (EIT) Fleischhauer et al. (2005); Pritchard et al. (2010). Moreover, due to the long-range nature of the Rydberg-Rydberg interactions, also the photon-photon interactions are long-range which makes this setup a promising candidate for creating and analyzing interesting many-body states with applications ranging from quantum computation Tiarks et al. (2014); Gorniaczyk et al. (2014); Maxwell et al. (2013); Distante et al. (2016); Hofmann et al. (2013) to quantum simulations Hafezi et al. (2013); Ningyuan et al. (2015); Zhang et al. (2012); Jachymski et al. (2016); Maghrebi et al. (2015).

In an EIT setup photons travel as massive quasi-particles, so-called dark-state polaritons (DSPs) with group velocities much smaller than the vacuum speed of light Fleischhauer et al. (2005). If the EIT coupling involves atomic Rydberg states (see Fig.1(a)), interactions between Rydberg atoms are transferred to polaritons. In particular, it has been proposed and observed experimentally that repulsive Rydberg interactions lead to an avoided volume of photons for small distances, i.e., anti-bunching Gorshkov et al. (2011); Peyronel et al. (2012). On the other hand also bunching of photons has been observed, when driving the atoms off-resonantly Firstenberg et al. (2013). This can result from photonic bound states ("photonic molecules"), but also from bunched continuum components (scattering states). The formation of bound states requires an interplay between interactions and dissipation. Numerical simulations show that bunching can only be observed in the regime of weak to moderate polariton interactions, quantified by a small optical depth per blockade, $\xi\ll 1$, which will be defined later-on, see Fig. 1(c)–LABEL:sub@fig:Figure1:d. In the present paper we investigate the properties of photonic bound states, analyze conditions for their formation, and discuss possible ways to distinguish bound- and continuum-state contributions.

Specifically, we apply a Green’s function approach to model one-dimensional systems of interacting Rydberg polaritons and investigate the creation of states leading to the bunching and the dynamics at large times. For short times, bound- and scattering-state contributions are equally important to explain photon bunching. The bound state contribution decays exponentially due to losses, and thus bunching at large times comes solely from scattering states.

We find that the scattering states have a phase which depends only on the ratio of the probe field detuning and the decay rate of the excited atomic state. Therefore, this phase is very robust and can be used to separate bunched photons resulting from bound states and scattering states by homodyne detection, see Fig.1(b).

II Model and two-excitation equations

II.1 Model

We consider a quantized probe field $\hat{\mathcal{E}}$ propagating under conditions of EIT in a medium consisting of $N$ atoms with three levels driven by two optical fields in a ladder scheme as illustrated in Fig. 1(a).

The atomic ground state $\ket{\mathrm{g}}$ and the excited state $\ket{\mathrm{e}}$ are coupled by the quantized probe field $\hat{\mathcal{E}}(\mathbf{r},t)$ with carrier frequency $\omega_{\mathrm{p}}$ and wave vector $\mathbf{k}_{\mathrm{p}}$. The probe field is detuned from the atomic transition by the single-photon detuning $\Delta=\omega_{\mathrm{eg}}-\omega_{\mathrm{p}}$. Furthermore, a classical control field with frequency $\omega_{\mathrm{c}}$ drives the transition $\ket{\mathrm{g}}\leftrightarrow\ket{\mathrm{e}}$ with Rabi frequency $\Omega$ and detuning $\Delta_{\mathrm{c}}$, which is chosen such that the resulting two-photon detuning vanishes, i.e., $\delta=\Delta+\Delta_{\mathrm{c}}=0$. The intermediate state is subject to spontaneous decay with rate $\gamma$. The atoms are described by spin flip operators $\hat{\sigma}_{\mu\nu}^{i}=\ket{\mu}_{ii}\bra{\nu}$ and interact via the van der Waals potential $V(\mathbf{r})=C_{6}/\left|\mathbf{r}\right|^{6}$ in the level $\ket{\mathrm{r}}$. Assuming a homogeneous distribution of atoms we can describe them by coarse-grained continuous operators $\hat{\sigma}_{\mu\nu}(\mathbf{r})$, which for negligible atomic saturation are bosonic fields $\hat{\mathcal{S}}(\mathbf{r})=\hat{\sigma}_{\mathrm{gr}}(\mathbf{r})$ and $\hat{\mathcal{P}}(\mathbf{r})=\hat{\sigma}_{\mathrm{ge}}(\mathbf{r})$.

Finally, from the atom-field coupling Hamiltonian in rotating wave approximation and Maxwell’s equations we obtain in linear response in $g\hat{\mathcal{E}}$ the paraxial Maxwell-Bloch equations,

[TABLE]

where we defined the complex detuning $\Gamma=\gamma+\mathrm{i}\Delta$ and the coupling strength $g=\wp\sqrt{n\omega_{\mathrm{p}}/2\hbar\epsilon_{0}}$ with $\wp$ being the dipole moment of the $\ket{\mathrm{g}}\leftrightarrow\ket{\mathrm{e}}$ transition, and $n$ being the atomic number density. $\hat{F}_{\mathrm{ge}}$ is a Langevin noise operator Louisell (1973), which we introduced to preserve the commutation relations. Under EIT driving conditions the occupation of the level $\ket{\mathrm{e}}$ stays small, and thus the Langevin noise can be neglected.

As shown in Moos et al. (2015), in experimentally relevant situations the interaction can be described by a one-dimensional model and we can neglect the transverse kinetic energy $\tfrac{c}{k_{\mathrm{p}}}\nabla_{\perp}^{2}\hat{\mathcal{E}}$. Finally, assuming the time evolution being slow on the time scale set by the complex detuning $\left|\Gamma\right|$, we adiabatically eliminate the optical polarization $\hat{\mathcal{P}}$, leading to

[TABLE]

which is a set of coupled nonlinear integro-differential equations for the operators $\hat{\mathcal{E}}$ and $\hat{\mathcal{S}}$.

II.2 Dark-state polaritons

Let us first briefly summarize the description of the noninteracting limit, i.e., $V(z)\equiv 0$, which also applies to the case of a single photon propagating through the Rydberg medium. In this case, Eqs. (2) form a set of linear equations that can be expressed as

[TABLE]

The eigenmodes of the Hamiltonian $\hat{H}_{0}$ in the long-wavelength limit ($k\approx 0$) correspond to quasi-particles composed of light and matter excitation, the so-called dark- and bright-state polaritons which can be written as $\hat{\psi}_{\mathrm{d}}=-\cos\theta\hat{\mathcal{E}}+\sin\theta\hat{\mathcal{S}}$ and $\hat{\psi}_{\mathrm{b}}=\sin\theta\hat{\mathcal{E}}+\cos\theta\hat{\mathcal{S}}$, respectively, see, e.g., Fleischhauer et al. (2005). Here the mixing angle $\theta$ is defined by $\tan\theta=g/\Omega$. Treating the momentum $k$ perturbatively one finds that the dark-state polariton propagates lossless with the group velocity $v_{\mathrm{g}}=c\cos^{2}\theta$. Furthermore, it forms a quasi-particle with an effective mass $m$. The mass is approximately $m\approx(2v_{\mathrm{g}}L_{\text{abs}})^{-1}$ under slow-light conditions and in an off-resonant driving scheme, where $L_{\text{abs}}=\left|\Delta\right|c/g^{2}$ is the off-resonant optical depth, In contrast, the bright-state polariton propagates with velocity $c\sin^{2}\theta\approx c$ and is subject to losses with the rate $\gamma\Omega_{\mathrm{e}}^{2}/\left|\Delta\right|^{2}$, where the effective Rabi frequency is defined by $\Omega_{\mathrm{e}}^{2}=g^{2}+\Omega^{2}$.

For large separations between excitations the Rydberg-Rydberg interaction can be included as a perturbation, see, e.g., Otterbach et al. (2013). However, this approach does not capture bound states and is thus not applicable in general. Instead the full scattering problem has to be considered as was done in Bienias et al. (2014); see also Gullans et al. (2016).

II.3 Effective model for two excitations

To analyze the dynamics of interacting excitations we now consider the time evolution of two particles, which can be done by using wave functions $\mathcal{EE}(z_{1},z_{2},t)=\bra{0}\hat{\mathcal{E}}(z_{1},t)\hat{\mathcal{E}}(z_{2},t)\ket{\phi}$, and analogously defined components $\mathcal{ES},\mathcal{SE}$ and $\mathcal{SS}$ that can be combined into the four-component vector $\bm{\Psi}_{2}=(\mathcal{EE},\mathcal{ES},\mathcal{SE},\mathcal{SS})^{T}$. The time evolution of $\bm{\Psi}_{2}$ in real space is governed by the equation

[TABLE]

with $\hat{\mathcal{H}}_{0}=\hat{H}_{0}(z_{1})\otimes\openone_{2}+\openone_{2}\otimes\hat{H}_{0}(z_{2})$. The operator $\hat{\mathrm{P}}_{\mathcal{SS}}=\ket{\varphi_{4}}\bra{\varphi_{4}}$ denotes the projector onto the $\mathcal{SS}$-component of the wave function, i.e., two Rydberg excitations, with $\ket{\varphi_{4}}=(0,0,0,1)^{T}$. This equation can be integrated numerically to find the time evolution of a two-photon wave packet. In particular we simulate the time evolution starting in free space and propagating according to equation (2) through a sharp boundary. We find qualitatively a very different behavior inside the medium depending on the strength of the interaction potential $V(z_{1}-z_{2})$ as can be seen in Figs. 1(c) and 1(d). In the weakly interacting regime, Fig. 1(c), we find a bunching of photons, while in the strongly interacting regime, Fig. 1(d), the photons avoid a volume given by $\left|z_{1}-z_{2}\right|<R_{\mathrm{B}}$, with $R_{\mathrm{B}}=\left(\left|\Delta\right|C_{6}/2\Omega^{2}\right)^{1/6}$ being the off-resonant blockade radius .

To gain analytical insight into these observations, we employ a Green’s function approach to solve the time evolution of the two-photon wave-function, Eq. (4), similar to Bienias et al. (2014). We transform to center-of-mass- and relative coordinates of the two excitations, $R=\tfrac{1}{2}(z_{1}+z_{2})$ and $r=z_{1}-z_{2}$, respectively. Subsequently, we perform a Fourier transform with respect to the center of mass $R$ according to $f(R)=\int\mathrm{d}K\,\mathrm{e}^{iKR}{\tilde{f}}(K)$. Specifically, we consider the initial state

[TABLE]

where $\ket{\varphi_{1}}=(1,0,0,0)^{T}$, i.e., we assume that only the photonic component is present at the beginning of the evolution. Our calculation can easily be generalized to other initial states. Furthermore, we restrict ourselves to the case of negative single-photon detuning, $\Delta<0$. The solution for positive detuning can be derived straightforwardly.

We are interested in the (asymptotic) behavior of the amplitude $\mathcal{EE}(K,r,t)$ at large times. In this limit, the low-frequency contributions are the dominant ones (see Appendix A for more details) and by simple algebraic calculations one obtains for the two-photon amplitude

[TABLE]

Here the Green’s function $G(r,r^{\prime},\omega)$ is the solution of the integral equation

[TABLE]

where $W(r)$ denotes an effective potential which is defined by Eq. (11).

Under the condition that $\Im\{\sqrt{2m(\omega-v_{\mathrm{g}}K)}\}>0$, the free Green’s function $G_{0}$ in Eq. (7) has the coordinate representation

[TABLE]

In the low-energy regime, $|\omega|\ll\Omega^{2}/|\Gamma|$, i.e., for frequencies well inside the EIT window the Green’s function $G$ describes the evolution of a particle with an effective Hamiltonian

[TABLE]

The complex mass is given by

[TABLE]

and the effective interaction potential reads

[TABLE]

In the limit of slow light $g\gg\Omega$ and large single photon detuning $|\Delta|\gg\gamma$ the effective mass reduces to the simpler expression $m\approx\operatorname{sign}(\Delta)(4v_{\mathrm{g}}L_{\text{abs}})^{-1}$, which coincides with the results derived in Zimmer et al. (2008). Likewise the coefficient $\alpha$ simplifies to $\alpha\approx-\Delta/2\Omega^{2}$.

In the following we assume slow-light conditions and set $\sin^{2}\theta\approx 1$.

In Fig. 2 we show the effective potential for positive and negative detuning. For distances larger $R_{\mathrm{B}}$ the potential decays like the bare van der Waals potential, for small distances the potential becomes flat. We can interpret the effective potential as complex susceptibility of a single photon in the presence of a fixed Rydberg excitation at the origin resulting in a space dependent two-photon detuning Otterbach (2011); Gorshkov et al. (2011).

[TABLE]

Note that the complex mass (10) always has a positive imaginary part, effectively describing the (small) polariton losses due to spontaneous decay of the intermediate level $\ket{\mathrm{e}}$, while the sign of its real part can be tuned depending on the sign of the single-photon detuning $\Delta$. The product of the real part of the effective potential and the effective mass is always negative at small distances, suggesting the existence of bound states, independent of the sign of $\Delta$.

III Weakly bound states - photonic molecules

For the interacting problem we have to solve Eq. (7) for the Green’s function. In the far-detuned limit, when $\left|\Delta\right|\gg\gamma$, the Green’s function $G(r,r^{\prime},\omega)$ can be written as a sum,

[TABLE]

of bound eigenstates, denoted by $\psi_{n}(r)$, and continuum eigenstates, denoted by $\psi_{E}(r)$, The binding energies of the molecular states are in general complex and increase with the optical depth per blockade distance, $\xi\equiv R_{\mathrm{B}}/L_{\text{abs}}$.

A sufficient condition for the existence of bound states $\psi_{n}$ in the spectrum of this Hamiltonian is Simon (2005)

[TABLE]

We note that in our case of negative single-photon detuning, the product of $m$ and $W(r)$ is negative, and thus this condition is met.

III.1 Properties of bound states

The bound eigenstates $\psi_{n}(r)$ can be computed by numerical diagonalization of the effective Hamiltonian, Eq. (9), for discretized spatial coordinates on a finite spatial interval. This allows us to get an approximate spectrum of the bound eigenstates as a function of the optical depth per blockade distance, $\xi$, which is shown in Fig. 3(a). For sufficiently small $\xi$ only a single bound state exists and with increasing $\xi$ the number of bound eigenstates grows, as does their energy. As the effective Hamiltonian is only applicable in the regime of small energies, we show only energies with an absolute value smaller than $\frac{\Omega^{2}}{2\left|\Gamma\right|}$, corresponding to frequencies inside the EIT transparency window.

The number of bound states $N$ can be estimated Simon (2005) by

[TABLE]

This leads to a condition for the existence of only one bound state:

[TABLE]

Hence, a unique bound state exists only for small optical depth per blockade $\xi$. Consequently, to observe the formation of sufficiently long-lived photonic molecules in an experiment one hast to operate in this regime. Deeply bound states have small spatial extent, i.e., they are strongly localized and hard to excite by a flat initial photon distribution. The higher-$n$ bound eigenstates $\psi_{n}$, which exist for $\xi\gg 1$, are also hard to excite, since they exhibit many oscillations and, furthermore, are subject to strong decay, as we will show later-on. This explains the behavior seen in Fig. 1(c),LABEL:sub@fig:Figure1:d. The excitation of a bound photon state is only effective if a single bound state close to the continuum exists, i.e., in the weakly interacting limit.

It is well known Abramov et al. (1999) that the bound state energies of a one-dimensional Schrödinger equation with a complex potential $W(r)$ are bounded by

[TABLE]

Making use of this inequality we then obtain the following estimate for the energy of the bound state

[TABLE]

where in the last step we assumed that only a single bound state exists. Hence the bound state energy is inside the low-frequency region of the EIT transparency window.

An estimate for the size of the bound state $\psi_{0}(r)$ can be obtained from the uncertainty of the relative momentum. Assuming that we are in the regime of a single bound state close to the continuum, i.e., $\xi\lesssim 1$, a simple calculation shows that the momentum width of the bound state $\psi_{0}(r)$ is given by

[TABLE]

Using Heisenberg’s uncertainty relation, we can derive an approximate expression for the size of the bound state

[TABLE]

which yields

[TABLE]

where we used $\xi\lesssim 1$. Thus in the parameter regime, where bound states can be excited, their spatial extend is rather large and exceeds the absorption length as well as the blockade radius. In Fig. 3(b) we show an eigenstate for $\xi=0.2$ calculated by numerical diagonalization of the two-photon Hamiltonian compared to an exponential function with the size $\frac{\pi}{3}R_{\mathrm{B}}/\xi^{2}$, showing a very good agreement.

III.2 Bound-state components

The internal structure of the bound state can be found numerically by calculating the time evolution of an initially broadly distributed wave function consisting of two dark-state polaritons for the case of vanishing center-of-mass momentum $K=0$. We assume that the $\mathcal{SS}$-component has initially no excitation inside the blockaded region $\left|r\right|\lesssim R_{\mathrm{B}}$. The results are shown in Fig 4. As expected, the $\mathcal{SS}$-component is strongly suppressed inside the blockade radius. Here the bound state has mainly photonic character. In the case of $\Delta>0$ one recognizes a sharp peak of the $\mathcal{SS}$-component close to the blockade radius. This coincides with the sharp minimum seen in the effective potential for positive detuning at this distance, cf. Fig. 2(a).

If condition (16) for a single bound state is fulfilled, the amplitude $\mathcal{EE}(r,t)$ of finding two photons at relative distance $r$ reads

[TABLE]

where $C_{0}$ and $C(E)$ are the overlap integrals between the initial state and the bound and continuum eigenstates, respectively. We here consider only the first part, corresponding to the bound state and will discuss the continuum states in the following section.The amplitudes of the remaining components $\mathcal{ES},\mathcal{SE}$ and $\mathcal{SS}$, of the bound state can be obtained by substituting the solution for $\mathcal{EE}$ into the two-particle Schrödinger equation. A direct calculation gives

[TABLE]

In obtaining this expression we have assumed that $g^{2}/|\Delta|\gg cK,|E_{0}|$. The calculation of the spin component $\mathcal{SS}(r,t)$ is more involved, but straightforward. After simple algebra we arrive at

[TABLE]

Analogous calculations can be performed for the antisymmetric component $\mathcal{ES}_{-}=\mathcal{ES}-\mathcal{SE}$, which becomes negligible if the size of the bound state is much larger than the off-resonant optical length, $r_{\mathrm{b}}\gg L_{\text{abs}}$.

[TABLE]

In this weak-interaction limit the amplitudes, Eqs. (21) –(23) and (24), can be combined in a compact form, $\mathbf{\Psi}=(\mathcal{EE},\mathcal{ES},\mathcal{SE},\mathcal{SS})^{T}$, yielding

[TABLE]

where the factor $\cos^{2}\theta$ in front appears as a result of projecting the initial state onto the state of two free polaritons and can be changed by choosing a specific initial state vector $\mathbf{\Psi}(r,0)$. Note that this result is only applicable, when the energy $\left|E_{0}\right|$ of the bound state is much smaller than other energies involved in the system, e.g. $2\Omega^{2}/\left|\Delta\right|$. The spatial size of the bound state $\psi_{0}(r)$ is in this case larger than $L_{\text{abs}}$. We observe that $\mathbf{\Psi}(r,t)$ in equation (25) describes a two-photon wave packet that, although subject to decay, propagates form-stable and exhibits bunching for small distances, i.e., a photonic molecule state.

We find relatively good agreement of the different asymptotic forms of the wave function amplitudes, Eq. (25), with the numerical results. There is only a small deviation, as we still find a finite remaining antisymmetric component $\mathcal{ES}_{-}(r,t)$ and correspondingly a slightly increased spin excitation.

IV Bound states and continuum

A major goal of this paper is to analyze the interplay between bound state and continuum contributions. In this section we employ an effective Schrödinger equation to derive asymptotic analytical solutions of Eq. (7) including both, bound and continuum parts, that we compare to numerical solutions. We show that an initial state evolves into a superposition of bound states and scattering states both of which contribute to the bunching signal.

IV.1 Approximate analytic solutions

For the relevant propagation distances the dispersive nature of the interaction potential, i.e., the frequency-dependence of $W(r,\omega)$, can be ignored for the dynamics of the two-particle wave function (see appendix A). In this case we can proceed in a standard way and reformulate the integral equation, Eq. (7), as a Schrödinger-type initial value problem for the propagator $G(r,r^{\prime},t)$,

[TABLE]

In obtaining this result we have omitted the kinematic term $cK\cos^{2}\theta$ for the center-of-mass motion, as it generates a trivial shift in time. In the following we thus assume $K=0$. The potential $W(r)$, (11), can be treated as an effective interaction between two photons. At short times the model (26) does not approximate the full dynamics well, since our derivations of Eq. (26) are based on the assumption that the evolution time should be long compared to all other characteristic time scales of the system (see: appendix A).

As can be seen from the full numerical solution in Fig. 5, for weakly interacting photons, i.e., $\xi\ll 1$, the range of spatial variation of the two-photon amplitude $\mathcal{EE}(r,t)$ is much greater than the range of the potential, i.e., the blockade radius $R_{\mathrm{B}}$ . This suggests that $W(r)$ can be approximated by a delta-like pseudopotential.

[TABLE]

Assuming that the initial two-photon amplitude is uniformly distributed in the relative coordinate $r$ one can show that the Schrödinger equation for $\mathcal{EE}(r,t)$ with the effective interaction potential $W_{\text{eff}}(r)$ admits analytical solutions in closed form at large times. For convenience we introduce dimensionless time and space coordinates that are measured in units of $(2\Omega^{2}/\left|\Delta\right|)^{-1}$ and $R_{\mathrm{B}}$, respectively. After a lengthy but straightforward calculation we find the following expression for the two-photon amplitude $\mathcal{EE}(r,t)$,

[TABLE]

where for convenience we defined the constants

[TABLE]

The term 2\exp\bigl{(}-\mathrm{i}\frac{\beta\eta^{2}}{2}t-\beta\eta\left|r\right|\bigr{)} in Eq. (28) corresponds to a single bound state wavefunction of the effective potential $W_{\text{eff}}(r)$, if the condition $\Re(\beta\eta)>0$ is fulfilled. This holds, if $\left|\Delta\right|>0.8665\,\gamma$, i.e., under off-resonant driving conditions. The size of the bound state (in units of $R_{\mathrm{B}}$) is equal to

[TABLE]

In Fig. 5 we show the bound and continuum-state contributions obtained from Eq. (28) and compare them to the full numerical solution. One recognizes very good agreement, which also shows that the approximation used to derive Eq. (26) is justified. One notices that the spatial structure of bound and continuum states near $r=0$ is the same.

The complex energy of the bound states can be read off from Eq. (28). Up to second order in $\gamma/\Delta$ it is given by

[TABLE]

From this we can also read off the decay rate of the bound state which is approximately

[TABLE]

Note that both $E_{0}$ and $\gamma_{\mathrm{b}}$ are in units of $2\Omega^{2}/|\Delta|$. One recognizes that long lifetimes of bound states require small optical depth per blockade volume, $\xi$.

IV.2 Bound-state and continuum contributions to bunching

There are two distinct features of bound and continuum states. First of all, while in the vicinity of $r=0$ the continuum states have the same spatial structure as the bound state, they are the dominant contribution at large relative distances $r$, see Fig.5. This is due to the exponential localization of the bound state on a length scale $r_{\mathrm{b}}$. Secondly, as can be seen from Eq. (28), bound and continuum contributions have a different time evolution. While the continuum states decay diffusively in time, i.e., $\propto 1/\sqrt{t}$, bound states decay exponentially. This is illustrated in Fig. 6, where we have plotted the amplitudes of bound and continuum state as function of time at vanishing relative distance $r=0$ along with the two-photon amplitude $\mathcal{EE}(0,t)$. The larger the detuning $|\Delta|$ the slower the decay of the bound state. Nevertheless for large times the continuum contributions become the dominant part also for small relative distances. The oscillatory behavior of $\mathcal{EE}$ is an interference effect between bound and continuum contributions, which will be discussed in more detail later.

At large times, the solution of (28) at $r=0$, can be further simplified to

[TABLE]

where again the first term on the right hand side corresponds to the bound state while the second term gives the contribution from the continuum.

Using Eq. (28) or, for sufficiently large $\left|\Delta\right|$, Eq. (33) we can identify a crossover time $t_{0}$ at which the contribution of the scattering states becomes the dominant one. This also means that for $t\gg t_{0}$ any observed bunching is solely due to the scattering states. A simple analysis shows that $t_{0}$ is minimal if

[TABLE]

i.e., when $\gamma/\left|\Delta\right|=\tan\frac{3\pi}{16}\approx 0.6681\approx 2/3$ and at this point for $t_{0}$ one has

[TABLE]

according to Eq. (28). For $\xi<1$ this time scale is much larger than one in units of the typical EIT time scale $\left|\Delta\right|/2\Omega^{2}$.

V Filtering of bound and continuum components

As can be seen in Fig. 6 the two-photon amplitude shows an oscillatory time dependence. These oscillations result from an interference between bound- and continuum-states contributions to the two-photon amplitude, due to the different phases of these terms. We will now investigate the time evolution of the phase in more detail and will argue that this can be used to filter out the bound-state components, allowing for an experimental investigation of the photonic molecules alone.

In order to employ the phase shift of the photonic molecule for its experimental separation, it should be spatially homogeneous and at the same time sufficiently large. In Fig. 7 we show amplitude and phase of $\mathcal{EE}(r,t)$ as functions of relative distance $r$ and time $t$ obtained from numerically solving the full two-particle evolution. One recognizes that the phase shift is large and constant in space over the whole extend of the localized two-photon component for fixed times.

This is fully consistent with the approximate analytic solutions obtained in Sec.IV.

As can be seen from Eq. (33), the bound state attains a dynamical phase

[TABLE]

At the point $\gamma/\left|\Delta\right|=\tan\frac{3\pi}{16}$, where $t_{0}$ is minimal, this dynamical phase vanishes, as $\Re(\beta\eta^{2})=0$.

At the same time, the phase of the continuum states approaches for large times the constant value

[TABLE]

for $\gamma/\left|\Delta\right|=\tan\frac{3\pi}{16}$. Going to larger detunings the bound state attains a nonvanishing dynamical phase and the phase of the continuum states increases with $\left|\Delta\right|$ up to a value of

[TABLE]

for very large $\left|\Delta\right|$. This is further illustrated in Fig. 8, where we have plotted the phases of the bound and continuum states as function of time for the detuning of $\Delta\approx-\tan\frac{3\pi}{16}\gamma$ and a larger detuning of $\Delta=-12\gamma$. We have verified the accuracy of the approximate solutions by comparison to numerical solutions of the Maxwell-Bloch equations in the case $K=0$.

We note that this phase of the continuum state is very robust, as it only depends on $\left|\Delta\right|/\gamma$ and can be tuned by changing the frequencies of the probe and control fields.

When combining our results about the amplitudes and phases of the bound state and continuum state components one can distinguish three regimes depending on the ratio $\left|\Delta\right|/\gamma$. First, for small detuning the continuum states dominate the dynamics at all times and attain a phase of $-\pi/2$. Secondly, for intermediate detuning the crossover time $t_{0}$ increases and the dynamics is governed by an interplay of bound and continuum state components, and finally, for large detuning the continuum state decays very quickly and the bunching is solely due to the bound state. To observe the photonic molecule, one could simply go to large detuning and wait until the continuum states are decayed. However, this would be experimentally challenging. Therefore it is better to work in the regime of intermediate detunings. As can be seen from Fig. 8, in this case the phase of the continuum contribution attains its robust asymptotic value long before the crossover point $t_{0}$. Thus, the continuum contribution can effectively be filtered out by interferometric techniques, as sketched in Fig. 1(b), allowing for an isolation and observation of the probe field component corresponding to the molecular state.

VI Conclusion

We discussed the bunching of dark-state polaritons propagating under conditions of electromagnetically induced transparency in a gas of atomic three-level atoms and interacting via van der Waals-type interactions mediated by Rydberg interactions of the atoms. By employing a Green’s function approach, we derived an effective model for two dark-state-polariton excitations and analyzed its spectral properties, showing the existence of bound eigenstates. We showed that for weak interactions, quantified by the optical depth per blockade, and in an off-resonant driving scheme the model has a single eigenstate close to the scattering continuum. We argued that, while the higher-$n$ bound states are difficult to excite, this low-energy single bound state can experimentally observed. We confirmed this by numerical integration of the full Maxwell-Bloch equations for two particles which shows bunching for sufficiently small values of the optical depth per blockade, but anti-bunching for larger values, as has also be shown in recent experiments Peyronel et al. (2012); Firstenberg et al. (2013).

By using the Green’s function approach we showed that this bunching feature cannot solely explained by the bound eigenstate, but rather comes about by an interplay of bound and continuum states. We derived closed analytic expressions for the bound state and continuum wave functions in the limit of weak interactions, where the effective interaction potential can be approximated by a $\delta$-potential. This expressions allowed us to investigate the time-dependence of the individual components. Specifically, we showed that the bound state decays exponentially in time, whereas the scattering states have an diffusive time-dependence. Thus, for small times the bunching has to be explained by a superposition of bound and continuum wave function, while for large times the polariton pair is dominated by the continuum. Moreover, we found that, after some time, the continuum component attains a robust and constant phase, while the bound state exhibits a dynamical phase. This allows to filter bound and continuum components by making use of a homodyne detection scheme. We here concentrated on an effective one-dimensional setting. In three spatial dimensions there is an additional constraint for the existence of a bound state, which we discuss in appendix B.

Appendix A Green’s function approach

Here we present the details of the Green’s function approach employed in the derivation of Eqs. (6) to (11). The noninteracting result, Eq. (8) follows then immediately from Eq. (6) when setting the effective potential $W$ equal to zero.

To simplify the derivation in the following we consider an initial vector $\ket{\Psi_{0}}\equiv\mathbf{\Psi}(K,r,t=0)$ given by

[TABLE]

where $\ket{\varphi_{1}}=(1,0,0,0)^{T}$ denotes a pure field excitation, i.e., $\mathcal{EE}$. The derivation with other initial conditions can be treated analogously. The time evolution of the two-photon amplitude $\mathcal{EE}(t)$ for admits a spectral Fourier-Laplace representation, given by

[TABLE]

for $t>0$. The full and free Green’s functions are defined by

[TABLE]

respectively, where $\hat{\mathcal{H}}=\hat{\mathcal{H}}_{0}+V(r)\hat{\mathrm{P}}_{\mathcal{SS}}$ as in the main text. We denote the Green’s function governing the time evolution of $\mathcal{EE}(t)$ by

[TABLE]

This operator usually has a branch cut and in the presence of an interaction between atoms it may have poles, which correspond to resonant states with negative imaginary parts ($t>0$). For sufficiently large times, larger than the decay time of resonances, branch cut singularities of the Green’s function contribute into the integral (40) only.

The Green’s function (42) satisfies the operator equation

[TABLE]

where the propagators in this equation can be written in the form

[TABLE]

with the quantities

[TABLE]

the effective mass

[TABLE]

and

[TABLE]

In the next step we find the equation for the Green’s function $\bra{\varphi_{4}}\hat{G}(\omega)\ket{\varphi_{1}}$. It can be derived in an analog manner to Eq. (43) which yields

[TABLE]

where $\alpha_{00}(\omega)=\frac{\mathrm{i}\Gamma}{2\Omega^{2}-\mathrm{i}\omega\Gamma}$. Absorbing the factor $(1+\alpha_{00}(\omega)V)$ into the Green’s function $G_{41}$ by defining

[TABLE]

we can write the operator equation (51) in the closed form

[TABLE]

where the effective potential $W(r,\omega)$ is defined by

[TABLE]

Combining Eq. (53) with Eq. (43), we arrive at the following equation

[TABLE]

for the required Green’s function $G_{11}(\omega)$. The evolution of the two-photon amplitude, in the coordinate representation,

[TABLE]

is obtained by substituting Eq. (55) into the integral (40).

We note that the first integral at large $t\gg\frac{\Delta^{2}}{2\gamma g^{2}}$ is negligible. We thus have

[TABLE]

where $G(r,r^{\prime},\omega)=G_{41}(\omega)/\gamma(\omega)$, which obeys the following integral equation,

[TABLE]

As was pointed out earlier, our interests are restricted to large times. Thus we can further simplify the expressions for the Green’s functions by considering only small frequencies and momentum in Fourier space. In particular, we assume low frequencies $\omega\ll\min(2\Omega^{2}/\left|\Gamma\right|,2g^{2}/\left|\Gamma\right|)$ and small center-of-mass momentum $cK\ll\min(2\Omega^{2}/\left|\Gamma\right|,2g^{2}/\left|\Gamma\right|)$. In this limit all quantities such as $m_{0}(\omega),\Lambda_{0}(\omega)$ and $\gamma(\omega)$ take much simpler forms, and, moreover, the effective potential, Eq. (54) becomes independent of $\omega$. This leads to the following expression for $\mathcal{EE}(r,t)$,

[TABLE]

where the Green’s function $G(r,r^{\prime},\omega)$ satisfies the following integral equation

[TABLE]

The free Green’s function

[TABLE]

describes a Schrödinger particle with the complex effective mass

[TABLE]

Note that the solution of the Schrödinger problem, in integral representation, Eq. (59), will be of little help for any practical purposes. However, the Eqs. (59) to (61) allow to show that the time evolution of the two-photon amplitude $\mathcal{EE}\left(r,t\right)$ obeys the following Schrödinger equation

[TABLE]

with the initial condition

[TABLE]

Appendix B Bound states in three dimensions

In the main text we assume that the system is one-dimensional and the results are strictly valid only in this case. As in experimental setups only an approximate confinement to one dimension can be achieved we analyze the influence of higher dimensions. Therefore we consider the effective Hamiltonian corresponding to the three dimensional problem,

[TABLE]

where $k_{\mathrm{p}}=\omega_{\mathrm{p}}/c$ denotes the carrier wave number of the probe field. In order to neglect the transversal kinetic energy terms

[TABLE]

compared to the longitudinal one, the following condition should be fulfilled,

[TABLE]

where $w$ denotes the probe beam waist. We can rewrite this as a condition for the parameter $\xi$, which yields

[TABLE]

i.e., imposing a lower bound on the interaction strength.