## Abstract

A scheme is proposed to generate stable light bullets (LBs) in a cold Rydberg atomic system with a parity-time (PT) symmetric potential, by utilizing electromagnetically induced transparency (EIT). Using an incoherent population pumping between two low-lying levels and spatial modulations of control and auxiliary laser fields, we obtain a two-dimensional (2D) periodic optical potential with PT symmetry. Based on PT symmetry potential and the long-range Rydberg-Rydberg atomic interaction, the system may support slow LBs with low light intensity. Further, it is found that the local and non-local nonlinear coefficients and PT-symmetric potential can be tuned and used to manipulate the behavior of LBs.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Spatiotemporal localization of waves is a topic that attracts broad interest in physics, ranging over diverse fields such as plasma physics, Bose-Einstein condensates (BECs), and nonlinear optics. In nonlinear optics, owing to their curious physics and technological applications, LBs attracted intensive theoretical [1–3] and experimental interest in recent years [4–6]. However, LBs generally suffer severe instabilities, typically propagating for just a few diffraction lengths. The overwhelming demand is the identification of optical media with tunable dispersion, diffraction, and nonlinearities that can support the generation of stable LBs. It is encouraging that a number of setups have been proposed where the stabilization of LBs is possible, such as in materials with saturable, competing or nonlocal nonlinearities [7,8], special nonlinear interactions [3,9], optical tandem systems [10,11], waveguide arrays [12,13], PT-symmetric optical lattices [14,15], and EITs in cold Rydberg atomic gases [16].

The introduction of PT-symmetric optical lattices and EIT in a cold Rydberg gas, have significantly promoted the stability of LBs in nonlinear optical systems [14–16]. Kartashov *et al*. first displayed the 3D fundamental and vortex-carrying solitons in PT-symmetric periodic lattices [14]. The imaginary part of lattice affect solitons stability. The stability of fundamental solitons (FSs) can be extended almost to the PT symmetry breaking point with the complex spectrum. Vortex solitons (VSs) display asymmetric profiles, but still found to be stable within narrow parameter regions. On the other hand, theoretical [17] and experimental [18] studies revealed that strong and long-range optical nonlinearities can be built with Rydberg atoms [19]. F. Maucher *et al*. proposed a scheme for the creation of stable 3D bright solitons in Bose-Einstein condensates. Off-resonant dressing to Rydberg nD states provides nonlocal attractive interactions, leading to self-trapping of mesoscopic atomic clouds by a collective excitation of a Rydberg atom pair [7]. Huang *et al*. realized stable single LBs and their active control in cold Rydberg gases [16]. Utilizing EIT, they succeeded in mapping the strong atom-atom interaction in Rydberg states to light fields, which resulted in a giant nonlocal Kerr nonlinearity and the formation of LBs with low power and slow propagation velocity.

Importantly, the Rydberg-Rydberg atomic interaction can be mapped to a nonlocal optical nonlinearity through EIT [17]. Rydberg atomic interactions find applications in quantum information, precision measurement, quantum simulation, and many-body physics [20–23]. It will be interesting if one can obtain realistic optical systems that not only possess PT symmetry, but also support nonlocal weak-light LBs. In this work, we derive a nonlinear (3 + 1)D light propagation equation in a cold Rydberg atomic system with PT symmetric potential, taking into account many-body correlations. Based on the strong interaction between Rydberg atoms through the EIT [24–27], we reveal that the degree of nonlocality of the Kerr nonlinearity and the real part and the imaginary part of the PT potential can be used to control the behavior of LBs. Further, we obtain stable different spatiotemporal distributions of LBs, including fundamental soliton (FSs), two-pole soliton (TPSs), four-pole soliton (FPSs) and VSs.

The paper is arranged as follows. In Sec. 2, we introduce the physical model. In Sec. 3, we display the dynamic characteristics of FSs in the system with the PT-symmetric optical potential. Second, we analyze multipole solitons with peaks located at the lattice maxima. Third, the VSs and TPSs with peaks located elsewhere are discussed as well. The last section, Sec.4, summarizes the main results obtained in this work.

## 2. Model

