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

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.


I. INTRODUCTION
High-dimensional spatiotemporal optical solitons, alias light bullets (LBs) [1], 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 [2] because of their rich nonlinear physics and technological applications [3][4][5][6][7][8]. Experimental signatures of LBs in different types of optical media, such as (2+1)D LBs in LiIO 3 crystals [9] and quasi-(3+1)D LBs in waveguide arrays [10,11], have been reported. Recently, trains of (3+1)D dark solitons have been observed in a photorefractive material by Lahav et al. [12]. 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 [13][14][15] and lead glass [16], in which the nonlocal nonlinearity (resulted from the reorientation [13][14][15] or thermal motion [16] 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 (∼ c) of optical pulses, stable LBs can only be realized in the form of pulse-train solitons (not single LBs) [12]. A different approach is based on electromagnetically induced transparency (EIT) in atomic gases [17], in which optical absorption can be largely suppressed due to quantum interference. Together with Kerr nonlinearities [18,19] induced by resonant laser fields, it has been shown that stable (1+1)D temporal [20][21][22] and spatial [23,24] 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 [25] and experimental [26] studies revealed that strong and longrange optical nonlinearities can be built with Rydberg atomic gases [27][28][29], which are in electronically highlying states with large principal quantum number n [30]. 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 [27-29, 31, 32]. Importantly, the Rydberg-Rydberg interaction can be mapped to a nonlocal optical nonlinearity through EIT, which is strong even at the single photon level [27][28][29]. Long lifetimes (in the order of tens of microsecond) guarantee that the induced photon-photon interaction is largely coherent during light propagation [33][34][35]. This provides a new platform to study quantum nonlinear optics [28] and develop new photon devices, such as single-photon switches and transistors [36][37][38], quantum memories and phase gates [39][40][41][42][43][44].
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 [45], which is complemented by the latter [46,47], 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 meanfield theory [48][49][50][51], 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) ] [52][53][54][55][56][57]. 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 laddertype level configuration, shown in Fig. 1(a). A weak, pulsed probe field of angular frequency ω p (half Rabi frequency Ω p ) couples the ground state |1 and inter- , and Rydberg state |3 are respectively driven by a pulsed probe field (with pulse duration τ 0 ) and a strong control field. State |2 has a large spontaneous decay rate Γ 12 ∼ MHz. The weak decay Γ 23 ∼ kHz from |3 to |2 is also taken into account. The van der Waals interaction V (r − r ′ ) between the two atoms in Rydberg states respectively located at r and r ′ shifts the Rydberg state energy. (b) Geometry of the system. The probe and control laser fields counter-propagate in the Rydberg gas. The involution of the weak probe pulse in the Rydberg gas is considered, while depletion of the strong control field is neglected. (c) Storage and retrieval of a (3+1)D light bullet, illustrated by an isosurface plot of the light intensity of the light bullet before storage (z = 0), at the beginning of the storage (z = 5.4 mm) and after the storage (z = 10.8 mm); see text for details. mediate state |2 , and a strong, continuous-wave control field of angular frequency ω c (half Rabi frequency Ω c ) couples state |2 and a Rydberg state |3 . The probe field has a pulse length τ 0 at the entrance of the medium. The electric-field vector of the system is E(r, t) = l=p,c e l E l exp[i(k l · r − ω l t)] + c.c., where e l (k l ) is the unit polarization vector (wavevector) of the electric-field component with envelope E l (l = p, c). A possible experimental arrangement of beam geometry is shown in Fig. 1 The system works at a ultracold temperature and the probe and control fields counter-propagate [i.e., k p = (0, 0, k p ), k c = (0, 0, −k c )], so that the center-ofmass motion of the atoms and dephasing due to the atomic collisions are negligible. Under the electricdipole approximation, the system Hamiltonian isĤ H = N a d 3 rĤ H (r, t) witĥ where N a is atomic density, ω α is the eignenergy of atomic state |α , Ω p = (e p · p 21 )E p / and Ω c = (e c · p 32 )E c / are respectively the half Rabi frequencies of the probe and control fields (with p αβ the electric dipole matrix element associated with the transition from |β to |α ),Ŝ αβ are transition operators (α, β = 1, 2, 3) satisfying the commutation relation Ŝ αβ (r, t),Ŝ µν (r ′ , t) = δ ανŜµβ (r, t) − δ µβŜαν (r ′ , t) δ rr ′ . In the last term on the right hand side (RHS), V (r − r ′ ) = − C 6 /|r − r ′ | 6 is the van der Waals (vdW) interaction potential with dispersion coefficient C 6 between the Rydberg atoms located at the positions r and r ′ , respectively. To be concrete, we will consider strontium atoms ( 88 Sr) 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 = |5s 2 1 S 0 , |2 = |5s5p 1 P 1 , |3 = |5sns 1 S 0 , with n the principle quantum number [58]. The spontaneous emission rates of the atoms are given by Γ 12 = 2π × 16 MHz, Γ 23 = 2π × 16.7 kHz, so one has γ 21 = Γ 12 /2, γ 31 = Γ 23 /2, γ 32 = (Γ 12 + Γ 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 [17], is the coherence between level |1 and |2 . Its dynamics is determined by the Bloch equations, where ρ αβ (r, t) ≡ Ŝ αβ (r, t) are the one-body correlators (one-body density matrix elements) [59], d αβ = ∆ α − ∆ β +iγ αβ (∆ 1 = 0; α, β = 1, 2, 3; α = β), ∆ 2 = ω p −(ω 2 − ω 1 ) and ∆ 3 = ω p +ω c −(ω 3 −ω 1 ) are respectively the onephoton and two-photon detunings, γ αβ = (Γ α + Γ β )/2 + γ col αβ with Γ β = α<β Γ αβ . Here Γ αβ denotes respectively the spontaneous emission decay rate from the state |β to the state |α , and γ col αβ represents the dephasing rate reflecting the loss of phase coherence between |α and |β due to, e.g., atomic motion and the interaction between the atoms in the ground and Rydberg states.
From Eq.

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 [21,22,60] widely applied for solving nonlinear wave problems [60]. Because our calculation is exact to third order (i.e. up to Ω 3 p ), higher-order nbody correlators (n ≥ 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 and make the asymptotic expansion [21,22]: where ǫ a parameter characterizing the order of Ω p . To obtain a divergencefree expansion, all quantities on the right hand side of the expansion are considered as functions of the multiscale variables z l = ǫ l z (l = 0, 1, 2), (x 1 , y 1 ) = ǫ (x, y), and t l = ǫ 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 ǫ 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 Ω [61]. 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(ω) is the linear dispersion relation, given by K(ω) = ω/c + κ 12 (ω + d 31 )/D(ω). At the second order, a solvability (i.e. divergence-free) condition requires ∂F/∂z 1 + (1/V g )∂F/∂t 1 = 0, with V g = (∂K/∂ω) −1 the group velocity of the envelope F . The explicit expression of the second-order approximate solution for the one-body correlators, i.e. ρ (2) αβ , is presented in Appendix A. Expansion equations for the twobody 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 : 33 )]/D(ω) is the Kerr coefficient describing the local nonlinear optical response contributed from the photon-atom interaction (which vanishes if the twophoton detuning ∆ 3 is taken to be zero). Explicit expressions of a 32 , and a (2) 33 are presented in Appendix A.
33,31 (r ′ −r)V (r ′ − r)dz ′ (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 where 33,31 (r ′ − 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 [49].
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 ≈ 0 when ∆ 3 = 0, it is necessary to take a non-zero ∆ 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. as functions of r ′ ⊥ /R b . We show the local regime (a) with R 0 = 300 µm, nonlocal regime (b) with R 0 = 10 µm, and strong nonlocal regime (c) with R 0 = 1 µm. In all situations, real parts (blue solid line) dominates over the imaginary part (orange dashed line). For a better visualization, G 2 has been amplified 10 8 times. We also show the intensity profile of the probe field |U These parameters guarantee that the system is in the dispersive nonlinearity regime (i.e. |∆ 2 | ≫ Γ 12 ).

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 [49]; 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 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., |∆ 2 | ≫ Γ 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 [29,49]: , for which different types of LB solutions can be obtained.

A. 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 ≪ 1. Fig. 2(a) shows (imaginary part of G 2 ), and |U/U 0 | 2 (the intensity profile of the probe field) as functions of r ′ ⊥ /R b by taking R 0 = 300µ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 Here . That is to say, for large probefield 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 Under the EIT condition (|Ω c | 2 > γ 21 γ 31 ) and a large one-photon detuning ∆ 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 and L A ≡ 1/α 0 the typical diffraction length, dispersion length, and absorption length, respectively. Note that we have taken L D = L NL [with L NL ≡ 1/(U 2 0W 1 + U 2 0W 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 ≡ (1/τ 0 )[|K 2 /(W 1 +W 2 )|] 1/2 . The tilde above corresponding quantities means taking their real parts. Due to the EIT effect and large ∆ 2 , we have d 0 ≪ 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 diff (≫ L D ) due to large R 0 , thus we have g diff ≪ 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 if K 2 < 0 (i.e. s d = −1), withK 0 =K| ω=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) × 10 −13 cm −1 s 2 , W 1 = (6.86 + 0.062i) × 10 −18 cm −1 s 2 , W 2 = (1.79 + 0.0076i) × 10 −14 cm −1 s 2 . We see that, as expected, the imaginary parts of these coefficients are indeed much smaller than their real parts. By taking τ 0 = 1.2 × 10 −7 s, we obtain U 0 = 2 × 10 7 s −1 , L D = L 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µm), one has L diff = 1226 mm and hence g 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Ṽ g = 9.8 × 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 χ , is estimated to have the value (3.03 + 0.0129i) × 10 −8 m 2 V −2 , which is 11 orders of magnitude higher than that obtained in conventional nonlinear optical media [62]. 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 [21,63], the maximum average power density to generate such ultraslow optical soliton is estimated to beP max = 1.2 µW. Thus, very low input power is needed for generating the optical soliton in the system.

B. Nonlocal response region
We now turn to consider nonlocal response region, corresponding to R b /R 0 ∼ 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 through new scalings . 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. [9][10][11][12][13][14][15][16], 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 µs) [45], this is very different from the slow nonlocal optical nonlinearity in the systems studied in Refs. [12][13][14][15][16] where the response time of the nonlocal optical nonlinearity is in the range of 1 s or even longer; (ii) The nonlocal optical nonlinearity (χ (3) p ∼ 10 −8 m 2 V −2 ) by the Rydberg-Rydberg interaction is much stronger and possesses a faster response speed than the local optical nonlinearity (χ (3) p ∼ 10 −11 m 2 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 µs]. Based on these important properties (absent in the systems considered in Refs. [9][10][11][12][13][14][15][16]), 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. [9][10][11][12][13][14][15][16] 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.
According to the above analysis, we can assume the LB solution takes the form 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 −(ξ 2 + η 2 )/(2w 2 s (s)) is based on the balance between the diffraction and the nonlocal Kerr nonlinearity, and the factor sech[σ/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 [64]. Through a Ritz optimization procedure, the LB energy E, defined by E = |u| 2 dξdηdσ = 2πA 2 w 2 s 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 ∂E/∂w s < 0 (curve C 2 ), the LB is quite stable. Yet, in regions ∂E/∂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 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 (ξ, η, σ, ϕ) as a test solution, with where L |m| p 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 sech(σ/w t ) with the normalization constant LG) m p mode carry orbital angular momentum m along the z direction [52].
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 N a = 3 × 10 10 cm −3 ; panel (b) is for the light vortex corresponding to mode (LG) 1 0 [N a = 4.95 × 10 10 cm −3 ]. In the simulation, the beam radius R 0 for both the two modes is taken to be 10 µm. From Fig. 4(a) and Fig. 4(b), we see that for the lowerorder modes (LG) 0 0 and (LG) 1 0 , the pulses have nearly no deformation after propagating to z = 8L diff (≈ 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) 2 1 by choosing R 0 = 10 µm, 5 µm, and 2.5 µ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 µm; Fig. 5(c)], the vortex is relatively stable compared with that having a large beam radius [R 0 = 10 µ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.  2 1 in the nonlocal response region, as a function of x/R b and y/R b when propagating to the distance at z/(2L diff ) = 0, 1, 2, 3, and 4, respectively. Parameters are R0 = 10 µm, Na = 9.9 × 10 10 cm −3 (a), R0 = 5 µm, Na = 1.46 × 10 11 cm −3 (b), and R0 = 2.5 µm, Na = 5.3 × 10 11 cm −3 (c). Other system parameters are the same as those used in Fig. 4.

C. 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 ≫ 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 ′ ⊥ /R b is shown by taking R 0 = 1 µ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 ′ ⊥ = r ⊥ , which means that [65,66], plotted by the purple dashed line in Fig. 2(c), agreeing well with the exact one near r ⊥ = 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 where U = U 0 u exp(ig 3 ), g 3 = 2L diff P 0 G 2 (0), and 2 (0), with P 0 = d 2 r ⊥ |U | 2 the power of the probe pulse (approximately a constant). The definitions of s, η, and ξ are the same as those in Sec. IV B. 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.
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 + εf R ) u mp (z = 0, x, y, t)exp(ig 3 ), with ε 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) 2 1 ) by taking ε = 0.1, where isosurface plots are illustrated of these pulses when propagating to distances s = z/(2L 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 much slower than the light speed in vacuum. The maximum average power density P max for generating such LBs and LVs can be obtained by using Poynting's vector [21,63], which is estimated to bē which corresponds to the maximum average peak inten-sityĪ max ≃ 3.6 × 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 [17,22,67], including ones performed in Rydberg atomic systems [39,41]. 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, where T off and T 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 on − T off .

A. 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 diff (at the beginning of the storage), and z = 8L diff (after the storage), with L diff = 1.36 mm, are plotted. We see that when the control field Ω 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 Ω 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) 1 0 mode. The black dashed line in the figure shows the process of the switching-on, switching-off, and re-switching-on of the . The black dashed line shows the switchingon, switching-off, and re-switching-on of the control field |Ωcτ0|. The curves 1, 2 and 3 are temporal profiles of the probe pulse |Ωpτ0| respectively at z = 0 (before the storage), z = 4L diff (at the beginning of the storage), and 8L diff (after the storage), with L diff = 1.36 mm; the corresponding isosurface plots for |Ωpτ0| = 0.1 are also shown. In the calculation, we have set R0 = 10 µm, C6 ≃ 2π × 81.6 GHz µm 6 , Ωpτ0 = 27 s −1 , ∆2τ0 = −1340, ∆3τ0 = −1.8, and Ωc0τ0 = 90 with τ0 = 9 × 10 −7 s. The control parameters are Ts = 0.2τ0, T off = 5τ0, Ton = 15τ0.
control field |Ω c τ 0 |. The curves 1, 2 and 3 are temporal profiles of the probe pulse |Ω p τ 0 | when the LV propagates at z = 0 (before the storage), z = 4L diff (at the beginning of the storage), and 8L diff (after the storage), respectively. Corresponding isosurface plots with |Ω p τ 0 | = 0.1 illustrate the storage of retrieval of the LV.
When used as optical memories, we need to know the efficiency η, which can be described by the energy ratio between the retrieved pulse and the input pulse [68], i.e., η = I out /I in , where I out = ∞ Ton dt dxdy |Ω out p (x, y, t)| 2 , dxdy |Ω in p (x, y, t)| 2 , Ω in p (x, y, t) = Ω p (x, y, z, t)| z=0 and Ω out p (x, y, t) = Ω p (x, y, z, t)| z=Lz , with L z the medium length in the propagation direction. Based on this, the memory efficiency of the LB (LV with (LG) 1 0 mode) is given by η = 93.12% (η = 92.38%) for L z = 8L diff = 10.9 mm.
The quality for the preservation of pulse waveshape in the memory can be characterized by the overlap integral [68], dxdy Ω out p (x, y, t + ∆T )Ω in p (x, y, t)|, ∆T is the time interval between the peaks of the input and the retrieved probe pulses. In the case shown in Fig. 7, ∆T = 28τ 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 = η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.

B. 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) 2 1 mode]. In the figure, the black dashed line is the time sequence of the switching-on and switchingoff of the control field |Ω c τ 0 |. Curves 1, 2 and 3 are temporal profiles of the probe pulse |Ω p τ 0 | respectively at z = 0 (before the storage), z = 6L diff (at the beginning of the storage), and 12L diff (after the storage), with L diff = 0.047 mm; the corresponding isosurface plots for |Ω p τ 0 | = 0.1 are also illustrated. Here system parameters used in the calculation are the same as those given in Sec. IV C but with R 0 = 1.86 µm, R b = 19 µm, and τ 0 = 2.7×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 η = 93.97% and F = 91.86% (η = 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 highquality 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 multidimensional 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 momentum [53][54][55][56][57]. The highly efficient method proposed in this work has potential applications in topological quantum information processing.