Stable single light bullets and vortices and their active control in cold Rydberg gases

TL;DR

This paper demonstrates the generation, stability, and active control of high-dimensional light bullets and vortices in cold Rydberg gases, with potential applications in optical information processing.

Contribution

It introduces a scheme to produce stable high-dimensional optical pulses in Rydberg gases using electromagnetically induced transparency and nonlocal Kerr nonlinearities.

Findings

Stable propagation of light bullets and vortices with low power

High-efficiency storage and retrieval of optical pulses

Control over high-dimensional nonlinear optical processes

Abstract

Realizing single light bullets and vortices that are stable in high dimensions is a long-standing goal in the study of nonlinear optical physics. On the other hand, the storage and retrieval of such stable high dimensional optical pulses may offer a variety of applications. Here we present a scheme to generate such optical pulses in a cold Rydberg atomic gas. By virtue of electromagnetically induced transparency, strong, long-range atom-atom interaction in Rydberg states is mapped to light fields, resulting in a giant, fast-responding nonlocal Kerr nonlinearity and the formation of light bullets and vortices carrying orbital angular momenta, which have extremely low generation power, very slow propagation velocity, and can stably propagate, with the stability provided by the combination of local and the nonlocal Kerr nonlinearities. We demonstrate that the light bullets and vortices obtained can be stored and retrieved in the system with high efficiency and fidelity. Our study provides a new route for manipulating high-dimensional nonlinear optical processes via the controlled optical nonlinearities in cold Rydberg gases.

Equations

$${\hat{\mathcal{H}}}_H(\mathbf{r},t)$$

$$i\left(\frac{\partial }{\partial z}+\frac{1}{c}\frac{\partial }{\partial t}\right){\Omega }_p+\frac{c}{2{\omega }_p}{\nabla }_{\perp }^2{\Omega }_p+{\kappa }_{12}{\rho }_{21}(\mathbf{r},t)=0,$$

$$i\frac{\partial }{\partial t}{\rho }_{11}-i{\Gamma }_{12}{\rho }_{22}-{\Omega }_p{\rho }_{12}+{\Omega }_p^{\ast }{\rho }_{21}=0,$$

$$i\frac{\partial }{\partial t}{\rho }_{22}-i{\Gamma }_{23}{\rho }_{33}+i{\Gamma }_{12}{\rho }_{22}+{\Omega }_p{\rho }_{12}-{\Omega }_p^{\ast }{\rho }_{21}$$

$$-{\Omega }_c{\rho }_{23}+{\Omega }_c^{\ast }{\rho }_{32}=0,$$

$$i\frac{\partial }{\partial t}{\rho }_{33}+i{\Gamma }_{23}{\rho }_{33}+{\Omega }_c{\rho }_{23}-{\Omega }_c^{\ast }{\rho }_{32}=0,$$

$$\left(i\frac{\partial }{\partial t}+d_{21}\right){\rho }_{21}-{\Omega }_p({\rho }_{22}-{\rho }_{11})+{\Omega }_c^{\ast }{\rho }_{31}=0,$$

$$\left(i\frac{\partial }{\partial t}+d_{31}\right){\rho }_{31}-{\Omega }_p{\rho }_{32}+{\Omega }_c{\rho }_{21}$$

$$-{\mathcal{N}}_a\int d^3{\mathbf{r}}^{\prime }V({\mathbf{r}}^{\prime }-\mathbf{r}){\rho }_{33,31}({\mathbf{r}}^{\prime },\mathbf{r},t)=0,$$

$$\left(i\frac{\partial }{\partial t}+d_{32}\right){\rho }_{32}-{\Omega }_p^{\ast }{\rho }_{31}-{\Omega }_c({\rho }_{33}-{\rho }_{22})$$

$$-{\mathcal{N}}_a\int d^3{\mathbf{r}}^{\prime }V({\mathbf{r}}^{\prime }-\mathbf{r}){\rho }_{33,32}({\mathbf{r}}^{\prime },\mathbf{r},t)=0,$$

$$i\left(\frac{\partial }{\partial z}+{\alpha }_0\right)U-\frac{1}{2}{K}_2\frac{{\partial }^{2}U}{\partial {\tau }^2}+\frac{c}{2{\omega }_p}{\nabla }_{\perp }^{2}U+W_1|U{|}^{2}U$$

$$a_{33,31}^{(3)}({\mathbf{r}}^{\prime }-\mathbf{r})$$

$$\simeq \frac{-2|{\Omega }_c{|}^2{\Omega }_c(2\omega +d_{21}+d_{31})/|D(\omega ){|}^2}{2(\omega +d_{21})|{\Omega }_c{|}^2-D_2(\omega )[2\omega +2d_{31}-V({\mathbf{r}}^{\prime }-\mathbf{r})]},$$

$$i\frac{\partial u}{\partial s}-s_d\frac{{\partial }^{2}u}{\partial {\sigma }^2}+2u|u{|}^2=-g_{\mathrm{diff}}\left(\frac{{\partial }^2}{\partial {\xi }^2}+\frac{{\partial }^2}{\partial {\eta }^2}\right)u+i{d}_{0}u,$$

$${\Omega }_p=\frac{1}{{\tau }_0}\sqrt{\frac{\widetilde{|K_2|}}{\tilde{W}}}\mathrm{sech}\left[\frac{1}{{\tau }_0}\left(t-\frac{z}{\widetilde{V_g}}\right)\right]e^{i{\tilde{K}}_{0}z-iz/2L_{\mathrm{D}}},$$

$$i\frac{\partial u}{\partial s}+\left(\frac{{\partial }^2}{\partial {\xi }^2}+\frac{{\partial }^2}{\partial {\eta }^2}\right)u+g_d\frac{{\partial }^{2}u}{\partial {\sigma }^2}+g_1|u{|}^{2}u$$

$$u=$$

$$\times \exp \left[-i{C}_s(s)\frac{{\xi }^2+{\eta }^2}{2w_s^2(s)}-i{C}_t(s)\frac{{\sigma }^2}{2}+i\varphi (s)\right],$$

$$u_{mp}=$$

$$i\frac{\partial u}{\partial s}+\left(\frac{{\partial }^2}{\partial {\xi }^2}+\frac{{\partial }^2}{\partial {\eta }^2}\right)u+g_d\frac{{\partial }^2}{\partial {\sigma }^2}u+g_1|u{|}^{2}u$$

$$-g_4({\xi }^2+{\eta }^2)u=i{d}_{0}u,$$

$$V_{\mathrm{LB}}\approx 1.25\times 10^{-6}c,$$

$${\bar{P}}_{\max }\approx 0.2\mathrm{nW},$$

$$a_{21}^{(2)}=\frac{i}{{\kappa }_{12}}\left(\frac{1}{V_{\mathrm{g}}}-\frac{1}{c}\right),$$

$$a_{31}^{(2)}=-\frac{i}{{\Omega }_c^{\ast }}\frac{\omega +d_{31}}{D(\omega )}-\frac{\left(\omega +d_{21}\right)}{{\Omega }_c^{\ast }}{a}_{21}^{(2)},$$

$$a_{11}^{(2)}=\frac{[i{\Gamma }_{23}-2|{\Omega }_c{|}^2{M}_1]M_2-i{\Gamma }_{12}|{\Omega }_c{|}^2{M}_3}{-{\Gamma }_{12}{\Gamma }_{23}-i{\Gamma }_{12}|{\Omega }_c{|}^2{M}_1},$$

$$a_{33}^{(2)}=\frac{1}{i{\Gamma }_{12}}\left(M_2-i{\Gamma }_{12}{a}_{11}^{(2)}\right),$$

$$a_{32}^{(2)}=\frac{1}{d_{32}}\left(-\frac{{\Omega }_c}{D(\omega )}+2{\Omega }_c{a}_{33}^{(2)}+{\Omega }_c{a}_{11}^{(2)}\right),$$

$$\left[\begin{array}{ccc}2\omega +2d_{21}& 0& 2{\Omega }_c^{\ast }\\ 0& 2\omega +2d_{31}-V& 2{\Omega }_c\\ {\Omega }_c& {\Omega }_c^{\ast }& 2\omega +d_{21}+d_{31}\end{array}\right]\left[\begin{array}{c}{\rho }_{21,21}^{(2)}\\ {\rho }_{31,31}^{(2)}\\ {\rho }_{31,21}^{(2)}\end{array}\right]$$

$$=\left[\begin{array}{c}-2\frac{\omega +d_{31}}{D(\omega )}\\ 0\\ \frac{{\Omega }_c}{D(\omega )}\end{array}\right]F^2{e}^{2i\theta },$$