We consider a cold, lifetime-broadened four level atomic system with an inverted-Y type configuration, as shown in Fig. 1(a). To realize an optical PT symmetry in a Rydberg-dressed EIT system, an incoherent population pumping (${\Gamma _{21}}$) coupled to the two low-lying states ${\kern 1pt} |1 \rangle$ and ${\kern 1pt} |2 \rangle$ is introduced [28]. In addition, the control and auxiliary fields are assumed to be periodically modulated. The dynamics of the system is described by the Hamiltonian ${\hat{H}_H}(t) = {N_a}\int_{ - \infty }^{ + \infty } {{d^3}} {\textbf r}{\hat{H}_H}({\textbf r},t)$, where ${N_a}$ is atomic density. Under electric-dipole and rotating-wave approximations, the Hamiltonian density in the interaction picture is:

From Eq. (1), we derive the optical Bloch equations for the density-matrix elements ${\rho _{jl}} \equiv \langle {\hat{S}_{jl}}\rangle $. The dynamics of the density matrix $\hat{\rho }$ is:

In the paper, we take the laser-cooled strontium (88Sr) atomic gas as a realistic candidate for our theoretical model described above. We choose ${R_ \bot } = 12\mu m$, ${\tau _0} = 1.1 \times {10^{ - 6}}s$, ${\Gamma _{21}} = 0.2\pi MHz$, ${\Gamma _3} = 2\pi \times {10^6}MHz$, ${\Gamma _4} = {\Gamma _{34}} = 2\pi \times 16.7kHz$, ${\Delta _2} = 1.67 \times {10^6}{s^{ - 1}}$, ${\Delta _3} = 9.67 \times {10^7}{s^{ - 1}}$, ${\Delta _4} = 1.36 \times {10^7}{s^{ - 1}}$, ${N_a} = 1.0 \times {10^{12}}c{m^{ - 3}}$, ${\Omega _c} = {\Omega _{c0}}\textrm{ = }1.2 \times {10^7}{s^{ - 1}}$, ${\Omega _a} = {\Omega _{a0}}\textrm{ = }5 \times {10^6}{s^{ - 1}}$, ${V_0} = 21.17$, and $\Omega = 4$.

## 3. Local and nonlocal light bullets

We now discuss possible LBs in the system. In 2D cases without nonlinearity, the corresponding Bloch modes are calculated by the Fourier collocation method [31]. The stationary solution of Eq. (4) in the form $\psi = q{e^{ibs}}$(*b* is propagation constant) is numerically found by the squared-operator method (SOM) or the Newton conjugate-gradient method (NCG) [31,32]. Numerically, integration of Eq. (4) is implemented by the spilt-step Fourier method [33]. In addition, soliton’s stability is analyzed by the evolution of perturbed state ${ \psi |_{s = 0}} = q(1 + \sigma f)$, where $\sigma$ is a small perturbation amplitude and $f$ is a random variable uniformly distributed within the interval [−1,1]. We fix $\sigma \textrm{ = }0.02$ and use the same array $f$ for all simulations. In order to monitor the process of evolution, we use the estimated width of the wave packet $W(s) = {\kappa ^{ - 1}}$, where the integral form-factor is $\kappa = \sqrt U /\int {|\psi |d\overrightarrow r }$ and the energy $U = \int {{{|\psi |}^2}d\overrightarrow r } $. Defining the widths of LBs (spatial ${W_1}^2 = \int {{\xi ^2}{{|q |}^2}d\overrightarrow r } /U$, temporal ${W_2}^2 = \int {{\tau ^2}{{|q |}^2}d\overrightarrow r } /U$), and the width ratio $\rho \textrm{ = }{{W(s)} \mathord{\left/ {\vphantom {{W(s)} {W(0)}}} \right.} {W(0)}}$, LBs is identified as stable one at a certain propagation distance ${s_{cr}}$ and a width ratio $\rho $(${s_{cr}} > 50$ and $\rho \textrm{ = }1.2$). Otherwise, LBs is unstable.

#### 3.1 Fundamental solitons

