Stable high-dimensional weak-light soliton molecules and their active control

We start by considering a laser-cooled, dilute three-level atomic gas, interacting with a weak probe laser field $\mathbf{E}_{p}$ with center frequency $\omega_{p}$ (wavenumber $k_{p}=\omega_{p}/c$), driving the transition $|1\rangle\leftrightarrow|2\rangle$, and a strong, continuous-wave control laser field $\mathbf{E}_{c}$ with frequency $\omega_{c}$ (wavenumber $k_{c}=\omega_{c}/c$), driving the transition $|2\rangle\leftrightarrow|3\rangle$; see Fig. 1(a).

Here $|1\rangle$, $|2\rangle$, and $|3\rangle$ denote, respectively, the ground, intermediate, and high-lying Rydberg states; $\Gamma_{12}$ and $\Gamma_{23}$ are spontaneous-emission decay rates from $|2\rangle$ to $|1\rangle$ and from $|3\rangle$ to $|2\rangle$, respectively. The interaction between the two Rydberg atoms respectively at positions $\mathbf{r}$ and $\mathbf{r}^{\prime}$ is described by the van der Waals (vdW) potential

with $C_{6}$ the dispersion coefficient [117, 118]. The initial atomic population is prepared in the ground state $|1\rangle$. The total electric-field vector in the system is given by $\mathbf{E}=\mathbf{E}_{c}+\mathbf{E}_{p}\equiv\sum_{l=c,p}{\mathbf{e}_{l}\mathcal{E}}_{l}\exp[i(\mathbf{k}_{l}\cdot\mathbf{r}-\omega_{l}t)]+\mathrm{c.c}.$, where c.c. stands for the complex conjugate, while $\mathbf{e}_{c}$ and $\mathbf{e}_{p}$ ($\mathcal{E}_{c}$ and $\mathcal{E}_{p}$) are, respectively, polarization unit vectors (envelopes) of the control and probe fields. To suppress the Doppler effect, the probe and control fields are assumed to counter-propagate along the $z$ direction, i.e., $\mathbf{k}_{p}=k_{p}\mathbf{e}_{z}$ and $\mathbf{k}_{c}=-k_{c}\mathbf{e}_{z}$, with $\mathbf{e}_{z}$ being the unit vector of the $z$ direction.

The Hamiltonian of the system is given by $\hat{H}=\mathcal{N}_{a}\int d^{3}{r}\hat{\mathcal{H}}_{0}({\bf r},t)+(\mathcal{N}_{a}/2)\int d^{3}{r}\hat{\mathcal{H}}_{1}({\bf r},t)$. Here $d^{3}r=dxdydz$, $\mathcal{N}_{a}$ is atomic density, $\hat{\mathcal{H}_{0}}({\bf r},t)$ is the Hamiltonian density describing the atoms and the coupling between the atoms and light fields, $\hat{\mathcal{H}_{1}}({\bf r},t)$ is the Hamiltonian density describing the Rydberg-Rydberg interaction. Under the electric-dipole and rotating-wave approximations, $\hat{\mathcal{H}_{0}}$ and $\hat{\mathcal{H}_{1}}$ have the forms

Here $\Delta_{2}=\omega_{p}-(E_{2}-E_{1})/\hbar$ and $\Delta_{3}=\omega_{p}+\omega_{c}-(E_{3}-E_{1})/\hbar$ are, respectively, the one- and two-photon detunings, with $E_{\alpha}$ being the eigenvalue of the energy of the state $|\alpha\rangle$; $\hat{S}_{\alpha\beta}\equiv\hat{\sigma}_{\beta\alpha}\exp\{i[(\mathbf{k}_{\beta}-\mathbf{k}_{\alpha})\cdot\mathbf{r}-(\omega_{\beta}-\omega_{\alpha}+\Delta_{\beta}-\Delta_{\alpha})t]\}$ is the atomic transition operator [124]; $\Omega_{p}={(\mathbf{e}_{p}\cdot\mathbf{p}_{12})\mathcal{E}_{p}}/(2\hbar)$ and $\Omega_{c}={(\mathbf{e}_{c}\cdot\mathbf{p}_{23})\mathcal{E}_{c}}/(2\hbar)$ are, respectively, half Rabi frequencies of the probe and control fields, with $\mathbf{p}_{\alpha\beta}$ the electric-dipole matrix elements associated with the transition $|\beta\rangle\leftrightarrow|\alpha\rangle$. The potential $V$ in the expression of $\hat{\mathcal{H}}_{1}$ is taken as per Eq. (1) [117, 118, 120, 121].