$$\left[\begin{array}{cccc}{d}_{21}+d_{12}& 0& -{\Omega }_c& {\Omega }_c^{\ast }\\ -{\Omega }_c^{\ast }& {\Omega }_c^{\ast }& d_{21}+d_{13}& 0\\ 0& d_{31}+d_{13}& {\Omega }_c& -{\Omega }_c^{\ast }\\ -{\Omega }_c& {\Omega }_c& 0& d_{21}^{\ast }+d_{13}^{\ast }\end{array}\right]\left[\begin{array}{c}{\rho }_{21,12}^{(2)}\\ {\rho }_{31,13}^{(2)}\\ {\rho }_{21,13}^{(2)}\\ {\rho }_{21,13}^{\ast (2)}\end{array}\right]$$

$$=\left[\begin{array}{c}\frac{\omega +d_{31}}{D}-\frac{\omega +d_{31}^{\ast }}{D(\omega {)}^{\ast }}\\ \frac{{\Omega }_c^{\ast }}{D(\omega {)}^{\ast }}\\ 0\\ \frac{{\Omega }_c}{D(\omega )}\end{array}\right]|F{|}^2{e}^{-2\bar{\alpha }{z}_2}.$$

$$\left[\begin{array}{c}\begin{array}{cccccccc}{M}_{31}& {\Omega }_c^{\ast }& -i{\Gamma }_{23}& 0& {\Omega }_c^{\ast }& -{\Omega }_c& 0& 0\\ {\Omega }_c& M_{32}& 0& -i{\Gamma }_{23}& 0& 0& {\Omega }_c^{\ast }& -{\Omega }_c\\ 0& 0& M_{33}& {\Omega }_c^{\ast }& -{\Omega }_c^{\ast }& {\Omega }_c& 0& 0\\ 0& 0& {\Omega }_c& M_{34}& 0& 0& -{\Omega }_c^{\ast }& {\Omega }_c\\ {\Omega }_c& 0& -{\Omega }_c& 0& M_{35}& 0& {\Omega }_c^{\ast }& 0\\ -{\Omega }_c^{\ast }& 0& {\Omega }_c^{\ast }& 0& 0& M_{36}& 0& {\Omega }_c^{\ast }\\ 0& {\Omega }_c& 0& -{\Omega }_c& {\Omega }_c& 0& M_{37}& 0\\ 0& -{\Omega }_c^{\ast }& 0& {\Omega }_c^{\ast }& 0& {\Omega }_c& 0& M_{38}\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{\rho }_{22,21}^{(3)}\\ {\rho }_{22,31}^{(3)}\\ {\rho }_{33,21}^{(3)}\\ {\rho }_{33,31}^{(3)}\\ {\rho }_{32,21}^{(3)}\\ {\rho }_{21,23}^{(3)}\\ {\rho }_{32,31}^{(3)}\\ {\rho }_{31,23}^{(3)}\end{array}\end{array}\right]$$

$$=\left[\begin{array}{c}\begin{array}{c}-a_{21,12}^{(2)}+a_{21,21}^{(2)}-a_{22}^{(2)}\\ -a_{31,12}^{(2)}+a_{21,31}^{(2)}\\ -a_{33}^{(2)}\\ 0\\ a_{21,31}^{(2)}-a_{32}^{(2)}\\ -a_{32}^{\ast (2)}-a_{21,13}^{(2)}\\ a_{31,31}^{(2)}\\ -a_{31,13}^{(2)}\end{array}\end{array}\right]|F({\mathbf{r}}^{\prime }){|}^{2}F(\mathbf{r})e^{-2\bar{\alpha }{z}_2^{\prime }}{e}^{i\theta },$$

$${\rho }_{33,31}^{(3)}=a_{33,31}^{(3)}|F({\mathbf{r}}^{\prime }){|}^{2}F(\mathbf{r})e^{-2\bar{\alpha }{z}_2^{\prime }}{e}^{i\theta },$$

$$a_{33,31}^{(3)}$$

$$\simeq \frac{-2|{\Omega }_c{|}^2{\Omega }_c(2\omega +d_{21}+d_{31})/|D(\omega ){|}^2}{2(\omega +d_{21})|{\Omega }_c{|}^2-D_2(\omega )[2\omega +2d_{31}-V({\mathbf{r}}^{\prime }-\mathbf{r})]}.$$

Full text

††preprint:

I INTRODUCTION

High-dimensional spatiotemporal optical solitons, alias light bullets (LBs) Silberberg , are solitary nonlinear waves localized in $n$-spatial dimensions and one time axis [($m$+1)D; $m=1,2,3$]. In recent years, the study of LBs has attracted intensive theoretical and experimental interests Kivshar because of their rich nonlinear physics and technological applications Malomed ; Mihalache ; Belic ; Miha ; Mateo ; Karta . Experimental signatures of LBs in different types of optical media, such as (2+1)D LBs in LiIO3 crystals Liu and quasi-(3+1)D LBs in waveguide arrays Minardi ; Eile , have been reported. Recently, trains of (3+1)D dark solitons have been observed in a photo-refractive material by Lahav et al. Lah . These LBs, however, generally suffer severe instability, which typically can propagate to only a few diffraction lengths. To levitate the rapid distortion of LBs, short pulses (in the order of femtosecond) with high laser powers were commonly used experimentally. Despite fast experimental progresses, it remains a challenge to create stable single LBs, as they can only be realized by carefully balancing their dispersion, diffraction and optical nonlinearities.

Theoretical efforts have attempted to study the formation of LBs with different mechanisms. A commonly used approach for creating stable (3+1)D LBs is to exploit local and nonlocal optical nonlinearities with quite different response times. These nonlinearities have been examined, e.g., in nematic liquid crystals Bur ; Pec1 ; Pec2 and lead glass Gur , in which the nonlocal nonlinearity (resulted from the reorientation Bur ; Pec1 ; Pec2 or thermal motion Gur of molecules) has a very long response time (typically in the order of second or even longer), while the local nonlinearity (resulted from the electronic Kerr effect) has a very short response time (in the order of femtosecond). Due to the mismatch of the response-time scales and the fast propagation velocity ($\sim c$) of optical pulses, stable LBs can only be realized in the form of pulse-train solitons (not single LBs) Lah . A different approach is based on electromagnetically induced transparency (EIT) in atomic gases FIM , in which optical absorption can be largely suppressed due to quantum interference. Together with Kerr nonlinearities Hau ; Pri induced by resonant laser fields, it has been shown that stable (1+1)D temporal Wu ; huang ; Chen and spatial Hong ; Michinel optical solitons can form. Though promising, the local nature of the Kerr nonlinearity doesn’t support stable LBs, which suffer unavoidable catastrophic collapse during propagation.

On the other hand, recent theoretical Friedler and experimental Moha studies revealed that strong and long-range optical nonlinearities can be built with Rydberg atomic gases REITreview ; Fir ; Mur , which are in electronically high-lying states with large principal quantum number $n$ Gallagher . Large dipole moments in Rydberg states render a strong, long-range Rydberg-Rydberg interaction between remote atoms. Rydberg interactions find applications in quantum information, precision measurement, quantum simulation based on Rydberg atoms Saffman ; Low ; REITreview ; Fir ; Mur . Importantly, the Rydberg-Rydberg interaction can be mapped to a nonlocal optical nonlinearity through EIT, which is strong even at the single photon level REITreview ; Fir ; Mur . Long lifetimes (in the order of tens of microsecond) guarantee that the induced photon-photon interaction is largely coherent during light propagation Gor0 ; BingHe ; LiuYang . This provides a new platform to study quantum nonlinear optics Fir and develop new photon devices, such as single-photon switches and transistors Baur ; Gorn1 ; Gorn2 , quantum memories and phase gates Maxwell2013 ; Distante ; LiL ; Busche ; Tia2 ; Liang .

In this work, we present a scheme for the generation and storage of stable (3+1)D LBs in highly tunable cold gases of Rydberg atoms. A key element is the co-existence of giant nonlocal and local optical nonlinearities. The former features a fast (sub-microsecond) response Zhang , which is complemented by the latter thomas ; Gullans , whose response is relative slower (in the order of microsecond). In conjunction with tunable dispersion and diffraction, this allows us to precisely control dynamics of LBs. To go beyond the commonly used mean-field theory Sche ; Sev ; Stan ; Bai , we derive a nonlinear (3+1)D light propagation equation taking into account of one-body and two-body correlations. We show that the synergetic nonlocal and local optical nonlinearities in the system permit us to obtain stable single (3+1)D LBs as well as LBs that carry definite orbital angular momentum [i.e. light vortices (LVs) ] Andrews ; Veissier ; Nicolas ; Parigi ; Oliveira ; ding . It is revealed that the stability of the single LBs and LVs is achieved via a two-step self-trapping mechanism: First, the transverse self-trapping of an input high-dimensional laser pulse is quickly reached through the balance between the diffraction and the fast nonlocal nonlinearity; then the laser pulse is further self-trapped in the longitudinal (propagation) direction through the interplay between the dispersion and the slow local nonlinearity. As a result, stable bright (3+1)D LBs can be generated with ultraslow propagation velocity, extremely low light power, and narrow bandwidth. More importantly, based on the active character of the system, we can manipulate the LBs and LVs by tuning system parameters, and can even store and retrieve LBs and LVs in the Rydberg medium with high efficiency and fidelity.