First, we turn to consider dynamic characteristic of FSs. Figure 2(a) shows the band gap spectrum with the 2D PT-symmetric lattice. One can see that the scopes of 3D solitons are closely connected to the domains of the linear spatial Bloch modes in 2D. Two upper bands are displayed with the real component of the PT potential ${p_r} = 6$. The maximal real part $b_{re}^{\max }$ of the propagation constant defines the lower edge of the semi-infinite forbidden gap. When ${p_i} = 3$, the two bands are separated by a forbidden gap. In Figs. 2(e) and 2(g), the existence of stable FSs is owing to minimum $b$ obtained from the lower cutoff ${b_{co}}$($b_{re}^{\max } = 3.1$). When ${p_i} = 7$, it is above the critical value ${p_i} = {p_r} = 6$, hence the real parts merge and the imaginary parts appear in the spectrum, and the broken PT-symmetric phase occurs.

Previous results demonstrated that the stable FSs could be supported by nonlocal nonlinearity [7] or a combination of local nonlinearities and PT potentials [14]. In this paper, we find that the stable FSs can exist with the nonlocal nonlinearity ($\alpha = 1$ and $\gamma = 0.01$) and the PT potential ${p_i} \ne 0$ in a Rydberg-EIT system. From Figs. 2(b)–2(e), one can see that the dynamic characteristics of solitons depend on four parameters (the real/imaginary part of the PT symmetry ${{{p_r}} \mathord{\left/ {\vphantom {{{p_r}} {{p_i}}}} \right.} {{p_i}}}$, and the nonlinear local/nonlocal parameters ${\gamma \mathord{\left/ {\vphantom {\gamma \alpha }} \right.} \alpha }$). Figures 2(b)–2(c) shows *U(b)* dependence at different ${p_r}$ and ${p_i}$. Unlike local nonlinearity case [15], we take $\alpha = 1$ and $\gamma = 0.01$, *U* doesn’t monotonically decrease with bigger *b*. Figure 2(b) indicates that ${p_r}$ increase linearly with *U* when *b* is constant, and decrease with *b* when *U* is constant. On the contrary, ${p_i}$ will lead to an increase of *U* at the same *b*, and an decrease of *b* at the same $U$[see Fig. 2(c)]. In Fig. 2(d), we fix $\alpha \textrm{ = }1$, one can see that an increase of $\gamma$ causes an decrease of the stable interval of *b*. Different with Fig. 2(d), for a fixed $\gamma$, it is shown that an increase of $\alpha$ causes an increase of the stable interval of *b*, thus the FSs are stable here [see Fig. 2(e)].

Based on the parameters given above, one finds propagating velocity [16,34,35] $v \approx 2.12 \times {10^{ - 6}}c$ and the maximum average power density by using Poynting’s vector, which is estimated to be ${\tilde{P}_{\max }} = 1.2n\textrm{W}$[see dot B in Fig. 2(d)]. Thus, this optical pulse moves very slowly and possesses very small power. Figure 2(f) shows the widths of FSs (spatial ${W_1}$, temporal ${W_2}$) as functions of $\alpha$. One finds that for a fixed *b*, the widths increase slowly as $\alpha$ increases, while ${W_2}$ increases faster with $\alpha$ than that of ${W_1}$. Therefore, ${W_1}/{W_2}$ decreases with the increase of $\alpha$. Further, we analyze soliton stability and find an interval of *b* where stable solitons can reside in Fig. 2(g). It is shown that the interval of stability becomes broader and broader with increasing $\alpha$.