The atomic dynamics is governed by the Heisenberg equation of motion for the operators $\hat{S}_{\alpha\beta}({\bf r},t)$, i.e. $i\hbar\frac{\partial}{\partial t}\hat{S}_{\alpha\beta}({\bf r},t)=[\hat{H},\hat{S}_{\alpha\beta}({\bf r},t)]$. Taking expectation values on the both sides of this equation, we obtain the optical Bloch equation involving one- and two-body reduced density matrices, with the form

where $\hat{\rho}({\bf r},t)$ is reduced one-body density matrix (DM) with matrix elements $\rho_{\alpha\beta}({\bf r},t)\equiv\langle\hat{S}_{\alpha\beta}({\bf r},t)\rangle$, $\Gamma$ is a $3\times 3$ relaxation matrix describing the spontaneous emission and dephasing. Due to the existence of the Rydberg-Rydberg interaction, two-body reduced DM, i.e. $\hat{\rho}_{\rm 2body}(\bf r^{\prime},{\bf r},t)$ [with DM elements $\rho_{\alpha\beta,\mu\nu}({\bf r}^{\prime},{\bf r},t)$], is involved in Eq. (3), represented by the last term $\hat{R}\,[\hat{\rho}_{\rm 2body}]$. The explicit form of Eq. (3) is given in Sec. 1 of Supporting Information.

From Eq. (3), we see that to get the solution of one-body DM elements $\rho_{\alpha\beta}$, equations for two-body elements $\rho_{\alpha\beta,\mu\nu}$ are needed, which, in turn, involve three-body DM element $\rho_{\alpha\beta,\mu\nu,\gamma\delta}$, and so on. As a result, one obtains a hierarchy of infinite equations for $N$-body DM elements ($N=1,\,2,\,3,\cdots$) that must be solved simultaneously. To solve such a chain of equations, a suitable treatment beyond mean-field approximation must be adopted. A powerful one is the reduced density-matrix expansion, by which the hierarchy of the infinite equations is truncated consistently and the problem is reduced to solving a closed system of equations for the one- and two-body DM elements, as elaborated recently [113, 122, 125].

The propagation of the probe field is governed by the Maxwell equation, which, under the slowly-varying-envelope approximation, is reduced to [122]

Here $\nabla_{\perp}^{2}=\partial_{x}^{2}+\partial_{y}^{2}$ and $\kappa_{12}=\mathcal{N}_{a}\omega_{p}|\mathbf{p}_{12}|^{2}/(2\epsilon_{0}c\hbar)$.

Since the probe field is much weaker than the control field, the depletion of the atomic population in the ground state is small and a standard perturbation method can be applied to solve the system of Maxwell-Bloch (MB) equations (3) and (4). To include the many-body correlations produced by the strong Rydberg-Rydberg interactions in a reasonable way, a beyond mean-field approximation [113, 122, 125] mentioned above must be used. Then, in the leading-order approximation, we obtain the solution for the probe field $\Omega_{p}=F(x,y,z,t)e^{iK(\omega)z-i\omega t}$, where $F$ denotes a slowly varying envelope function and $K(\omega)$ stands for the linear dispersion relation, $K(\omega)=\omega/c+\kappa_{12}(\omega+d_{31})/[|\Omega_{c}|^{2}-(\omega+d_{21})(\omega+d_{31})]$.

At the third-order approximation (details are given in Sec. 2 of the Supporting Information), we obtain the nonlinear equation for the probe-field envelope

where $T=t-z/V_{g}$, with $V_{g}=(\partial K/\partial\omega)^{-1}$ being the group velocity of the envelope, and $K_{2}=\partial^{2}K/\partial\omega^{2}$ defines the group-velocity dispersion. The last two terms in Eq. (5) contributed, respectively, from the local and nonlocal optical Kerr nonlinearities in the system. Explicit expressions of the coefficients $W$ and $G$ (nonlocal nonlinear response function) are given in Sec. 2 of the Supporting Information. The local optical Kerr nonlinearity is contributed to by the short-range interaction between photons and atoms, proportional to the atomic density $\mathcal{N}_{a}$, whereas the nonlocal optical Kerr nonlinearity is contributed by the long-range Rydberg-Rydberg interaction, which scaled quadratically with the atomic density (i.e., proportional to $\mathcal{N}_{a}^{2}$).