The paper is arranged as follows. In Sec. II, we describe the physical model of the Rydberg-EIT system. In Sec. III, we derive the (3+1)D nonlocal nonlinear envelope equation based on the multiple-scale method. In Sec. IV, we discuss the physical property of the nonlocal and local Kerr nonlinearities and present LB and LV solutions in various regimes. Based on numerical calculations, we illustrate in Sec. V that LBs can be used as optical memories such that they can be dynamically stored and read out in the Rydberg medium. Finally, in the last section (Sec. VI) we summarize the main results obtained in this work.

II Model

We consider a three-level atomic system with a ladder-type level configuration, shown in Fig. 1(a).

A weak, pulsed probe field of angular frequency $\omega_{p}$ (half Rabi frequency $\Omega_{p}$) couples the ground state $|1\rangle$ and intermediate state $|2\rangle$, and a strong, continuous-wave control field of angular frequency $\omega_{c}$ (half Rabi frequency $\Omega_{c}$) couples state $|2\rangle$ and a Rydberg state $|3\rangle$. The probe field has a pulse length $\tau_{0}$ at the entrance of the medium. The electric-field vector of the system is ${\bf E}({\bf r},t)=\sum_{l=p,c}{\bf e}_{l}\,\mathcal{E}_{l}\,\exp[i({\bf k}_{l}\cdot{\bf r}-\omega_{l}t)]+{\rm c.c.}$, where ${\bf e}_{l}$ $({\bf k}_{l})$ is the unit polarization vector (wavevector) of the electric-field component with envelope $\mathcal{E}_{l}\,\,(l=p,c)$. A possible experimental arrangement of beam geometry is shown in Fig. 1(b).

The system works at a ultracold temperature and the probe and control fields counter-propagate [i.e., ${\bf k}_{p}=(0,0,k_{p})$, ${\bf k}_{c}=(0,0,-k_{c})$], so that the center-of-mass motion of the atoms and dephasing due to the atomic collisions are negligible. Under the electric-dipole approximation, the system Hamiltonian is $\hat{H}_{H}={\cal N}_{a}\int{\rm d}^{3}\mathbf{r}\,\hat{{\cal H}}_{H}(\mathbf{r},t)$ with

[TABLE]

where $\mathcal{N}_{a}$ is atomic density, $\hbar\omega_{\alpha}$ is the eignenergy of atomic state $|\alpha\rangle$, $\Omega_{p}=(\mathbf{e}_{p}\cdot\mathbf{p}_{21})\mathcal{E}_{p}/\hbar$ and $\Omega_{c}=(\mathbf{e}_{c}\cdot\mathbf{p}_{32})\mathcal{E}_{c}/\hbar$ are respectively the half Rabi frequencies of the probe and control fields (with $\mathbf{p}_{\alpha\beta}$ the electric dipole matrix element associated with the transition from $|\beta\rangle$ to $|\alpha\rangle$), $\hat{S}_{\alpha\beta}$ are transition operators $(\alpha,\beta=1,2,3)$ satisfying the commutation relation $\left[\hat{S}_{\alpha\beta}(\mathbf{r},t),\hat{S}_{\mu\nu}(\mathbf{r}^{\prime},t)\right]=\left(\delta_{\alpha\nu}\hat{S}_{\mu\beta}(\mathbf{r},t)-\delta_{\mu\beta}\hat{S}_{\alpha\nu}(\mathbf{r}^{\prime},t)\right)\delta_{\mathbf{r}\mathbf{r}^{\prime}}$. In the last term on the right hand side (RHS), $\hbar V({\bf r}-{\bf r}^{\prime})=-\hbar C_{6}/|{\bf r}-{\bf r}^{\prime}|^{6}$ is the van der Waals (vdW) interaction potential with dispersion coefficient $C_{6}$ between the Rydberg atoms located at the positions ${\bf r}$ and ${\bf r}^{\prime}$, respectively.

To be concrete, we will consider strontium atoms (88Sr) in this work, although our approach is valid for other Rydberg atomic gases. The energy-levels shown in Fig. 1(a) are selected as $|1\rangle=|5s^{2}\,{}^{1}S_{0}\rangle$, $|2\rangle=|5s5p^{1}P_{1}\rangle$, $|3\rangle=|5sns^{1}S_{0}\rangle$, with $n$ the principle quantum number Mauger . The spontaneous emission rates of the atoms are given by $\Gamma_{12}=2\pi\times 16\,\,{\rm MHz},$ $\Gamma_{23}=2\pi\times 16.7\,{\rm kHz}$, so one has $\gamma_{21}=\Gamma_{12}/2$, $\gamma_{31}=\Gamma_{23}/2$, $\gamma_{32}=(\Gamma_{12}+\Gamma_{23})/2$. For this choice, the vdW interaction in ${}^{1}S_{0}$ states is isotropically attractive ($C_{6}>0$), which is important to realize a self-focusing nonlinearity.

To study propagation of the probe field, we assume the size of the atomic gas is much larger than the Rydberg blockade radius (given later). Depletion of the control field will be neglected (except in the discussion of LB memories in Sec. V). Additionally assuming that the weak probe pulse is still a classical field (consisting of many photons), so that a semi-classical approach can be used. With the slowly varying envelope and rotating wave approximations, propagation of the probe field is governed by the coupled Maxwell-Bloch equations FIM ,

[TABLE]

where $\nabla_{\bot}^{2}=\partial^{2}/\partial x^{2}+\partial^{2}/\partial y^{2}$, $\kappa_{12}={\cal N}_{a}\omega_{p}|\mathbf{p}_{12}|^{2}/(2\epsilon_{0}c\hbar)$ with $\epsilon_{0}$ the vacuum dielectric constant. $\rho_{21}({\bf r},t)\equiv\langle{\hat{S}}_{21}({\bf r},t)\rangle$ is the coherence between level $|1\rangle$ and $|2\rangle$. Its dynamics is determined by the Bloch equations,

[TABLE]

where $\rho_{\alpha\beta}({\bf r},t)\equiv\langle{\hat{S}}_{\alpha\beta}({\bf r},t)\rangle$ are the one-body correlators (one-body density matrix elements)  note00 , $d_{\alpha\beta}=\Delta_{\alpha}-\Delta_{\beta}+i\gamma_{\alpha\beta}$ ($\Delta_{1}=0$; $\alpha,\beta=1,2,3;\alpha\neq\beta)$, $\Delta_{2}=\omega_{p}-(\omega_{2}-\omega_{1})$ and $\Delta_{3}=\omega_{p}+\omega_{c}-(\omega_{3}-\omega_{1})$ are respectively the one-photon and two-photon detunings, $\gamma_{\alpha\beta}=(\Gamma_{\alpha}+\Gamma_{\beta})/2+\gamma_{\alpha\beta}^{\rm col}$ with $\Gamma_{\beta}=\sum_{\alpha<\beta}\Gamma_{\alpha\beta}$. Here $\Gamma_{\alpha\beta}$ denotes respectively the spontaneous emission decay rate from the state $|\beta\rangle$ to the state $|\alpha\rangle$, and $\gamma_{\alpha\beta}^{\rm col}$ represents the dephasing rate reflecting the loss of phase coherence between $|\alpha\rangle$ and $|\beta\rangle$ due to, e.g., atomic motion and the interaction between the atoms in the ground and Rydberg states.

The last terms on the left hand side of Eq. 2(e) and Eq. 2(f), i.e. the two-body correlators $\rho_{33,3\alpha}(\mathbf{r}^{\prime},\mathbf{r},t)\equiv\langle\hat{S}_{33}(\mathbf{r}^{\prime},t)\hat{S}_{3\alpha}(\mathbf{r},t)\rangle$ ($\alpha=1,2$), result from the interaction between Rydberg atoms, where attractive (repulsive) atomic interaction lead to nonlocal self-focusing (self-defocuing) nonlinearities. Besides, there is another type of (local) nonlinearity when two-photon detuning is non-zero (i.e. $\Delta_{3}\neq 0$) Bai , contributed solely from the resonant coupling between the probe field and atoms.