Figures 3(a)–3(b) display the distribution of FSs in the PT-symmetric lattice with the local nonlinearity or the nonlocal nonlinearity, respectively [see dots A and B in Figs. 2(d) and 2(e)]. It is shown that FSs are stable with different transverse modulus and vertical modulus width distributions. Figures 3(c)–3(e) show the evolution distributions of FSs with stationary solution corresponding to dots C-E in Figs. 2(d) and 2(e). The width ratios of FSs $\rho (s)$ in Figs. 3(c)–3(f) are plotted with different color in Fig. 2(h). Figure 3(c) shows the distribution of FSs under an unbroken PT-symmetry (${p_i} < {p_r}$) and the local nonlinearity with $\alpha = 0$ and $\gamma = 1$. One can see that the distributions of FSs are localized in transverse direction [Fig. 3(c4)]. In Fig. 2(h), we can see only line 2 corresponds to stable soliton. Due to introduction of imaginary part of the potential ${p_i} \ne 0$, *U* isn’t conserved in evolution. Our simulation shows that even though a numerically accurate stable soliton could also lose *U* and deform its shape after long enough distance evolution by this algorithm. This means, for stable solitons, the shape deformation after long distance comes from numeric error, which explains $\rho $’s increase at big *s* for line 2. In addition, we found $\rho $’s increase rate is different for different $\Delta s$. At relatively small *b*, line 1 in Fig. 2(h) shows that FSs are unstable ($\rho $ exceeds 1.2 at distance $s = 11$). The FSs disperse along $\tau$-axis, but don’t expand along $\xi - \eta$ plane in Fig. 3(c4). When the local nonlinearity is replaced by the nonlocal one with $\alpha = 1$ and $\gamma = 0.01$, we find that the FSs won’t disperse along $\tau$-axis, only expanding along $\xi - \eta$ plane slowly after a long distance $s = 265$ in Fig. 3(d) and line 2 in Fig. 2(h), thus the FSs are stable. When *b* is further increased, FSs split and shift to transverse boundaries quickly [Fig. 3(e)]. Further, the bigger *b* is, the shorter propagation distance ${s_{cr}}$ is. From Fig. 3(d) and Fig. 3(e), one can see that the deformation of spatial and phase distribution in Fig. 3(d) is slower than that in Fig. 3(e). On the other hand, under the broken PT-symmetry (${p_i} > {p_r}$) and the nonlocal nonlinearity, it is shown that the deformation of spatial and phase distribution is quicker even with less perturbation in Fig. 3(f) than that in Figs. 3(d)–3(e). This point is clearly seen from the lines 3 and 4 in Fig. 2(h).

#### 3.2 Multipole solitons with peaks located at lattice maxima

Unlike Ref [14]., we find that the real part of the potential ${p_i} = 0$ and nonlocal nonlinearity can support TPSs and FPSs in a cold Rydberg atomic system. Figures 4(a) and 4(b) show *U* of the function of *b* for TPSs in-phase and FPSs with $l = 1$, respectively. One can see that the stability intervals of the local TPSs and FPSs ($\gamma = 1$ and $\alpha = 0$) are smaller than that of the nonlocal ones ($\gamma = 0.01$ and $\alpha = 1$). Further, Fig. 4(c) shows the evolution of the mass center of the in-phase TPSs for the local nonlinearity ($\gamma = 1$ and $\alpha = 0$) and the nonlocal case ($\gamma = 0.01$ and $\alpha = 1$), respectively. One finds that compared to the local nonlinearity, the oscillation of mass center in nonlocal nonlinearity is smaller, and the duration propagation distances is longer. The solitons oscillate as a “breather” which are also been reported in 2D case [36]. It indicates the nonlocal TPSs is more stable with its shape keeping undestroyed. When *b* is relatively small, this oscillation will last within bounds for longer propagation distances and the energy will transfer between the two poles. Nevertheless, at bigger *b* [$b = 5$ in Fig. 4(c)], the unstable local TPSs will concentrate on one pole and shift to transverse boundaries. In one word, local TPSs have a smaller stability interval of *b*.