Equation (5) can be further written into the non-dimensional form

where $u=\Omega_{p}/U_{0}$, $\zeta=z/(2L_{\mathrm{diff}})$, $(\xi,\eta)=(x,y)/R_{0}$, $\tau=T/\tau_{0}$, $d=-\mathrm{sgn}(K_{2})L_{\mathrm{diff}}/L_{\mathrm{disp}}$, $w=2L_{\mathrm{diff}}|U_{0}|^{2}W$, and $g=2L_{\mathrm{diff}}R_{0}^{2}|U_{0}|^{2}G(\xi^{\prime}-\xi,\eta^{\prime}-\eta)$, with the diffraction and dispersion lengths respectively given by $L_{\mathrm{diff}}=\omega_{p}R_{0}^{2}/c$ and $L_{\mathrm{disp}}=\tau_{0}^{2}/|K_{2}|$. Here $U_{0}$, $R_{0}$, and $\tau_{0}$ are typical half Rabi frequency, beam radius in the transverse plane ($x$,$y$), and temporal duration of the probe-field envelope, respectively.

To address a typical example, we consider laser-cooled strontium ${}^{88}\mathrm{Sr}$ atoms, with atomic levels $|1\rangle=\left|5s^{2}~{}^{1}S_{0}\right\rangle$, $|2\rangle=\left|5s5p~{}^{1}P_{1}\right\rangle$, and $|3\rangle=\left|5sns~{}^{1}S_{0}\right\rangle$. For the principal quantum number $n=60$, the dispersion parameter $C_{6}\approx 2\pi\times 81.6\,\mathrm{GHz}\cdot\mu\mathrm{m}^{6}$ (which implies the Rydberg-Rydberg interaction is attractive). The spontaneous emission decay rates are $\Gamma_{12}\approx 2\pi\times 32$ MHz and $\Gamma_{23}\approx 2\pi\times 16.7$ kHz, while the detunings are taken to be $\Delta_{2}=-2\pi\times 240$ MHz and $\Delta_{3}=-2\pi\times 0.16$ MHz. The density of the atomic gas is $\mathcal{N}_{a}=9\times 10^{10}~{}\mathrm{cm}^{-3}$, and the half Rabi frequency of the control field is $\Omega_{c}=2\pi\times 5$ MHz. Since $\Delta_{2}\gg\Gamma_{12},\,\Delta_{3}$, which makes the system works in a dispersive nonlinearity regime, the imaginary parts of coefficients in Eq. (5) are much smaller than the corresponding real parts, and hence Eq. (6) can be approximately considered as a real-coefficient one.

The nonlocal nonlinear response function $g(\xi^{\prime}-\xi,\eta^{\prime}-\eta)$ has a very complicated expression. For the convenience of the subsequent variational calculation for the interaction force between solitons, we approximate it by a Gaussian function, i.e.,

(details of relations between $g$ and $g_{F}$ are given in Sec. 3 of the Supporting Information), where $g_{0}=\iint d\xi d\eta g(\xi^{\prime}-\xi,\eta-\eta)$ is a constant and $\sigma$ characterizes the nonlocality degree of the nonlinearity, defined as

Here $R_{b}=(|C_{6}/\delta_{\mathrm{EIT}}|)^{1/6}$ denotes the Rydberg blockade radius, where $\delta_{\mathrm{EIT}}\approx|\Omega_{c}|^{2}/|\Delta_{2}|$ is the linewidth of the EIT transmission spectrum for $|\Delta_{2}|\gg\Gamma_{12}$ [120, 121]. With the values of parameters adopted above, we get $R_{b}\approx 9.6\,\mu\mathrm{m}$. If $R_{b}\ll R_{0}$ (i.e., $\sigma\rightarrow 0$, the local limit), the nonlocal response reduces to the delta function, i.e., $g_{0}\delta(\xi-\xi^{\prime},\eta-\eta^{\prime})$. In this limit, the nonlocal Kerr nonlinearity is reduced to the usual local term, $g_{0}|u|^{2}u$. If $R_{b}\gg R_{0}$ ($\sigma\rightarrow\infty$, the strongly nonlocal limit), the nonlocal response, given by Eq. (7), reduces to a linear term $g_{0}Pu$, where $P=\iint d\xi d\eta\,|u|^{2}$ is the power of the probe field (this limit corresponds to the so-called “accessible-soliton” model, which is actually a limit one [98, 100]).