From Eq. (2), we see that to solve the one-body correlators $\rho_{\alpha\beta}({\bf r},t)$, we need to know values of the two-body correlators $\rho_{33,3\alpha}({\bf r}^{\prime},{\bf r},t)\equiv\langle{\hat{S}}_{33}({\bf r}^{\prime},t){\hat{S}}_{3\alpha}({\bf r},t)\rangle$ ($\alpha=1,2$), for which three-body correlators, i.e., $\rho_{\alpha\beta,\mu\nu,\zeta\eta}({\bf r}^{\prime\prime},{\bf r}^{\prime},t)\equiv\langle{\hat{S}}_{\alpha\beta}({\bf r}^{\prime\prime},t){\hat{S}}_{\mu\nu}({\bf r}^{\prime},t){\hat{S}}_{\zeta\eta}({\bf r},t)\rangle$, will be involved, etc. As a result, we obtain an infinite hierarchy of equations of motion for the correlators of one-body, two-body, three-body, and so on. Therefore approximations have to be used in order to truncated the hierarchy.

III The (3+1)D nonlinear envelope equation

Due to the nonlinear coupling originated by the Rydberg-Rydberg interaction, it is difficult to solve the hierarchy of the equations for many-body correlators by employing conventional techniques. Fortunately, since in our consideration the probe field is relatively weak, we can solve the hierarchy of equations by using the method of multiple-scales huang ; Chen ; Jeff widely applied for solving nonlinear wave problems Jeff . Because our calculation is exact to third order (i.e. up to $\Omega_{p}^{3}$), higher-order $n$-body correlators ($n\geq 3$) involved in the hierarchy play no significant role and hence can be neglected.

We assume that the atoms are initially populated in state $|1\rangle$ and make the asymptotic expansion huang ; Chen : $\rho_{\alpha 1}=\sum_{l=0}\epsilon^{2l+1}\rho_{\alpha 1}^{(2l+1)}$, $\rho_{32}=\sum_{l=1}\epsilon^{2l}\rho_{32}^{(2l)}$, $\rho_{\beta\beta}=\sum_{l=0}\epsilon^{2l}\rho_{\beta\beta}^{(2l)}$ with $\rho_{\beta\beta}^{(0)}=\delta_{\beta 1}\delta_{\beta 1}$($\alpha=2,3;\beta=1,2,3$), $\Omega_{p}=\epsilon\Omega_{p}^{(1)}+\epsilon^{2}\Omega_{p}^{(2)}+\epsilon^{3}\Omega_{p}^{(3)}+\cdots$, where $\epsilon$ a parameter characterizing the order of $\Omega_{p}$. To obtain a divergence-free expansion, all quantities on the right hand side of the expansion are considered as functions of the multiscale variables $z_{l}=\epsilon^{l}z$ $(l=0,1,2)$, $(x_{1},y_{1})=\epsilon\,(x,y)$, and $t_{l}=\epsilon^{l}t$ $(l=0,1)$. Substituting the above expansion into the equations of the one-body and two-body correlators (omitted here for saving space) and the Maxwell Eq. (1), and comparing the coefficients of $\epsilon^{l}$ $(l=1,2,3...)$, we obtain a set of linear but inhomogeneous equations which can be solved order by order.

At the first order, we obtain the solution $\Omega_{p}^{(1)}=F\,{\rm exp}[i(Kz_{0}-\omega t_{0})]$, $\rho_{21}^{(1)}=(\omega+d_{31})\Omega_{p}^{(1)}/D(\omega)$, and $\rho_{31}^{(1)}=-\Omega_{c}\Omega_{p}^{(1)}/D(\omega)$, with $D(\omega)=|\Omega_{c}|^{2}-(\omega+d_{21})(\omega+d_{31})$ and other $\rho_{\alpha\beta}^{(1)}$=0 note0 . In above expressions, $F=F(x_{1},y_{1},z_{1},t_{1},z_{2})$ is an envelope function to be determined yet, and $K(\omega)$ is the linear dispersion relation, given by $K(\omega)=\omega/c+\kappa_{12}(\omega+d_{31})/D(\omega)$. At the second order, a solvability (i.e. divergence-free) condition requires $\partial F/\partial z_{1}+(1/V_{g})\partial F/\partial t_{1}=0$, with $V_{g}=(\partial K/\partial\omega)^{-1}$ the group velocity of the envelope $F$. The explicit expression of the second-order approximate solution for the one-body correlators, i.e. $\rho_{\alpha\beta}^{(2)}$, is presented in Appendix A. Expansion equations for the two-body correlators are given in Appendix B.

With the first-order and the second-order solutions given above, we can go to the third order approximation. A solvability condition at this order yields the (3+1)D equation for the envelope function $F$:

[TABLE]

after returning to original variables, where ${\bf r}_{\bot}=(x,y,0)$, $U=\epsilon F$exp$(-\alpha_{0}z)$ (with $\alpha_{0}={\rm Im}(K)$ a decay constant) $\tau=t-z/V_{g}$ is travelling coordinate, $K_{2}=\partial^{2}K/\partial\omega^{2}$ describes group velocity dispersion, $W_{1}=\kappa_{12}[\Omega_{c}^{\ast}a_{32}^{(2)}+(\omega+d_{31})(2a_{11}^{(2)}+a_{33}^{(2)})]/D(\omega)$ is the Kerr coefficient describing the local nonlinear optical response contributed from the photon-atom interaction (which vanishes if the two-photon detuning $\Delta_{3}$ is taken to be zero). Explicit expressions of $a_{11}^{(2)}$, $a_{32}^{(2)}$, and $a_{33}^{(2)}$ are presented in Appendix A. $G_{2}({\bf r_{\bot}^{\prime}}-{\bf r_{\bot}})=(\kappa_{12}N_{\alpha}\Omega_{c}^{\ast})/D(\omega)\int_{-\infty}^{+\infty}a_{33,31}^{(3)}(\mathbf{r}^{\prime}-\mathbf{r})V(\mathbf{r}^{\prime}-\mathbf{r})dz^{\prime}$ (the kernel in the integral) is a reduced effective interaction potential describing the nonlocal nonlinear optical response contributed from the Rydberg-Rydberg interaction. The expression of $a_{33,31}^{(3)}(\mathbf{r}^{\prime}-\mathbf{r})$ is given by

[TABLE]

where $D_{2}(\omega)=(\omega+d_{21})(2\omega+d_{21}+d_{31})-|\Omega_{c}|^{2}$. The calculation detail for obtaining $a_{33,31}^{(3)}(\mathbf{r}^{\prime}-\mathbf{r})$ is presented in Appendix B. For simplicity, we have assumed that the spatial length of the probe pulse in the propagation (i.e., $z$) direction is much larger than the range of Rydberg-Rydberg interaction, so that a local approximation along the $z$ direction can be made Sev .

In Eq. (3) the second and third terms represent respectively the dispersion and diffraction, the fourth and fifth terms are the local and nonlocal nonlinearities of the system. Because $W_{1}\approx 0$ when $\Delta_{3}=0$, it is necessary to take a non-zero $\Delta_{3}$ for obtaining a finite local Kerr nonlinearity. It is just the synthetic Kerr nonlinearities (i.e. the local and nonlocal ones combined together) that make the system support stable LBs, as shown below.

IV Ultraslow weak light bullets and vortices

The form of the solution of the (3+1)D envelope Eq. (3) depends heavily on the property of the nonlocal nonlinearity, which is characterized by the nonlocality degree of the system. The nonlocality degree can be described by using the parameter $R_{b}/R_{0}$. Here $R_{b}$ is the spatial width of the effective interaction potential $G_{2}({\bf r_{\bot}}-{\bf r}_{\bot}^{\prime})$, i.e. the Rydberg blockade radius, defined by $R_{b}=[|C_{6}d_{21}|/(2|\Omega_{c}|^{2})]^{1/6}$ Sev ; $R_{0}$ is the transverse spatial width of the probe-pulse envelope $U$, which can be measured by the transverse beam radius of the incident probe pulse. For example, a Gaussian-type incident pulse $U=U_{0}\,{\rm exp}[-(r_{\bot}/R_{0})^{2}]$, with $r_{\bot}=[x^{2}+y^{2}]^{1/2}$. Hence, the degree of nonlocality can be easily tuned by adjusting $R_{0}$ or $R_{b}$.

In the following, we take the system parameters in the dispersive regime (i.e., $|\Delta_{2}|\gg\Gamma_{12}$). Exact values of parameters will be given later. This allows us to divide the degree of nonlocality of the system into three typical regions Sev ; Mur : (i) $R_{b}/R_{0}\ll 1$ (local response region); (ii) $R_{b}/R_{0}\sim 1$ (nonlocal response region); (iii) $R_{b}/R_{0}\gg 1$ (strong nonlocal response region), for which different types of LB solutions can be obtained.