Before presenting the intensity distributions of multipole solitons, we display the real part distribution of PT symmetric potential with ${p_r} = 6$ in Fig. 4(d). The lattice maxima is located at ($2{k_1}\pi /\Omega $, $2{k_2}\pi /\Omega $) for ${p_r} > 0$ and ($2{k_1}\pi /\Omega + \pi /\Omega $, $2{k_2}\pi /\Omega + \pi /\Omega $) for ${p_r} < 0$(${k_{1,2}}$ are integers), respectively. Figure 5 shows the intensity distribution of TPSs and FPSs. Comparing Figs. 5(a2)–5(d2) and 5(a3) with Fig. 4(d), we find that soliton’s poles are exactly resided in the lattice maxima. With $\Omega = 4$ and $q = \pi /4$, the coordinates are (${\pm} q$, ${\pm} q$) in Figs. 5(a2)–5(a3) and Fig. 5(b2), (${\pm} 2q$, ${\pm} 2q$) in Fig. 5(c2), and (${\pm} 2q$,0) in Fig. 5(d2), respectively. Figure 5(a1)–5(a2) show a typical TPSs in-phase [see dot $\textrm{A}^{\prime}$ in Fig. 4(a), and Fig. 5(a3) shows another out-of-phase one. The intensity distribution of these TPSs indicate that the latter is less stable than the former. The perturbed state of the two poles respectively shifts along the ${\pm} \tau$ axis, the distant of the two poles are increasing with the propagation distance.

Figures 5(b)–5(d) show three types of the FPSs with $l = 1$, whose peaks are located at different lattice maxima. Figure 5(b) shows FPSs, peaks located at four neighboring lattice maxima ($2 \times 2$), exist with the optical lattice (${p_r} ={-} 6$ and ${p_i} = 0$) and the nonlocal nonlinearity ($\alpha = 1$), [see dot $\textrm{B}^{\prime}$ in Fig. 4(b)]. Unlike the FSs, the FPSs could only be stable in a less interval of *b*. The unstable FPSs will concentrate on one pole and shift to transverse boundaries. Furthermore, we consider a square constructed lattice ($3 \times 3$). One find that solitons with four-poles located at the vertices in Fig. 5(c) or at the midpoint of the square in Fig. 5(d) with small stable interval. Like FSs, ${s_{cr}}$ of all these multipole nonlocal solitons ($\gamma = 0.01$ and $\alpha = 1$) in this section are bigger than that of local one ($\gamma = 1$ and $\alpha = 0$). Hence, they become more stable.

#### 3.3 Vortex and two-pole solitons with peaks located at other places

Now, we investigate the behavior of VSs and TPSs by adjusting ${p_r}$, ${p_i}$, and $\gamma$, respectively. The corresponding $U(b)$ dependences are shown in Figs. 6(a)–6(c). If ${p_r} \ne 0$, both the VSs with $l = 1$ and the TPSs could exist. The VSs with $l = 2$ could not exist. In Fig. 6(a), it is shown that the energy threshold of the VSs increase as ${p_r}$ increases (the red dot line); On the contrary, the energy threshold of the TPSs decreases (the blue solid line). In Fig. 6(b), it is shown that stable TPSs could exist with ${p_i} \ne 0$, while ${p_i}$ results in a linear increase of *U* at the same *b* and an expansion of stable interval of *b*. Furthermore, one can observe a point ${b_{cr}}$ at which ${{dU} \mathord{\left/ {\vphantom {{dU} {db = 0}}} \right.} {db = 0}}$, and ${b_{cr}}$ decreasing with ${p_i}$. However, it is found that VSs with $l = 1$ cannot exist. Figure 6(c) shows $U$ as a function of *b* with changing $\gamma$ for both VSs and TPSs. One finds that *U* decreases faster for TPSs (blue lines) than that for VSs (red dashed lines).

Without PT-symmetric potential, ${p_r}\textrm{ = }{p_i}\textrm{ = 0}$, Eq. (4) becomes a homogenous one, which admits Galilean invariance. Solitons can move in any direction with different propagation constant *b*. We find the perturbation give the VSs a turbulence with a certain factor ${e^{i{v_x}x + i{v_y}y}}$, for different b, the turbulences are different. Figure 6(d) displays the VSs’s trajectory ($l = 1$) of the mass center with perturbation in $\xi - \eta$ plane with $\alpha = 1$, $\gamma = 0.01$ and ${p_r}\textrm{ = }{p_i}\textrm{ = 0}$. One can see that the VSs shift to the boundary along certain direction with its shape almost unchanged at $b = 5$. Further, the velocity keeps constant and the soliton achieves bigger velocity and different deviation angles at a bigger *b*. Figure 6(e) shows some typical evolution of $\rho $ for VSs with ${p_i} = 0$, $\alpha = 1$ and $\gamma = 0.01$. One can see that the corresponding $\rho $ won’t exceed 1.2 for a long propagation distance with ${p_r} = 0$ and $b = 5$, thus VSs are stable. On the other hand, it is unstable with ${p_r} = 0$ and $b = 1$. If more bigger $b = 20$ is chosen, the VSs won’t shift and keep localized for a long propagation distance with ${p_r} = 6$ in Fig. 6(e). The stable distribution of VSs ($\xi - \eta$ plane) is shown in Fig. 7(g) at $b = 7$. Similarly, the system could also support stable TPSs by adjusting ${p_r}$, ${p_i}$ and *b*. Figures 6(f) and 6(g) show the variation $\rho $ of TPSs with $\alpha = 1$ and $\gamma = 0.01$ for different *b*, ${p_{_r}}$ and ${p_i}$, respectively. In Fig. 6(f), it is shown that $b$(smaller $b = 8$ or larger $b = 12.5$) results in ${s_{cr}}$ smaller than that of the intermediate *b* ($b = 9.35$), which is needed to obtain stable TPSs. Hence, TPSs keeps stable only in an specific interval of *b*. In Fig. 6(g), one find that the stable TPSs could not exist with ${p_r} = 0$. Further, Fig. 6(h) shows the stability domains of TPSs in the plane (${p_i}$,b) with $\alpha = 1$, $\gamma = 0.01$ and ${p_r} = 6$. The upper 2 lines display the interval of *b*, one can see that the stability domains increases with ${p_i}$. The lowest line represents a smallest ${b_{cr}}$, where TPSs can exist. It is shown that ${b_{cr}}$ decreases slowly as ${p_i}$ increases.

Next, we display the density profile and the phase distributions of VSs with $l = 1$ and TPSs in Fig. 7. VSs with $l = 1$ are shown in Figs. 7(a)–7(c) and TPSs (out-of-phase) in Figs. 7(d)–7(e) [see A’’-F’’ in Figs. 6(a)–7(b)]. One can see that VSs become more localized with ${p_r} \ne 0$ and a big *b* [see Figs. 6(a), 6(c), 6(d), 6(e) and Figs. 7(a)–7(c)], and the distance between the two poles of TPSs decreases as b increases[see Figs. 7(d), 7(e)]. An interesting feature of the perturbed VS’s evolution is that the soliton will split at small *b* [Fig. 7(h)] and the soliton will shift to the boundary at big b, being beyond the critical value of 4.5∼5 [Fig. 7(g)]. The similar dynamic characteristic in 2D case can be found in Ref [37]. In addition, for ${p_i} \ne 0$, another type of TPSs can exists with peaks located at different places in Fig. 7(i). 2D periodic PT-symmetric potential is the key point to the existence of the LBs in this system. First, stable fundamental solitons can exist under the PT-symmetric potential and nonlocal nonlinearity. Second, the real part of the PT potential and nonlocal nonlinearity can support multipole solitons, VSs, and two-pole solitons.

## 4. Conclusion

In this work, we have proposed a scheme for realizing PT symmetry and nonlocal 3D LBs (fundamental, two-pole, four-pole and vortex solitons) via the Rydberg-dressed EIT. Based on the PT-symmetric lattice and the giant nonlocal Kerr nonlinearity, the system could support PT nonlocal 3D LBs, which have very low light intensity and low velocity. In particular, we have found that the PT symmetry and the degree of nonlocal nonlinearity, which can be actively tuned in our system, produce significant effects on the behavior of LBs. Research presented here opens a route for developing non-Hermitian nonlinear optics, in particular for manipulating LBs with controllable PT-symmetric optical potentials and nonlocal Kerr nonlinearity, which may find applications in optical information processing and transmission.

## Funding

Research and Development Fund, Hubei University of Science and Technology (2019-21GP03); Hubei Technological Innovation Special Fund (2018ABA076, 2019BEC206); National Natural Science Foundation of China (11847103, 11975172, 51976198, 51976199); Qatar National Research Fund (NPRP 8-028-1-001).

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **D. Mihalache, D. Mazilu, F. Lederer, and Y. S. Kivshar, “Collisions between discrete surface spatiotemporal solitons in nonlinear waveguide arrays,” Phys. Rev. A **79**(1), 013811 (2009). [CrossRef]

**2. **C. Q. Dai, Y. Fan, and Y. Y. Wang, “Three-dimensional optical solitons formed by the balance between different-order nonlinearities and high-order dispersion/diffraction in parity-time symmetric potentials,” Nonlinear Dyn. **98**(1), 489–499 (2019). [CrossRef]

**3. **Y. V. Kartashov, B. A. Malomed, and L. Torner, “Solitons in nonlinear lattices,” Rev. Mod. Phys. **83**(1), 247–305 (2011). [CrossRef]

**4. **S. Minardi, F. Eilenberger, Y. V. Kartashov, A. Szameit, U. Röpke, J. Kobelke, K. Schuster, H. Bartelt, S. Nolte, L. Torner, F. Lederer, A. Tünnermann, and T. Pertsc, “Three-dimensional light bullets in arrays of waveguides,” Phys. Rev. Lett. **105**(26), 263901 (2010). [CrossRef]

**5. **Y. S. Kivshar and G. P. Agrawal, “Optical Solitons: From Fibers to Photonic Crystals,”(Academic, 2006).

**6. **O. Lahav, O. Kfir, P. Sidorenko, M. Mutzafi, A. Fleischer, and O. Cohen, “Three-dimensional spatiotemporal pulse-train solitons,” Phys. Rev. X **7**(4), 041051 (2017). [CrossRef]

**7. **F. Maucher, N. Henkel, M. Saffman, W. Królikowski, S. Skupin, and T. Pohl, “Rydberg-induced solitons: three-dimensional self-trapping of matter waves,” Phys. Rev. Lett. **106**(17), 170401 (2011). [CrossRef]

**8. **L. Song, Z. Yang, X. Li, and S. Zhang, “Controllable Gaussian-shaped soliton clusters in strongly nonlocal media,” Opt. Express **26**(15), 19182 (2018). [CrossRef]

**9. **E. L. Falcão-Filho, C. B. de Araújo, G. Boudebs, H. Leblond, and V. Skarka, “Robust two-dimensional spatial solitons in liquid carbon disulfide,” Phys. Rev. Lett. **110**(1), 013901 (2013). [CrossRef]

**10. **L. Torner, S. Carrasco, J. P. Torres, L.-C. Crasovan, and D. Mihalache, “Tandem light bullets,” Opt. Commun. **199**(1-4), 277–281 (2001). [CrossRef]

**11. **L. Torner and Y. V. Kartashov, “Light bullets in optical tandems,” Opt. Lett. **34**(7), 1129 (2009). [CrossRef]

**12. **B. B. Baizakov, B. A. Malomed, and M. Salerno, “Multidimensional solitons in a low-dimensional periodic potential,” Phys. Rev. A **70**(5), 053613 (2004). [CrossRef]

**13. **D. Zhao, S. Ke, Y. Hu, B. Wang, and P. Lu, “Optical bistability in parity-time-symmetric dielectric multilayers incorporated with graphene,” J. Opt. Soc. Am. B **36**(7), 1731 (2019). [CrossRef]

**14. **Y. V. Kartashov, C. Hang, G. Huang, and L. Torner, “Three-dimensional topological solitons in PT-symmetric optical lattices,” Optica **3**(10), 1048 (2016). [CrossRef]

**15. **S. L. Xu, N. Z. Petrovic, M. R. Belic, and Z. L. Hu, “Light bullet supported by parity-time symmetric potential with power-law nonlinearity,” Nonlinear Dyn. **84**(4), 1877–1882 (2016). [CrossRef]

**16. **Z. Bai, W. Li, and G. Huang, “Stable single light bullets and vortices and their active control in cold Rydberg gases,” Optica **6**(3), 309 (2019). [CrossRef]

**17. **I. Friedler, D. Petrosyan, M. Fleischhauer, and G. Kurizki, “Long-range interactions and entanglement of slow single-photon pulses,” Phys. Rev. A **72**(4), 043803 (2005). [CrossRef]

**18. **A. K. Mohapatra, T. R. Jackson, and C. S. Adams, “Coherent optical detection of highly excited Rydberg states using electromagnetically induced transparency,” Phys. Rev. Lett. **98**(11), 113003 (2007). [CrossRef]

**19. **O. Firstenberg, C. S. Adams, and S. Hofferberth, “Nonlinear quantum optics mediated by Rydberg interactions,” J. Phys. B **49**(15), 152003 (2016). [CrossRef]

**20. **A. V. Gorshkov, J. Otterbach, M. Fleischhauer, T. Pohl, and M. D. Lukin, “Photon-photon interactions via Rydberg blockade,” Phys. Rev. Lett. **107**(13), 133602 (2011). [CrossRef]

**21. **S. Baur, D. Tiarks, G. Rempe, and S. Dürr, “Single-photon switch based on Rydberg blockade,” Phys. Rev. Lett. **112**(7), 073901 (2014). [CrossRef]

**22. **E. Distante, A. Padrón-Brito, M. Cristiani, D. Paredes-Barato, and H. de Riedmatten, “Storage enhanced nonlinearities in a cold atomic Rydberg ensemble,” Phys. Rev. Lett. **117**(11), 113001 (2016). [CrossRef]

**23. **Q. Y. Liang, A. V. Venkatramani, S. H. Cantu, T. L. Nicholson, M. J. Gullans, A. V. Gorshkov, J. D. Thompson, C. Chin, M. D. Lukin, and V. Vuletić, “Observation of three-photon bound states in a quantum nonlinear medium,” Science **359**(6377), 783–786 (2018). [CrossRef]

**24. **C. Hang and G. Huang, “Parity-time symmetry along with nonlocal optical solitons and their active controls in a Rydberg atomic gas,” Phys. Rev. A **98**(4), 043840 (2018). [CrossRef]

**25. **Y.-W. Guo, S.-L. Xu, J.-R. He, P. Deng, M. R. Belić, and Y. Zhao, “Transient optical response of cold Rydberg atoms with electromagnetically induced transparency,” Phys. Rev. A **101**(2), 023806 (2020). [CrossRef]

**26. **H. Schempp, G. Günter, C. S. Hofmann, C. Giese, S. D. Saliba, B. D. DePaola, T. Amthor, and M. Weidemüller, “Coherent population trapping with controlled interparticle interactions,” Phys. Rev. Lett. **104**(17), 173602 (2010). [CrossRef]

**27. **S. Sevinçli, N. Henkel, C. Ates, and T. Pohl, “Nonlocal nonlinear optics in cold Rydberg gases,” Phys. Rev. Lett. **107**(15), 153001 (2011). [CrossRef]

**28. **P. Bienias and H. P. Büchler, “Quantum theory of Kerr nonlinearity with Rydberg slow light polaritons,” New J. Phys. **18**(12), 123026 (2016). [CrossRef]

**29. **R. W. Boyd, “Nonlinear Optics,” 3rd Ed. (Academic/Elsevier, 2008)

**30. **W. Demtröder, “Laser Spectroscopy: Basic Concepts and Instrumentation,” 3rd ed. (Springer, Berlin, 2003), Chap. 10.

**31. **J. Yang, “Nonlinear Waves in Integrable and Nonintegrable Systems,” (Society for Industrial and Applied Mathematics, Philadelphia, 2010).

**32. **J. Yang and T. I. Lakoba, “Universally-convergent squared-operator iteration methods for solitary waves in general nonlinear wave equations,” Stud. Appl. Math. **118**(2), 153–197 (2007). [CrossRef]

**33. **W. Bao and D. Jaksch, “An explicit unconditionally stable numerical method for solving damped nonlinear Schrödinger equations with a focusing nonlinearity,” SIAM J. Numer. Anal. **41**(4), 1406–1426 (2003). [CrossRef]

**34. **G. Huang, L. Deng, and M. G. Payne, ““Dynamics of ultraslow optical solitons in a cold three-state atomic system,” Phys. Rev. E **72**(1), 016617 (2005). [CrossRef]

**35. **S. Xu, Q. Zhou, D. Zhao, M. R. Belić, and Y. Zhao, “Spatiotemporal solitons in cold Rydberg atomic gases with Bessel optical lattices,” Appl. Math. Lett. **106**, 106230 (2020). [CrossRef]

**36. **J. Huang, Y. Weng, and H. Wang, “Stability and internal interaction of multipole solitons in nonlocal PT-symmetric lattices,” Opt. Express **26**(9), 11667 (2018). [CrossRef]

**37. **M. Shen, Y. Lin, C. Jeng, and R. K. Lee, “Vortex pairs in nonlocal nonlinear media,” J. Opt. **14**(6), 065204 (2012). [CrossRef]