The susceptibility of the probe field is defined by $\chi_{p}=\mathcal{N}_{a}(\mathbf{e}_{p}\cdot\mathbf{p}_{12})^{2}\rho_{21}/(\varepsilon_{0}\hbar\Omega_{p})$, which can be further expressed as $\chi_{p}\approx\chi^{(1)}+\chi_{\mathrm{loc}}^{(3)}|\mathcal{E}_{p}|^{2}+\chi_{\mathrm{nloc}}^{(3)}|\mathcal{E}_{p}|^{2}$. In this expansion, $\chi^{(1)}$ represents the linear susceptibility; $\chi_{\mathrm{loc}}^{(3)}$ and $\chi_{\mathrm{nloc}}^{(3)}$ are respectively the local and nonlocal third-order nonlinear susceptibilities, associated to the coefficients of Eq. (5) as

With the system’s parameters introduced above, we obtain $\chi_{\mathrm{loc}}^{(3)}\sim 10^{-11}~{}\mathrm{m^{2}V^{-2}}$ and $\chi_{\mathrm{nloc}}^{(3)}\sim 10^{-8}~{}\mathrm{m^{2}V^{-2}}$, i.e., the nonlocal optical Kerr nonlinearity is three orders of magnitude stronger than the local one due to the Rydberg-Rydberg interaction.

We now address the interaction force between two nonlocal (2+1)D [126] spatial solitons. By setting $d=0$ (which implies that time duration $\tau_{0}$ of the probe field is large enough to make the dispersion of the system negligible, valid for $L_{\mathrm{disp}}\gg L_{\mathrm{diff}}$), Eq. (6) reduces to the form

The Lagrangian corresponding to this equation is $L=\iint_{-\infty}^{\infty}d\xi d\eta\,\mathcal{L}$, with the density $\mathcal{L}=\frac{i}{2}(uu_{\zeta}^{\ast}-u^{\ast}u_{\zeta})+|u_{\xi}|^{2}+|u_{\eta}|^{2}-\frac{w}{2}|u|^{4}-\frac{1}{2}|u|^{2}\iint g(\xi-\xi^{\prime},\eta-\eta^{\prime})\left|u\left(\xi^{\prime},\eta^{\prime},\zeta\right)\right|^{2}d\xi d\eta$.

The bound state of two solitons (i.e., the two-soliton “molecule”) is sought by means of the ansatz

with each soliton approximated by a Gaussian, i.e.,

It includes two identical Gaussian beams, located at different positions $(-\mathcal{D}/2,0)$ and $(\mathcal{D}/2,0)$, with a $\pi$ phase difference. Variational $z$-dependent parameters are $A$ (amplitude), $\mathcal{D}$ (spatial separation), $b$ (chirp), and $\phi$ (phase), while the total power $P=2\pi a^{2}A^{2}[1-e^{-\mathcal{D}^{2}/(4a^{2})}]$ is a conserved quantity. Due to the $\pi$ phase difference, there is a repulsive interaction between the two solitons. Such a repulsive interaction is expected to balance the attractive interaction induced by the self-focusing nonlocal optical Kerr nonlinearity (induced by the attractive Rydberg-Rydberg interactions), and thus help to form a stable two-soliton molecule [127].

Following the standard procedure of the variational approximation [113, 127, 128, 129], one can derive evolution equations for the variational parameters $A$, $\mathcal{D}$, $b$, and $\phi$, and hence the equation of motion for the spatial separation $\mathcal{D}$ between center-of-mass positions of the interacting solitons

This equation can be cast in the form of the equation of motion for the Newtonian particle $M_{s}(d^{2}\mathcal{D}/d\zeta^{2})=-\partial U/\partial\mathcal{D}$, where $M_{s}$ is the effective mass of each soliton and $U=U(\mathcal{D})$ denotes the effective potential which accounts for the interaction between the two solitons.

We recall that in locally nonlinear media, the interaction between two solitons is determined by the overlap between their wave functions, which quickly decays as the separation between them increases. Normally, the interaction becomes negligible if the separation between the solitons is $2\sim 3$ times greater than their widths [114]. Nevertheless, for nonlocally nonlinear media the interaction between two solitons takes place when they have no tangible overlap, which may be called contactless interaction. To demonstrate this clearly, we introduce an overlap parameter