IV.1 Local response region

We first consider the case when the range of the Rydberg-Rydberg interaction (equivalently the spatial width of the effective atomic interaction potential $G_{2}$) is much smaller than the that of the beam radius of the probe pulse), i.e. $R_{b}/R_{0}\ll 1$. Fig. 2(a)

shows Re[$G_{2}({\bf r_{\bot}^{\prime}}-{\bf r_{\bot}})$] (real part of $G_{2}$), Im[$G_{2}({\bf r_{\bot}^{\prime}}-{\bf r_{\bot}})$] (imaginary part of $G_{2}$), and $|U/U_{0}|^{2}$ (the intensity profile of the probe field) as functions of $r^{\prime}_{\bot}/R_{b}$ by taking $R_{0}=300{\rm\mu m}$ (and hence $R_{b}/R_{0}=0.019$). From the figure we see that $|U|^{2}$ varies very slowly compared with $G_{2}$, and hence the integral $\int{d^{2}\mathbf{r_{\bot}^{\prime}}G_{2}({\bf r}_{\bot}^{\prime}-{\bf r}_{\bot})|U(\mathbf{r_{\bot}^{\prime}},z,\tau)|^{2}}U(\mathbf{r},\tau)$ can be approximated by $[\int d^{2}{\bf r}_{\bot}^{\prime}G_{2}({\bf r}_{\bot}^{\prime}-{\bf r}_{\bot})]\,|U({\bf r},\tau)|^{2}U(\mathbf{r},\tau)=W_{2}|U(\mathbf{r},\tau)|^{2}U(\mathbf{r},\tau)$. Here $W_{2}=\int{d^{2}\mathbf{r_{\bot}^{\prime}}G_{2}(\bf{r_{\bot}^{\prime}})}=P/Q$, with $P\simeq-i4\pi^{2}\kappa_{12}{\cal N}_{a}\sqrt{C_{6}}|\Omega_{c}|^{4}(2\omega+d_{21}+d_{31})$, $Q\simeq 3D(\omega)|D(\omega)|^{2}\sqrt{D_{3}(\omega)}$, and $D_{3}(\omega)=D_{2}[2(\omega+d_{21}|\Omega_{c}|^{2})-D_{2}(2\omega+d_{31})]$. That is to say, for large probe-field beam width the nonlocal Kerr effect contributed by the Rydberg-Rydberg interaction reduces into a local Kerr nonlinearity. As a result, Eq. (3) is simplified into a local nonlinear Schrödinger (NLS) equation with the Kerr coefficient given by $W_{1}+W_{2}$.

Under the EIT condition ($|\Omega_{c}|^{2}>\gamma_{21}\gamma_{31}$) and a large one-photon detuning $\Delta_{2}$, imaginary parts of the coefficients in Eq. (3) are very small. Taking these small imaginary parts as perturbations, Eq. (3) can be written in the dimensionless form

[TABLE]

with $s=z/(2L_{D})$, $\sigma=\tau/\tau_{0}$, $(\xi,\eta)=(x,y)/R_{0}$, $u=U/U_{0}$, $g_{\rm diff}=L_{\rm D}/L_{\rm diff}$, $d_{0}=-2L_{D}/L_{A}$, and $s_{d}={\rm sgn}(K_{2})=\pm 1$. Here $L_{\rm diff}\equiv(\omega_{p}R_{0}^{2})/c$, $L_{D}\equiv{\tau_{0}}^{2}/|\tilde{K}_{2}|$ and $L_{A}\equiv 1/\alpha_{0}$ the typical diffraction length, dispersion length, and absorption length, respectively. Note that we have taken $L_{D}=L_{\rm NL}$ [with $L_{\rm NL}\equiv 1/(U_{0}^{2}\tilde{W_{1}}+U_{0}^{2}\tilde{W_{2}})$ being a typical nonlinear length], i.e., a balance of dispersion and nonlinearity is assumed to favor the formation of solitons, thus the typical Rabi frequency of the probe field is given by $U_{0}\equiv(1/\tau_{0})[|\tilde{K}_{2}/(\tilde{W}_{1}+\tilde{W}_{2})|]^{1/2}$. The tilde above corresponding quantities means taking their real parts. Due to the EIT effect and large $\Delta_{2}$, we have $d_{0}\ll 1$, hence the dissipation in the system plays a negligible role; furthermore, in the present local response region the system has a large diffraction length $L_{\rm diff}$ ($\gg L_{\rm D}$) due to large $R_{0}$, thus we have $g_{\rm diff}\ll 1$, which means that the diffraction in the system can be neglected. Ignoring the terms on RHS of Eq. (5) and converting to the original variables, we obtain the bright soliton solution

[TABLE]

if $K_{2}<0$ (i.e. $s_{d}=-1$), with $\tilde{K}_{0}=\tilde{K}|_{\omega=0}$.

We now give a practical example for the formation of the optical soliton given above. Choosing the parameters the same as those used in Fig. 2(a), we obtain the numerical values of the coefficients in Eq. (5), given by $K_{2}=(-1.03+0.06i)\times 10^{-13}{\rm cm}^{-1}{\rm s}^{2}$, $W_{1}=(6.86+0.062i)\times 10^{-18}{\rm cm}^{-1}{\rm s}^{2}$, $W_{2}=(1.79+0.0076i)\times 10^{-14}{\rm cm}^{-1}{\rm s}^{2}$. We see that, as expected, the imaginary parts of these coefficients are indeed much smaller than their real parts. By taking $\tau_{0}=1.2\times 10^{-7}$ s, we obtain $U_{0}=2\times 10^{7}$ s*-1*, $L_{D}=L_{\rm NL}=1.4$ mm, $L_{A}=560$ mm, and $d_{0}=-0.005$, which means that the dissipation effect of the system is indeed small. In addition, because the beam radius of the probe pulse has taken to be large ($R_{0}=300{\rm\mu m}$), one has $L_{\rm diff}=1226$ mm and hence $g_{\rm diff}=0.001$, so the diffraction effect [i.e. the first term on the RHS of Eq. (5)] in the system can indeed be neglected.

Now let’s discuss properties of the soliton. With the system parameters it is easy to obtain $\tilde{V_{g}}=9.8\times 10^{-5}c$, i.e. the soliton obtained has an ultraslow propagation velocity compared with $c$, the light speed in vacuum. Furthermore, the third-order nonlinear optical susceptibility, given by $\chi_{p}^{(3)}=2c|\mathbf{p}_{12}|^{2}(W_{1}+W_{2})/(\hbar^{2}\omega_{p})$, is estimated to have the value $(3.03+0.0129i)\times 10^{-8}\,{\rm m}^{2}\,{\rm V}^{-2}$, which is 11 orders of magnitude higher than that obtained in conventional nonlinear optical media Boyd . The physical reasons for such large third-order nonlinear optical susceptibility are due to the strong Rydberg-Rydberg interaction and the EIT effect in the system. By using Poynting’s vector huang ; Newell , the maximum average power density to generate such ultraslow optical soliton is estimated to be $\bar{P}_{\rm max}=1.2\,{\rm\mu W}$. Thus, very low input power is needed for generating the optical soliton in the system.

IV.2 Nonlocal response region

We now turn to consider nonlocal response region, corresponding to $R_{b}/R_{0}\sim 1$, which can be achieved by decreasing the probe-beam radius $R_{0}$ or increasing the principle quantum $n$ (thus $R_{b}$). One example is shown in Fig. 2(b). From the figure we see that $|U|^{2}$ varies as the same way as $G_{2}$. In this case, the nonlinear optical response contributed by the Rydberg-Rydberg interaction depends not only on a particular point on the light intensity, but also on a certain neighborhood of this point.

To seek possible LBs in the system, we write Eq. (3) into the dimensionless form

[TABLE]

through new scalings $s=z/(2L_{\rm diff})$, $(\xi,\eta)=(x,y)/R_{0}$, $g_{d}=-L_{\rm diff}\tilde{K}_{2}/\tau_{0}^{2}$, $d_{0}=-2L_{\rm diff}/L_{A}$, $g_{1}=2\tilde{W}_{1}|U_{0}|^{2}L_{\rm diff}$, and $N(\xi,\eta,s)=2L_{\rm diff}R_{0}^{2}|U_{0}|^{2}\tilde{G}_{2}(\xi,\eta)$. Note that in Eq. (7) there are four main effects, i.e., the diffraction, dispersion, local nonlinearity, and the nonlocal nonlinearity (the dissipation denoted by $id_{0}u$ is a small quantity). In general, for such equation no stable LB is possible because a balance between these four effects is hard to be achieved.

However, as indicated in the first section, different from the systems considered in Refs. Liu ; Minardi ; Eile ; Lah ; Bur ; Pec1 ; Pec2 ; Gur , stable, single LBs bounded in all spatial directions and in time can be realized in the present Rydberg-EIT system. The reasons are the following: (i) The nonlocal optical nonlinearity coming from the Rydberg-Rydberg interaction [described by the last term on the left hand side of Eq. (7)] has much fast response speed (with response time only in the order of $0.1\,\mu$s) Zhang , this is very different from the slow nonlocal optical nonlinearity in the systems studied in Refs. Lah ; Bur ; Pec1 ; Pec2 ; Gur where the response time of the nonlocal optical nonlinearity is in the range of 1 s or even longer; (ii) The nonlocal optical nonlinearity ($\chi_{p}^{(3)}\sim 10^{-8}{\rm m}^{2}{\rm V}^{-2}$) by the Rydberg-Rydberg interaction is much stronger and possesses a faster response speed than the local optical nonlinearity ($\chi_{p}^{(3)}\sim 10^{-11}{\rm m}^{2}{\rm V}^{-2}$) by the photon-atom interaction [described by the fourth term on the left hand side of Eq. (7), which has a response time in the order of $1\,\mu$s]. Based on these important properties (absent in the systems considered in Refs. Liu ; Minardi ; Eile ; Lah ; Bur ; Pec1 ; Pec2 ; Gur ), a single (3+1)D LB can form in the system via the following two-step mechanism for self-trapping: When a single (3+1)D probe pulse is incident into the system, it is firstly self-trapped in the two transverse (i.e. $x$, $y$) dimensions via the balance between the diffraction and the fast nonlocal optical nonlinearity; then it is further self-trapped in the longitudinal (i.e. $z$) direction by the balance between the dispersion and the slow local Kerr nonlinearity. Through such self-trapping processes, a stable single LB very different from those obtained in Refs. Liu ; Minardi ; Eile ; Lah ; Bur ; Pec1 ; Pec2 ; Gur can be realized by the synthetic nonlocal and local optical nonlinearities in the Rydberg atomic gas. Such LB can form in a very short distance and extremely low light power due to the giant optical nonlinearities and ultraslow propagation velocity of the probe pulse resulted from the Rydberg-based EIT effect.

Before presenting the LB solution, we make an estimation on the numerical values of the coefficients in Eq. (7). we choose $R_{0}=10$ ${\rm\mu m}$, $\tau_{0}=9\times 10^{-7}$ s, $U_{0}=3\times 10^{7}$ s*-1*, $\Delta_{2}=-2\pi\times 238\,{\rm MHz}$, $\Delta_{3}=-2\pi\times 0.32\,{\rm MHz}$, ${\cal N}_{a}=8.2\times 10^{10}$ cm*-3* and $\Omega_{c}=2\pi\times 16\,{\rm MHz}$. Then we obtain $R_{b}/R_{0}=0.6$, $g_{d}=0.134$, and $g_{1}=0.27$. We thus have the diffraction length $L_{\rm diff}=1.36$ mm and the dispersion length $L_{\rm D}(=L_{\rm LN})=10$ mm. As expected, $L_{\rm diff}\ll L_{\rm D}$, which supports the two-step mechanism for self-trapping described above (i.e., the diffraction plays its role earlier than the dispersion).

According to the above analysis, we can assume the LB solution takes the form

[TABLE]

where the parameter $w_{s}$ is the transverse beam width, $w_{t}$ is the pulse duration, $C_{s}$ is the wavefront curvature, and $C_{t}$ is the temporal chirp of the probe pulse. All the four parameters depend on variable $s$. In the ansatz (8), the factor $\exp\left[-(\xi^{2}+\eta^{2})/(2w_{s}^{2}(s))\right]$ is based on the balance between the diffraction and the nonlocal Kerr nonlinearity, and the factor ${\rm sech}[\sigma/w_{t}(s)]$ is based on the balance between the dispersion and the local Kerr nonlinearity.

We employ a variational method to solve Eq. (7) by taking the ansatz (8) to be a LB solution Rag . Through a Ritz optimization procedure, the LB energy $E$, defined by $E=\int\int\int|u|^{2}d\xi d\eta d\sigma=2\pi A^{2}w_{s}^{2}w_{t}$, is calculated as a function of the transverse beam width $w_{s}$, with the result shown in Fig. 3(a).

We observe that there are three branches for the LB solution, i.e., curves $C_{1}$, $C_{2}$, $C_{3}$. In the region where $\partial E/\partial w_{s}<0$ (curve $C_{2}$), the LB is quite stable. Yet, in regions $\partial E/\partial w_{s}>0$ (curves $C_{1}$ and $C_{3}$), the LB is unstable. This conclusion is verified by linearizing variational equations around the LB solution and examining their eigenvalues and eigenfunctions; it is also checked by using a numerical simulation on Eq. (7) directly. Fig. 3(b)-Fig. 3(d) show results for the pulse energy $E$, the beam width $w_{s}$, the pulse duration $w_{t}$ as functions of $z/(2L_{\rm diff})$, obtained respectively by choosing initial conditions from the curve $C_{1}$, $C_{2}$, and $C_{3}$ in Fig. 3(a), i.e. $(w_{s},w_{t},E)=(0.08,0.06,1.15)$ [Fig. 3(b)], $(w_{s},w_{t},E)=(0.66,0.41,0.44)$ [Fig. 3(c)], $(w_{s},w_{t},E)=(1.5,1.1,0.55)$ [Fig. 3(d)]. We see that only in the case of Fig. 3(c), the LB’s beam width $w_{s}$, pulse duration $w_{t}$, and energy $E$ keep almost unchanged, which means that the LB in the region of the curve $C_{2}$ is indeed stable during propagation. Note that the stable single LB solution obtained here is localized in all three spatial dimensions and also in time.

Moreover, Eq. (7) admits stable (3+1)D LBs carrying with orbital angular momentums (i.e., LVs). To demonstrate this, we take $u=u_{mp}(\xi,\eta,\sigma,\varphi)$ as a test solution, with

[TABLE]