where $u_{+}$ and $u_{-}$ are introduced in Eq. (12). Figure 2(a)

shows $J$ as a function of $\sigma$ for the ansatz (11), with system’s parameters $A=2$, $b=0$, $a=1$, and $\phi=0$. We see that $J$ is non-vanishing only for $\sigma\leqslant 0.2$.

For locally nonlinear media (for which $\sigma\sim 0$), the interaction between two solitons is negligible when there is no overlap between the solitons. However, the interaction between the two solitons is non-zero even they has no overlap if the nonlocality degree of the Kerr nonlinearity reaches to a critical value [i.e., $\sigma\geqslant 0.2$ in the present Rydberg gas]. To characterize the size of the SM, we define the separation-width ratio

i.e., the ratio of the separation of the two solitons, $\mathcal{D}$, to the width of a single soliton, $2a$ [see Fig. 2(a)]. Shown in Fig. 2(b) is the effective interaction potential $U$  (calculated following Ref. [127]) as a function of $\lambda$ for the ansatz (11) with $\sigma=1.8$, $A=2$, $b=0$, $a=1$, and $\phi=0$. It is seen that $U$ has a minimum $U_{\mathrm{min}}$ at $\lambda\approx 4.6$, which corresponds to separation $\mathcal{D}=\mathcal{D}_{0}\approx 9.2$ between the two solitons. Thus, it is possible to expect the existence of the SM with size close to $\mathcal{D}_{0}$. Because the width of each soliton in the SM is $2a=2$, and the separation between them at the equilibrium position is $\mathcal{D}_{0}=9.2$, the effective soliton interaction is indeed contactless. The reason for the occurrence of such a contactless interaction between the solitons is due to the giant nonlocal Kerr nonlinearity, which induces significant interaction between the solitons while their wave functions exhibit no overlap in space. On the other hand, in media with the local Kerr optical nonlinearity, the interaction between solitons, and hence the formation of SMs, is determined by the overlap of “tails” of the wave functions of adjacent solitons, and becomes negligible when there is no overlap between the solitons. Consequently, SMs in locally nonlinear media usually have only a small separation-width ratio.

Shown in the inset to Fig. 2(b) is the equilibrium separation $\mathcal{D}_{0}$ between the two solitons as a function of the degree of the nonlocality of the Kerr nonlinearity $\sigma$. It is seen that, as $\sigma$ increases, $\mathcal{D}_{0}$ grows at first, reaches its maximum, and decreases, eventually saturating to a small value. Thus, by changing $\sigma$, one can control $\mathcal{D}_{0}$ and thus the size of the SM. The maximum value of $\mathcal{D}_{0}$ is obtained at $\sigma\approx 1.6$, which is about 13.6, and the corresponding size of the SM in physical units is $13.6\,R_{0}\approx 82$ $\mu$m. Such a value is a realistic one for the experimental observation in Rydberg gases [130].

Stable (2+1)D spatial SMs may also be composed of vortex solitons. As an example, we consider a bound state of two vortex solitons built as per ansatz (11), with each component $u_{\pm}$ taken as the Laguerre-Gaussian beam

Here $A$ and $a$ are the soliton’s amplitude and radius, $r_{\pm}=\sqrt{(\xi\pm\mathcal{D}/2)^{2}+\eta^{2}}$, $L_{p}^{|l|}$ is the generalized Laguerre-Gaussian polynomial, with the azimuthal (radial) index $l$ ($p$), and $\varphi$ is the azimuthal angle. The ansatz based on Eqs. (11) and (19) introduces the superposition of two Laguerre-Gaussian beams with identical shapes, opposite signs, and centers placed at points $(\pm\mathcal{D}/2,0)$.

Figure 5(a) [see also Fig. 1(b)] shows the propagation of a typical (2+1)D two-vortex molecule. Here, we fix $l=1$ and $p=0$ in the Laguerre-Gaussian polynomial in Eq. (19). The input, composed as per Eqs. (11) and (19), includes the perturbation factor (16) too. Parameters of the input are $A=3.9$, $a=1.45$, $\mathcal{D}=10$, and $\epsilon=0.05$. The two-vortex molecule has a large size, with the equilibrium separation between pivots of the two vortices $\mathcal{D}_{0}=10$, corresponding to 60 $\mu\mathrm{m}$ in physical units (in experiments with BEC, “large” is usually a size which is essentially larger than 10 $\mu\mathrm{m}$ [133]). Such a contactless interaction between the two vortices and the formation of the vortex molecule is also due to the nonlocal Kerr nonlinearity contributed by the long-range Rydberg-Rydberg interaction between the atoms.