where $L_{p}^{|m|}$ are the generalized Laguerre-Gauss (LG) polynomials, with $m$ and $p$ radial and azimuthal indices, respectively. The ansatz (9) (in the absence of the factor ${\rm sech}(\sigma/w_{t})$ with the normalization constant $C_{mp}=\sqrt{2^{|m|+1}p!/[\pi(p+|m|)!]}$ is called (LG)${}^{m}_{p}$ mode. Since $\hat{L}_{z}\,{\rm(LG)}^{m}_{p}=m\hbar\,{\rm(LG)}^{m}_{p}$, here $\hat{L}_{z}=-i\hbar\partial/\partial\varphi$, (LG)${}^{m}_{p}$ mode carry orbital angular momentum $m\hbar$ along the $z$ direction Andrews .

In Fig. 4 we illustrate the evolution of $|u|^{2}$. Panel (a) is for the LB [i.e., the fundamental mode (LG)${}_{0}^{0}$] with atomic density ${\cal N}_{a}=3\times 10^{10}$ cm*-3*; panel (b) is for the light vortex corresponding to mode (LG)${}_{0}^{1}$ [${\cal N}_{a}=4.95\times 10^{10}$ cm*-3*]. In the simulation, the beam radius $R_{0}$ for both the two modes is taken to be 10 $\mu$m. From Fig. 4(a) and Fig. 4(b), we see that for the lower-order modes (LG)${}_{0}^{0}$ and (LG)${}_{0}^{1}$, the pulses have nearly no deformation after propagating to $z=8L_{\rm diff}$ ($\approx 10.9$ mm).

However, pulses corresponding to the higher-order LG modes may not keep its shape in the nonlocal response region. To show this, we carried out the simulation on the evolution of the LV corresponding to (LG)${}_{1}^{2}$ by choosing $R_{0}=10\,\mu$m, $5\,\mu$m, and $2.5\,\mu$m, with the result shown in the panels (a), (b) and (c) of Fig. 5, respectively. We see that for a small beam radius [$R_{0}=2.5\,\mu$m; Fig. 5(c)], the vortex is relatively stable compared with that having a large beam radius [$R_{0}=10\,\mu$m; Fig. 5(a)]. This result means that if the degree of nonlocality in the system increases (i.e., $R_{0}$ is reduced for fixed blockade radius $R_{b}$), the lifetime of the vortex pulse may be increased significantly. When degree of nonlocality becomes very large, the system enters into a strongly nonlocal response region, in this case all the LVs are stable during propagation; see the discussion given below.

IV.3 Strongly nonlocal response region

We lastly consider the situation when the range of the Rydberg-Rydberg interaction is much larger than that of the beam radius of the probe pulse, i.e. $R_{b}/R_{0}\gg 1$, which corresponds to the case shown in Fig. 2(c), where Re($G_{2}$), Im($G_{2}$), and $|U/U_{0}|^{2}$ as functions of $r_{\bot}^{\prime}/R_{b}$ is shown by taking $R_{0}=1\,{\rm\mu m}$ (and hence $R_{b}/R_{0}=5.8$). We see that, compared with $|U/U_{0}|^{2}$, the response function $G_{2}$ is very flat near $r_{\bot}^{\prime}=r_{\bot}$, which means that $G_{2}({\bf r}_{\bot}^{\prime}-{\bf r}_{\bot})\approx G_{2}(r_{\bot})\approx G_{2}(0)+G_{2}^{\prime\prime}(0)r_{\bot}^{2}/2$ snyder ; Krolikowski , plotted by the purple dashed line in Fig. 2(c), agreeing well with the exact one near $r_{\bot}=0$ [i.e., Re($G_{2}$) (blue solid line); Im($G_{2}$) is very small]. In this case, Eq. (3) can be reduced into the dimensionless form

[TABLE]

where $U=U_{0}u\,{\rm exp}(ig_{3})$, $g_{3}=2L_{\rm diff}P_{0}G_{2}(0)$, and $g_{4}=-L_{\rm diff}R_{0}^{2}P_{0}G_{2}^{(2)}(0)$, with $P_{0}=\int d^{2}{\bf r}_{\bot}|U|^{2}$ the power of the probe pulse (approximately a constant). The definitions of $s$, $\eta$, and $\xi$ are the same as those in Sec. IV.2. Thus, in the strongly nonlocal response region the nonlocal Kerr nonlinearity contributed by the Rydberg-Rydberg interaction reduces into a parabolic “external potential”. The physical reason for this reduction is that all the photons in the probe pulse experience an almost alike (effective) potential due to the very narrow probe beam radius and the very wide distribution of the potential.

Before presenting LB and LV solutions. we first give an example of numerical values of the coefficients in Eq. (10), we choose $\Delta_{3}=-2\pi\times 48\,{\rm kHz}$, $\Delta_{2}=-2\pi\times 160\,{\rm MHz}$, $R_{0}=1.86$ ${\rm\mu m}$, $\tau_{0}=2.6\times 10^{-7}$ s, ${\cal N}_{a}=5.8\times 10^{12}$ cm*-3*, $U_{0}=1.5\times 10^{7}$ s*-1*, and $C_{6}=\simeq 2\pi\times 167\,{\rm THz\,\mu m}^{6}$ (for the principle quantum number $n=120$). With these parameters, we obtain diffraction length $L_{\rm diff}=0.047$ ${\rm mm}$, and hence the blockade radius $R_{b}=19$ ${\rm\mu m}$ ($R_{b}/R_{0}=10.2$), $g_{d}=1.0-0.07i$, $g_{1}=2.0-0.014i$, $d_{0}=-0.03$, $g_{4}=2+0.1i$. Thus, the system is within the strongly nonlocal region with very small dissipation.

Similarly, a variational method is employed to solve Eq. (10) by using the ansatz (9) as test LB and LV solutions; their stability is investigated through a linear stability analysis. Then a numerical simulation starting directly from Eq. (7) working in this strong nonlocal response region is carried out by taking the solutions obtained via the variational method as an initial condition, together with a random disturbance to it. Specifically, we take $U(z=0,x,y,t)=U_{0}\,(1+\varepsilon f_{R})\,u_{mp}(z=0,x,y,t){\rm exp}(ig_{3})$, with $\varepsilon$ being a typical amplitude of the perturbation, and $f_{R}$ being a random variable uniformly distributed in the interval [0,1]. We find that the system allows indeed stable LBs and LVs. Shown in Fig. 6(a)

and Fig. 6(b) are time evolution of the (3+1)D LB and LV (corresponding to the fundamental mode (LG)${}_{0}^{0}$ and the higher-order mode (LG)${}_{1}^{2}$) by taking $\varepsilon=0.1$, where isosurface plots are illustrated of these pulses when propagating to distances $s=z/(2L_{\rm diff})=0,1,2,3,4$, respectively. We observe that these (3+1)D nonlinear optical pulses relax to self-cleaned forms quite close to the unperturbed ones, and their shapes undergo no apparent change during propagation.

Based on the solution (9) and using the parameters given above, we get the propagating velocity of the LBs and LVs

[TABLE]

much slower than the light speed in vacuum. The maximum average power density $P_{\rm max}$ for generating such LBs and LVs can be obtained by using Poynting’s vector huang ; Newell , which is estimated to be

[TABLE]

which corresponds to the maximum average peak intensity $\bar{I}_{\rm max}\simeq 3.6\times 10^{-4}$ W cm*-2*. Consequently, the (3+1)D LBs and LVs obtained in the present Rydberg-EIT system have ultraslow propagating velocity and extremely low generation power, which are very different from other schemes.

V Storage and retrieval of (3+1)D weak light bullets and vortices

One of main advantages of EIT is the possibility for realizing an active manipulation of optical pulses by tuning system parameters. Especially, optical pulses can be stored and retrieved through the switching off and on of the control field. In recent years, a large amount of studies on light memory by using EIT have been carried out FIM ; Chen ; Novikova , including ones performed in Rydberg atomic systems Maxwell2013 ; LiL . However, it is generally difficult to realize the memory of high-dimensional nonlinear optical pulses via conventional EIT because of the catastrophic collapse during propagation. Nevertheless, below we show that the storage and retrieval of the LBs and LVs with high efficiency and high fidelity are possible in the present Rydberg-EIT system.

To this end, we investigate the evolution of the LBs and LVs described above through solving the MB Eqs. (2) and (1) by using a control field that is switched on and off adiabatically, which can be described by the switching function with the form, $\Omega_{c}(t)=\Omega_{c0}\{1-1/2\tanh[(t-T_{\rm off})/T_{s}]+1/2\tanh[t-T_{\rm on}/T_{s}]\}$, where $T_{\rm{off}}$ and $T_{\rm{on}}$ are times at which the control field is of switched off and on. The duration of the switching time is $T_{s}$ and the storage time of the probe pulse is approximately given by $T_{\textrm{on}}-T_{\textrm{off}}$.

V.1 Memory in the nonlocal response region

We first study the storage and retrieval of the (3+1)D LB in the nonlocal response regime. The spatial waveshape of the input probe pulse is taken as a fundamental LG mode [i.e., (LG)${}_{0}^{0}$ mode], and the temporal profile is assumed as a hyperbolic secant one. Fig. 1(c) shows the numerical result of the probe-pulse intensity $|U/U_{0}|^{2}$ during the process of storage and retrieval. The pulse waveshapes for $z=0$ (before the storage), $z=4L_{\rm diff}$ (at the beginning of the storage), and $z=8L_{\rm diff}$ (after the storage), with $L_{\rm diff}=1.36$ mm, are plotted. We see that when the control field $\Omega_{c}$ is switched on, the (3+1)D LB is created; by switching off the control field, the LB is stored in the atomic medium; then the LB is retrieved when $\Omega_{c}$ is switched on again. We see that the retrieved LB has nearly the same waveshape as that before the storage.

It is possible to store LVs in the nonlocal response regime in the same way. In Fig. 7 we show the storage and retrieval of a (3+1)D LV of (LG)${}_{0}^{1}$ mode. The black dashed line in the figure shows the process of the switching-on, switching-off, and re-switching-on of the control field $|\Omega_{c}\tau_{0}|$. The curves 1, 2 and 3 are temporal profiles of the probe pulse $|\Omega_{p}\tau_{0}|$ when the LV propagates at $z=0$ (before the storage), $z=4L_{\rm diff}$ (at the beginning of the storage), and $8L_{\rm diff}$ (after the storage), respectively. Corresponding isosurface plots with $|\Omega_{p}\tau_{0}|=0.1$ illustrate the storage of retrieval of the LV.

When used as optical memories, we need to know the efficiency $\eta$, which can be described by the energy ratio between the retrieved pulse and the input pulse Novikova3 , i.e., $\eta=I_{\rm out}/I_{\rm in}$, where $I_{\rm out}=\int_{T_{\rm on}}^{\infty}dt\iint dxdy\,|\Omega_{p}^{\rm out}(x,y,t)|^{2}$, $I_{\rm in}=\int^{T_{\rm off}}_{-\infty}dt\iint dxdy\,|\Omega_{p}^{\rm in}(x,y,t)|^{2}$, $\Omega_{p}^{\rm in}(x,y,t)=\Omega_{p}(x,y,z,t)|_{z=0}$ and $\Omega_{p}^{\rm out}(x,y,t)=\Omega_{p}(x,y,z,t)|_{z=L_{z}}$, with $L_{z}$ the medium length in the propagation direction. Based on this, the memory efficiency of the LB (LV with (LG)${}_{0}^{1}$ mode) is given by $\eta=93.12\%$ ($\eta=92.38\%$) for $L_{z}=8L_{\rm diff}=10.9\,{\rm mm}$.

The quality for the preservation of pulse waveshape in the memory can be characterized by the overlap integral Novikova3 , $J^{2}=I_{\rm inter}^{2}/(I_{\rm in}I_{\rm out})$, where $I_{\rm inter}=|\int^{T_{\rm off}}_{-\infty}dt\iint dxdy\,\Omega_{p}^{\rm out}(x,y,t+\Delta T)\Omega_{p}^{\rm in}(x,y,t)|$, $\Delta T$ is the time interval between the peaks of the input and the retrieved probe pulses. In the case shown in Fig. 7, $\Delta T=28\tau_{0}$, and hence we obtain $J^{2}=97.26\%$ ($J^{2}=96.35\%$) for the LB (LV). Combined the above two parameters, one can determine the fidelity of the memory, defined by $F=\eta J^{2}$. From these results, we obtain $F=90.57\%$ ($F=89.01\%$) for the LB (LV). Thus the storage and retrieval of the LB and LV in the Rydberg atomic gas have not only high efficiency but also high fidelity.

For comparison, Fig. 8

shows the result of the memory for the LB when the Rydberg-Rydberg interaction is absent. We see that the LB suffers a significant deformation after the storage, with the fidelity of memory only $6.3\%$ (for the memory of LVs without the Rydberg-Rydberg interaction, the fidelity is even lower), which cannot be applied for light information processing because the information is lost during the storage.

V.2 Memory in the strong nonlocal response region

We now consider the storage and retrieval of the (3+1)D LBs and LVs in the strong nonlocal response region. Shown in Fig. 9(a)

[Fig. 9(b)] is the numerical result of the storage and retrieval of a (3+1)D LB [LV with the (LG)${}_{1}^{2}$ mode]. In the figure, the black dashed line is the time sequence of the switching-on and switching-off of the control field $|\Omega_{c}\tau_{0}|$. Curves 1, 2 and 3 are temporal profiles of the probe pulse $|\Omega_{p}\tau_{0}|$ respectively at $z=0$ (before the storage), $z=6L_{\rm diff}$ (at the beginning of the storage), and $12L_{\rm diff}$ (after the storage), with $L_{\rm diff}=0.047\,{\rm mm}$; the corresponding isosurface plots for $|\Omega_{p}\tau_{0}|=0.1$ are also illustrated. Here system parameters used in the calculation are the same as those given in Sec. IV.3 but with $R_{0}=1.86\,{\rm\mu m}$, $R_{b}=19\,{\rm\mu m}$, and $\tau_{0}=2.7\times 10^{-7}$ s. From the figure we see that the LB (LV) after the storage keeps nearly the same waveshape as that before the storage. The efficiency and fidelity of the LB memory (LV memory) in this strong nonlocal region can reach to $\eta=93.97\%$ and $F=91.86\%$ ($\eta=93.14\%$ and $F=89.11\%$). Even higher values can be reached by tuning systemic parameters (i.e., transverse beam radius $R_{0}$ or principle quantum $n$) to increase the degree of nonlocality.

The calculation on the storage and retrieval of other (higher-order) LVs in the strong nonlocal response region is also carried out (not illustrated here). The result shows that the memory of the higher-order LVs can also be realized and have also high efficiency and fidelity. All these results indicate that the Rydberg medium supports high-quality memory for various LG modes. Since LG modes carry definite orbital angular momentum and hence more information, the acquirement of high-quality memory by using the strong, nonlocal Rydberg-Rydberg interaction provides the possibility to realize a high-quality multi-dimensional light memory.

VI Summary

In this work, we have carried out a detailed investigation on the formation, propagation, and storage of ultraslow weak-light bullets and vortices via Rydberg-EIT in a cold atomic gas. By using an approach beyond mean-field theory, we have shown that the system may acquire two types of optical nonlinearities, i.e., the giant fast-responding nonlocal optical nonlinearity (contributed by the Rydberg-Rydberg interaction) and the relatively weak, slow-responding local Kerr nonlinearity (contributed by the photon-atom interaction), both of them are enabled by the Rydberg-EIT. We have derived a (3+1)D nonlinear envelope equation governing the spatiotemporal evolution of the high-dimensional probe pulse and present various (3+1)D light bullet and vortex solutions. The light bullets and vortices obtained have very slow propagating velocity and extremely low generation power, and can be stabilized by the interplay between the synergetic local and nonlocal optical nonlinearities. These (3+1)D nonlinear optical pulses can be dynamically manipulated based on the active character of the system. Our study expands the breadth in the study of higher dimensional, nonlocal nonlinear optics with Rydberg atomic gases. Recent experiments have reported the proof-of-principle demonstration of storing light modes with orbital angular momentumVeissier ; Nicolas ; Parigi ; Oliveira ; ding . The highly efficient method proposed in this work has potential applications in topological quantum information processing.

Acknowledgments

This work was supported by the NSF-China under Grants No. 11174080 and No. 11474099, by the China Postdoctoral Science Foundation under Grant No. 2017M620140, and by the Shanghai Sailing Program under Grant No. 18YF1407100. W.L. acknowledges support from the UKIERI-UGC Thematic Partnership No. IND/CONT/G/16-17/73 and EPSRC Grant No. EP/M014266/1 and EP/R04340X/1.

Appendix A EXPRESSIONS OF THE SECOND-ORDER SOLUTION

Expressions of the second-order solution are given by

[TABLE]

where $\rho_{\alpha 1}^{(2)}=a_{\alpha 1}^{(2)}\partial F/\partial t_{1}{\rm exp}(i\theta)$, $\rho_{32}^{(2)}=a_{32}^{(2)}|F|^{2}{\rm exp}(-2\bar{\alpha}z_{2})$, $\rho_{\beta\beta}^{(2)}=a_{\beta\beta}^{(2)}|F|^{2}{\rm exp}(-2\bar{\alpha}z_{2})$($\alpha=2,3;\beta=1,2,3$) with $\theta=K(\omega)z_{0}-\omega t_{0}$, $\bar{\alpha}=\epsilon^{-2}\alpha=\epsilon^{-2}$Im$[K(\omega)]$, $M_{1}=1/d_{32}-1/d_{32}^{\ast}$, $M_{2}=(\omega+d_{31}^{\ast})/D(\omega)^{\ast}-(\omega+d_{31})/D(\omega)$, and $M_{3}=1/[D(\omega)^{\ast}d_{32}^{\ast}]-1/[D(\omega)d_{32}]$.

Appendix B Expansion equations of the two-body correlators

Equations of motion for the two-body correlators at second-order and third-order approximations are given in the following.

(i) Second-order approximation: we have

[TABLE]

[TABLE]

The second-order solution for two-body correlators is given by $\rho_{\alpha 1,\beta 1}^{(2)}=a_{\alpha 1,\beta 1}^{(2)}F^{2}e^{2i\theta}$, $\rho_{\alpha 1,1\beta}^{(2)}=a_{\alpha 1,1\beta}^{(2)}|F|^{2}e^{-2\bar{\alpha}z_{2}}$ $(\alpha,\beta=2,3)$. Note that due to the weak probe field the vdW potential $\hbar V$ has no contribution to $\rho_{\alpha 1,1\beta}^{(2)}$, thus $\rho_{\alpha 1,1\beta}^{(2)}=\rho_{\alpha 1}^{(1)}\rho_{1\beta}^{(1)}$. Explicit expression of $\rho_{\alpha 1,\beta 1}^{(2)}$ can be directly obtained by solving Eq. (B).

(ii) Third-order approximation. At this order, we have

[TABLE]

where $M_{31}=\omega+i\Gamma_{12}+d_{21}$, $M_{32}=\omega+i\Gamma_{12}+d_{31}$, $M_{33}=\omega+i\Gamma_{23}+d_{21}$, $M_{34}=\omega+d_{31}+i\Gamma_{23}-V$, $M_{35}=\omega+d_{32}+d_{21}$, $M_{36}=\omega+d_{23}+d_{21}$, $M_{37}=\omega+d_{32}+d_{31}-V$ and $M_{38}=\omega+d_{23}+d_{31}$. Form these equations we obtain the third order solution

[TABLE]

with

[TABLE]

Here $P_{n}$ and $Q_{n}~{}(n=0,1,2,3)$ are functions of the spontaneous emission decay rate $\gamma_{\mu\nu}$, detunings $\Delta_{\mu}$, and half Rabi frequency $\Omega_{c}$.