We have also simulated the propagation of a three-vortex molecule, as shown in Fig. 5(b). As well as its two-vortex counterpart, it is found to be stable. The fourth column of Fig. 5(a) and (b) show the same outcome of the evolution, except that we have added initial tangential velocities $v=0.6$ to the vortex, hence the vortex molecule exhibits persistent rotation, keeping its stability in the course of the propagation. The rightmost panel shows the phase distribution in the (2+1)D spatial vortex molecules at $\zeta=2$.

The realization of (3+1)D [134] spatiotemporal solitons is a long-standing challenging goal of optical physics [82, 87, 88]. As mentioned above (3+1)D spatiotemporal solitons are strongly unstable in conventional optical media with the local Kerr nonlinearity. In a recent work, it has been shown that stable (3+1)D spatiotemporal solitons may exist in a cold Rydberg atomic gas, being supported by a two-step self-trapping mechanism [113].

To proceed, we first demonstrate that stable (3+1)D spatiotemporal soliton/vortex molecules are available in the Rydberg atomic gas. To form such states, dispersion of the probe field is necessary, which can be secured by using a probe pulse with a short time duration. To this end, we adopt a new set of system’s parameters: $\mathcal{N}_{a}=10^{11}$ cm^-3, $\Delta_{2}=-2\pi\times 240$ MHz, $\Delta_{3}=-2\pi\times 0.03$ MHz, $\Omega_{c}=2\pi\times 8$ MHz, and $\tau_{0}=0.1\,\mu\mathrm{s}$. With these parameters, the scaled coefficients in Eq. (6) are $d\approx 0.19$ and $w\approx 0.25$. The (3+1)D spatiotemporal two-soliton molecules can be sought by using ansatz (12) for a single soliton, multiplying it by the temporal-localization factor, $\mathrm{sech}(\tau/a_{\tau})\exp\left(ib_{\tau}\tau^{2}\right)$, where $a_{\tau}$ and $b_{\tau}$ stand for the temporal width and chirp of the probe pulse.

Shown in Fig. 6(a)

is the propagation of a typical (3+1)D spatiotemporal two-soliton molecule. The initial condition used in the numerical simulation again includes the small random perturbation. Parameters of the input are $A=1.6$, $\mathcal{D}=8$, $a=1.2$, $a_{\tau}=1$, $b=b_{\tau}=\phi=0$, and $\epsilon=0.05$. The (3+1)D spatiotemporal SM is found to be stable in the course of the propagation.

With the parameters given above, the propagation velocity of the (3+1)D spatiotemporal SM, produced by the formula $V_{g}=(\partial K/\partial\omega)^{-1}$ at $\omega=0$, is

and the required generation power is estimated to be $P_{\mathrm{gen}}=6.8\,\mu\mathrm{W}$. Thus, the SM travels indeed with an ultraslow velocity (in comparison to $c$) and may be created by a very low power, which is due to the interplay of the EIT and giant nonlocal Kerr nonlinearity induced by the Rydberg-Rydberg interaction between the atoms.

We have also carried out a numerical simulation for the propagation of a nonlocal (3+1)D spatiotemporal two-vortex molecule, with individual vortex solitons taken as per ansatz (19) times $\mathrm{sech}(\tau/a_{\tau})\exp\left(ib_{\tau}\tau^{2}\right)$. Fig. 6(b) shows the propagation of a typical two-vortex SM with small perturbations. The parameters used in the simulation are $A=2.5$, $a=1.3$, $a_{\tau}=1$, $\mathcal{D}=8$, $b=b_{\tau}=\phi=0$, and $\epsilon=0.05$. The two-vortex SM is also found to be quite stable, as well as the zero-vorticity spatiotemporal SM.

Keeping memory of optical pulses in atomic gases (i.e., storage of incident pulses, with ability to retrieve them), provided by the EIT technique, has attracted much interest [135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146]. Here we demonstrate that the storage and retrieval of (3+1)D spatiotemporal soliton/vortex molecules are possible in the present Rydberg-EIT system. To this end, we investigate the evolution of the (3+1)D spatiotemporal molecules, considered above, by solving the MB Eqs. (3) and (4) numerically, using a time-dependent control field

which provides switching action for the probe field. Here $T_{\mathrm{off}}$ and $T_{\mathrm{on}}$ are, respectively, the times at which the control field is switched off and on. The switching duration is $T_{s}$ and the storage time is $T_{\mathrm{on}}-T_{\mathrm{off}}$. Importantly, the switching speed of the control laser has a marginal effect on the quality of the SM storage and retrieval. This is because the (3+1)D spatiotemporal SM in the present system propagates with an ultraslow velocity, due to the EIT effect [147].

Shown in Fig. 7(a)

is the evolution of the probe-pulse amplitude $|\Omega_{p}\tau_{0}|$ in the course of the storage and retrieval process. The shapes of the probe pulse at $z=0$ (before storage), $z=4L_{\mathrm{diff}}$ (at the beginning of storage), and $z=8L_{\mathrm{diff}}$ (after the retrieval) are shown for $L_{\mathrm{diff}}=0.87$ mm. It is seen that switching off the control field provides for the storage of the (3+1)D spatiotemporal SM in the atomic medium, which is retrieved when the control field is switched on again. Further, the retrieved spatiotemporal SM has nearly the same shape as the original one prior to the storage. In the course of the storage, the information carried by the SM is converted into that kept in the atomic spin wave (i.e., the coherence matrix element $\rho_{13}$). A slight deformation affecting the optical memory is due to dissipation, including spontaneous emission and dephasing, as well as weak imbalance between diffraction, dispersion, and nonlinearity. We have also explored the storage and retrieval of (3+1)D spatiotemporal vortex molecules. Similar results are obtained for the vortex bound states, as shown in Fig. 7(b).

The quality of the storage and retrieval of nonlocal (3+1)D spatiotemporal SMs can be characterized by efficiency $\eta$ and fidelity $\eta{\mathcal{J}}$, where $\eta$ and ${\mathcal{J}}$ are defined as

Based on the results obtained in Fig. 7(a) and (b), we obtain $\eta=90.39\%$, ${\mathcal{J}}=98.81\%$, and $\eta{\mathcal{J}}=89.31\%$ for the (3+1)D spatiotemporal SM, and $\eta=90.32\%$, ${\mathcal{J}}=97.48\%$, and $\eta{\mathcal{J}}=88.05\%$ for the vortex molecule.

The strength of the nonlocal Kerr nonlinearity has a significant effect on the quality of the optical memory. Figure 7(c) shows fidelity $\eta{\mathcal{J}}$ of the retrieved (3+1)D spatiotemporal SM as a function of the probe-pulse amplitude. For this purpose, a set of probe-pulse isosurfaces, $|\Omega_{p}\tau_{0}|=3$, $11$, and $19$ at $z=8L_{\mathrm{diff}}\approx 7$ mm are displayed. For the moderate amplitude, $|\Omega_{p}\tau_{0}|\approx 11$, the fidelity reaches its maximum, with the retrieved spatiotemporal SM having nearly the same shape as the original one prior to the storage. For small and large amplitudes, the fidelity features small values, i.e., the retrieved SM is distorted greatly. This happens because, for the weak and strong probe-pulse amplitudes, the Kerr nonlinearity is either too weak or too strong to balance the diffraction and dispersion.

The nonlocality degree of the Kerr nonlinearity also has a significant effect on the memory quality for the spatiotemporal SMs. Figure 7(d) shows fidelity $\eta{\mathcal{J}}$ of the retrieved (3+1)D SM as a function of the nonlocality degree, $\sigma$. The probe-pulse isosurfaces ($|\Omega_{p}\tau_{0}|=0.5$) are displayed for $\sigma=0.5$, $1.4$, and $5.0$ at $z=8L_{\mathrm{diff}}\approx 7$ mm. The fidelity reaches its maximum in the case of the moderate nonlocality degree, $\sigma\approx 1.4$, letting the retrieved SM keep nearly the same shape as the original one had. In the cases of small and large nonlocality degrees, the fidelity may have only small values, greatly distorting the retrieved SM. This happens because, in the limit of local response ($\sigma\rightarrow 0$), the Kerr nonlinearity becomes local, making all (3+1)D solitons unstable, as mentioned above. On the other hand, in the limit of the strongly nonlocal response ($\sigma\rightarrow\infty$), the nonlocal Kerr nonlinearity reduces to a linear potential (as in the above-mentioned “accessible-soliton” model [98]), which cannot support stable (3+1)D SMs either.