Optical Pattern Formation in a Rydberg-Dressed Atomic Gas with Non-Hermitian Potentials

Abstract

Spontaneous pattern formation from a spatially homogeneous background of nonlinear systems driven out of equilibrium is a widespread phenomenon in nature. However, similar phenomena and their physical realization in nonlinear systems with external potentials of gain and loss remain a challenge. We propose a scheme to realize a new type of spatial pattern formation through the self-organization of laser light in a Rydberg-dressed atomic gas with self-defocusing Kerr nonlinearity as well as non-Hermitian optical potentials. We show that by a suitable design of control and assistant laser fields, non-Hermitian optical potentials with or without parity-time ($\mathcal{PT}$) symmetry for the probe laser field can be created. We find that through the nonlocal Kerr nonlinearity contributed by the long-range atom–atom interaction, a constant-intensity wave (CIW) may undergo modulation instability and induce spontaneous symmetry breaking, resulting in the emergence of various self-organized optical structures, which can be actively manipulated by tuning the nonlocality degree of the Kerr nonlinearity and by designing the non-Hermitian optical potentials. The results reported here open a door for developing non-Hermitian nonlinear optics.

1. Introduction

The dynamical instability of spatially continuous systems is one of the most fascinating scientific research topics. It has been found that many systems display fascinating spatiotemporal structures. Very well-known instability mechanisms have been established over the last few decades, such as Rayleigh–Benard convection, Taylor–Couette flow, and Benjamin–Feir instability [1,2,3], occurring in a broad range such as biology, chemistry, hydrodynamics, and soft-matter physics [4,5,6,7,8,9]. One such dynamical instability is modulation instability (MI), which is studied widely in parity-time ($\mathcal{PT}$) symmetric optical systems [10,11,12] and local and/or nonlocal Kerr nonlinear, even competing, nonlinear systems [13,14,15,16,17,18,19,20].

MI induces a spontaneous symmetry breaking (SSB) in both classical and quantum mechanics, widely happening in dual-core Bragg gratings [21,22], Bose–Einstein condensates [23], quantum field theory [24], and spin waves [25]. It not only is a destabilization mechanism to produce a self-induced breakup of an initially continuous wave into localized structures (i.e., soliton) in self-focusing nonlinear medium [26,27,28,29], but also provides an extended structure (ordered pattern) in self-defocusing optical medium from the homogeneous state [30,31,32].

In addition to SSB in general Hamiltonian systems, SSB in a non-Hermitian system also has been focused on the research in various physical systems. The primary motivation for such research is to develop non-Hermitian quantum theory [33,34,35,36]. With the advance of the research, it has been found that $\mathcal{PT}$-symmetric optics, or more generally a non-Hermitian flat-band system [37,38], may have many practical applications, such as the realization of unidirectional light propagations [39], coherent perfect absorbers [40,41,42], giant light amplification [43], novel lasers [44,45], precision measurement [46,47], and quantum computation [48], and etc.

Furthermore, in the past decade, a large amount of research was focused on the investigation of cold Rydberg atomic gases [49,50,51,52,53,54,55,56,57,58,59] working under a condition of electromagnetically induced transparency (EIT), such as quantum information [50], quantum many-body simulation [51], Rydberg spectroscopy [52], Rydberg blockade effect [53], Rydberg self-organization criticality [54,55], Rydberg sensing experiments for multi-frequency microwave recognition [56], highly sensitive measurement [57], and spectrometer [58], etc. EIT is an important quantum destruction interference effect typically occurring in resonant three-level atomic systems, by which the absorption of a probe laser field can be largely eliminated by using a control laser field [60]. In an interesting work, Sevincli et al. [61] reported a self-organized hexagonal optical pattern via an MI of plane-wave probe beam in a cold Rydberg atomic gas with a repulsive Rydberg–Rydberg interaction. Recently, it was demonstrated that a structural phase transition of optical patterns from a hexagonal lattice to two types of square lattices may occur in an EIT-based Rydberg gas with a microwave dressing between two Rydberg states [62].

In this work, combining the $\mathcal{PT}$ symmetry with widely concerned Rydberg atomic gas, we propose an inverted-$\mathrm{Y}$ type scheme to realize the self-organization of laser light with physically realistic non-Hermitian optical potentials by using a coherent Rydberg atomic gas through a spatial modulation of control laser field and the inclusion of the nonlocal Kerr nonlinearity of the probe laser field. We first show the MI with both $\mathcal{PT}$-symmetric and non-$\mathcal{PT}$-symmetric potentials in detail. Based on the MI, SSB happens and subsequently results in the self-organization of laser light. The rich structures appear from a constant-intensity wave (CIW), such as ring-like waves, and tetragonal and hexagonal pattern waves, all of which can be manipulated by the appropriate system parameters.

Before proceeding, we notice that in the literature there were many studies devoted to MI in general local and (or) nonlocal nonlinear optical systems and $\mathcal{PT}$-symmetric optical systems [11,12,13,14,15,16,17,21,22,36,41]. However, in many aspects, our work is different from those. First, in most of those Refs., no realistic physical scheme was provided on how to construct and realize non-Hermitian potentials; our work presents a practical scheme for realizing the non-Hermitian optical potentials, which is useful to guide realistic experiments. Second, no results on nonlocal Rydberg–Rydberg interaction are given. In our work, such interaction, especially self-defocusing nonlocal nonlinearity, is studied, and the result reveals that the property of the nonlocal structure (ordered pattern) is strongly dependent on the shape of the nonlocal nonlinear potential. Third, there is no report about two-dimensional (2D) MI, or the combination of optical potential with nonlocal nonlinearity is not considered. In our work, such a situation is investigated in detail, and interesting results are obtained. We consider all kinds of optical potential which include periodic potentials, parabolic potential, and ring-like optical potential. The results reported here open a route for exploring the SSB and optical self-organization based on cold Rydberg gases.

The article is organized as follows. In Section 2, the physical model of the Rydberg-EIT for inverted-Y configuration is considered. In Section 3, the 2D $\mathcal{PT}$-symmetry and non-$\mathcal{PT}$-symmetry potential are realized. The detailed MI analysis and optical pattern formations of the laser light for $\mathcal{PT}$-symmetric and non-$\mathcal{PT}$-symmetric potential are given in Section 4 and Section 5, respectively. Our conclusions are given in Section 6.

2. Model

We consider an ensemble of lifetime-broadened four-level atomic gas with an inverted-Y type level configuration, shown schematically in Figure 1a. The weak probe laser field with central angular frequency ${\omega }_p$, wavevector ${\mathbf{k}}_p$, and half Rabi frequency ${\Omega }_p$ drives the transition $|1\rangle \leftrightarrow |3\rangle$. The strong control (dressed) laser field with central angular frequency ${\mathsf{\omega }}_c$ (${\mathsf{\omega }}_d$), wavevector ${\mathbf{k}}_c$ (${\mathbf{k}}_d$), and half Rabi frequency ${\Omega }_c$ (${\Omega }_d$) drives the transition $|2\rangle \leftrightarrow |3\rangle$ ($|3\rangle \leftrightarrow |4\rangle$), respectively. Therefore, the total electric fields acting in the atomic system can be written as $\mathbf{E}(\mathbf{r},t)={\sum }_{\alpha }{\mathbf{e}}_{\alpha }{\mathcal{E}}_{\alpha }exp[i({\mathbf{k}}_{\alpha }·\mathbf{r}-{\omega }_{\alpha }t)+\mathrm{c}.\mathrm{c}.](\alpha =p,c,d)$, with ${\mathbf{e}}_{\alpha },{\mathcal{E}}_{\alpha }$ the unit vector of the polarization direction and the envelope of the corresponding laser fields, respectively. States $|1\rangle$ and $|2\rangle$ are the two ground states, and state $|3\rangle$ is the ordinary excited state. The $\mathsf{\Lambda }$-type EIT (levels $|1\rangle$, $|2\rangle$, and $|3\rangle$) is dressed by highly excited Rydberg state $|4\rangle$. Under the electric-dipole and rotating-wave approximations, the Hamiltonian density of the system reads

$$\begin{array}{cc}\hfill {\mathcal{H}}_H=& \sum _{\alpha =2}^4{\mathsf{\Delta }}_{\alpha }{\widehat{S}}_{\alpha \alpha }(\mathbf{r},t)-\hslash \left({\Omega }_p{\widehat{S}}_{13}+{\Omega }_c{\widehat{S}}_{23}+{\Omega }_d{\widehat{S}}_{34}+\mathrm{h}.\mathrm{c}.\right)\hfill \\ \hfill & +\frac{1}{2}{\mathcal{N}}_a\int d^3{\mathbf{r}}^{\prime }{\widehat{S}}_{44}(\mathbf{r},t)\hslash V_{\mathrm{vdw}}({\mathbf{r}}^{\prime }-\mathbf{r}){\widehat{S}}_{44}({\mathbf{r}}^{\prime },t),\hfill \end{array}$$

(1)

where ${\mathcal{N}}_a$ is the atomic density, and ${\mathsf{\Delta }}_{\alpha }$ is the detuning of state $|\alpha \rangle (\alpha =2,3,4)$. The half Rabi frequencies of the probe, control, and dressed fields are, respectively, ${\Omega }_p={\mathbf{e}}_p·{\mathbf{p}}_{31}{\mathcal{E}}_p/\hslash$, ${\Omega }_c={\mathbf{e}}_c·{\mathbf{p}}_{32}{\mathcal{E}}_c/\hslash$, and ${\Omega }_d={\mathbf{e}}_d·{\mathbf{p}}_{43}{\mathcal{E}}_d/\hslash$, with ${\mathbf{p}}_{\alpha \beta }$ the electric dipole matrix element associated with the transition from $|\beta \rangle$ to $|\alpha \rangle$, and ${\widehat{S}}_{\alpha \beta }=|\beta \rangle \langle \alpha |exp\left\{i[({\mathbf{k}}_{\beta }-{\mathbf{k}}_{\alpha })·\mathbf{r}-({\omega }_{\beta }-{\omega }_{\alpha }+{\mathsf{\Delta }}_{\beta }-{\mathsf{\Delta }}_{\alpha })t]\right\}$ the transition operator. The last term on the right-hand side of Equation (1) is the contribution of the Rydberg–Rydberg interaction with $\hslash V_{\mathrm{vdw}}({\mathbf{r}}^{\prime }-\mathbf{r})\equiv -\hslash C_6/{|{\mathbf{r}}^{\prime }-\mathbf{r}|}^6$ the long-range interaction potential and $C_6$ the dispersive parameter. Typically, for ${}^{87}\mathrm{Rb}$, the negative $C_6$ results in the repulsive interaction.

According to the Hamiltonian Equation (1), the dynamics of atoms in the system are governed by the optical Bloch equation, given by [63,64,65]

$$\frac{\partial \rho }{\partial t}=-\frac{i}{\hslash }\left[{\mathcal{H}}_H,\rho \right]-\Gamma \left[\rho \right],$$

(2)

where $\rho$ is a $4\times 4$ density matrix (with density matrix elements ${\rho }_{\alpha \beta }$; $\alpha ,\beta =1,2,3,4$) describing the atomic population and coherence, and $\Gamma$ is a $4\times 4$ relaxation matrix describing the spontaneous emission and dephasing. The explicit expression of Equation (2) is presented in Appendix A. Additionally, under the condition of slowly varying envelope approximation, the probe field is described by the Maxwell equation with the form

$$\begin{array}{c}\hfill i\left(\frac{\partial }{\partial z}+\frac{1}{c}\frac{\partial }{\partial t}\right){\Omega }_p+\frac{c}{2{\omega }_p}{\nabla }_{\perp }^2{\Omega }_p+\frac{{\omega }_p}{2c}{\chi }_p{\Omega }_p=0,\end{array}$$

(3)

where ${\nabla }_{\perp }^2=\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2},{\chi }_p={\mathcal{N}}_a{({\mathbf{e}}_p·{\mathbf{p}}_{13})}^2{\rho }_{31}/({ϵ}_0\hslash {\Omega }_p)$ is the probe-field susceptibility, and c and ${ϵ}_0$ are the vacuum speed and dielectric coefficient of vacuum, respectively. Without loss of generality, the probe field is assumed propagating along z direction, i.e., ${\mathbf{k}}_p=(0,0,{\omega }_p/c)$. Meanwhile, we set the strong control and dressed laser fields along with ${\mathbf{k}}_c=(0,0,k_c),{\mathbf{k}}_d=(0,0,-k_d)$ to decrease the first-order Doppler effect. In addition, since we consider a stationary state of the system, the time derivatives in the Maxwell–Bloch Equations (2) and (3) can be neglected (i.e., $\partial /\partial t=0$), which is valid if the probe, control, and dressed fields have long enough time durations.

Since the probe field is weaker than the control and dressed fields, the standard asymptotic perturbation expansion can be applied to solve the Bloch Equation (2) by taking ${\Omega }_p$ as a small parameter [63]. The solution of Equation (2) up to the third-order is presented in Appendix B. With these solutions, we obtain the expression of the optical susceptibility of the probe field

$$\begin{array}{c}\hfill {\chi }_p={\chi }_p^{\left(1\right)}+{\chi }_{p1}^{\left(3\right)}|{\Omega }_p{|}^2+\int d^3{\mathbf{r}}^{\prime }{\chi }_{p2}^{\left(3\right)}\left({\mathbf{r}}^{\prime }-\mathbf{r}\right){\left|{\Omega }_p\left({\mathbf{r}}^{\prime }\right)\right|}^2,\end{array}$$

(4)

where ${\chi }_p^{\left(1\right)}={\mathcal{N}}_a{|{\mathbf{e}}_p·{\mathbf{p}}_{13}|}^2{a}_{31}^{\left(1\right)}/({ϵ}_0\hslash )$, ${\chi }_{p1}^{\left(3\right)}={\mathcal{N}}_a{|{\mathbf{e}}_p·{\mathbf{p}}_{13}|}^2{a}_{31,1}^{\left(3\right)}/({ϵ}_0\hslash )$, and ${\chi }_{p2}^{\left(3\right)}\left({\mathbf{r}}^{\prime }-\mathbf{r}\right)={\mathcal{N}}_a^2{|{\mathbf{e}}_p·{\mathbf{p}}_{13}|}^2{a}_{31,2}^{\left(3\right)}a{a}_{44,41}^{\left(3\right)}\left({\mathbf{r}}^{\prime }-\mathbf{r}\right)$$V_{\mathrm{vdw}}\left({\mathbf{r}}^{\prime }-\mathbf{r}\right)/({ϵ}_0\hslash )$ are, respectively, the linear, the local, and nonlocal nonlinear susceptibilities with explicit expressions of $a_{31}^{\left(1\right)}$, $a_{31,1}^{\left(3\right)}$, $a_{31,2}^{\left(3\right)}$, and $a{a}_{44,41}^{\left(3\right)}\left({\mathbf{r}}^{\prime }-\mathbf{r}\right)$, respectively, given in Appendix B. A local approximation along the z direction can be made under the condition of the spatial length of the probe pulse in the z direction is much larger than the range of atom–atom interactions. Hence, the last term on the right side of Equation (4) can be reduced as $\int d^3{\mathbf{r}}^{\prime }{\chi }_{p2}^{\left(3\right)}\left({\mathbf{r}}^{\prime }-\mathbf{r}\right)|{\Omega }_p\left({\mathbf{r}}^{\prime }\right){|}^2\approx \int \int d{x}^{\prime }d{y}^{\prime }{\tilde{\chi }}_{p2}^{\left(3\right)}(x^{\prime }-x,y^{\prime }-y){\left|{\Omega }_p(x^{\prime },y^{\prime },z)\right|}^2$ with ${\tilde{\chi }}_{p2}^{\left(3\right)}(x^{\prime }-x,y^{\prime }-y,z)=\int d{z}^{\prime }{\chi }_{p2}^{\left(3\right)}\left({\mathbf{r}}^{\prime }-\mathbf{r}\right)$. Finally, we obtain

$$\begin{array}{cc}\hfill i\frac{\partial {\Omega }_p}{\partial z}& +\frac{c}{2{\omega }_p}{\nabla }_{\perp }^2{\Omega }_p+\frac{{\omega }_p}{2c}{\chi }_p^{\left(1\right)}{\Omega }_p+\frac{{\omega }_p}{2c}[{\chi }_{p1}^{\left(3\right)}{\left|{\Omega }_p\right|}^2\hfill \\ & +\int \int d{x}^{\prime }d{y}^{\prime }{\tilde{\chi }}_{p2}^{\left(3\right)}(x^{\prime }-x,y^{\prime }-y){\left|{\Omega }_p(x^{\prime },y^{\prime },z)\right|}^2]{\Omega }_p=0.\hfill \end{array}$$

(5)

We take cold ${}^{87}$Rb atomic gas as an example. The assigned atomic levels are $|1\rangle =\left|5S_{1/2},F=1\right.〉$, $|2\rangle =\left|5S_{1/2},F=2\right.〉$, $|3\rangle =\left|5P_{3/2},F=3\right.〉$, and $|4\rangle =\left|n{S}_{1/2}\right.〉.$ For principle quantum number $n=60$, $C_6=-140\times 2\pi \mathrm{GHz}\mathsf{\mu }{\mathrm{m}}^6$ [66]. Typical system parameters are chosen as ${\mathsf{\Delta }}_1=0$, ${\mathsf{\Delta }}_2=3.17\times {10}^6{\mathrm{s}}^{-1}$, ${\mathsf{\Delta }}_3=1.53\times {10}^8{\mathrm{s}}^{-1}$, ${\mathsf{\Delta }}_4=1.32\times {10}^7{\mathrm{s}}^{-1}$, ${\Gamma }_{21}=0.1\times 2\pi \mathrm{MHz}$, ${\Gamma }_3=6.1\times 2\pi \mathrm{MHz}$, ${\Gamma }_4=0.0167\times 2\pi \mathrm{MHz}$, ${\Omega }_c=2.0\times {10}^7{\mathrm{s}}^{-1}$, ${\Omega }_d=1.5\times {10}^7{\mathrm{s}}^{-1}$, and atomic density ${\mathcal{N}}_a=1\times {10}^{11}{\mathrm{cm}}^{-3}$. With those parameters, we obtain ${\chi }_{p1}^{\left(3\right)}\approx (2.14+0.017i)\times {10}^{-17}{\mathrm{Hz}}^{-2}$, and ${\tilde{\chi }}_{p2,\max }^{\left(3\right)}\approx (1.12+0.11i)\times {10}^{-13}{\mathrm{Hz}}^{-2}$. Note that there is not only the local Kerr nonlinearity contributed by light-atom interaction, but also the nonlocal Kerr nonlinearity resulting from the atom–atom interaction. However, the later interaction plays the leading role and is three orders of magnitude larger than the local case [63]. We also note that with the above parameters we have ${\mathsf{\Delta }}_3+{\mathsf{\Delta }}_4\gg {\Omega }_d$, so the system works in the regime of Rydberg-dressed EIT [67]. We eventually obtain the dimensionless nonlocal nonlinear Schrödinger equation (NNLSE) for the probe field envelope

$$\begin{array}{cc}\hfill & i\frac{\partial u}{\partial s}=-{\nabla }_{\perp }^{2}u+V_{\mathrm{op}}(\xi ,\eta )u+\int \int d{\xi }^{\prime }d{\eta }^{\prime }\tilde{R}({\xi }^{\prime }-\xi ,{\eta }^{\prime }-\eta )|u({\xi }^{\prime },{\eta }^{\prime },s){|}^{2}u,\hfill \end{array}$$

(6)

where $u={\Omega }_p/U_0$ is the dimensionless Rabi frequency with $U_0$ the typical Rabi frequency of the weak probe field, $s=z/\left(2L_{\mathrm{diff}}\right)$ with $L_{\mathrm{diff}}={\omega }_p{R}_{\perp }^2/c$ the typical diffraction length, $(\xi ,\eta )=(x,y)/R_{\perp },[({\xi }^{\prime },{\eta }^{\prime })=(x^{\prime },y^{\prime })/R_{\perp }]$ are dimensionless spatial coordinate with $R_{\perp }$ the typical radius of the probe beam, and $V_{\mathrm{op}}=-{\omega }_p^2{R}_{\perp }^2{\chi }_p^{\left(1\right)}/c^2$ the optical potential. Here, $\tilde{R}({\xi }^{\prime }-\xi ,{\eta }^{\prime }-\eta )=L_{\mathrm{diff}}{\left|U_0\right|}^2{R}_{\perp }^2{\omega }_p{\tilde{\chi }}_{p2}^{\left(3\right)}/c\equiv \alpha a_{44,41}^{\left(3\right)}[({\xi }^{\prime }-\xi )R_{\perp },({\eta }^{\prime }-\eta )R_{\perp }]V_{\mathrm{vdw}}[({\xi }^{\prime }-\xi )R_{\perp },({\eta }^{\prime }-\eta )R_{\perp }]$ with $\alpha =|{\mathbf{e}}_p·{\mathbf{p}}_{13}{|}^2{a}_{31,2}^{\left(3\right)}{R}_{\perp }^4{\omega }_p^2{U}_0^2{\mathcal{N}}_a^2/({ϵ}_0\hslash c^2)$. Since the local nonlinearity is much smaller than the nonlocal nonlinearity due to atom–atom interaction, the local term has been canceled in Equation (6).

Figure 1c shows a cross-sectional diagram ($\eta =0$) of the response function $\tilde{R}\left(\xi \right)$ as a function of $\xi$. We also give the response function $\mathcal{R}\left({\beta }_1\right)$ in momentum space as a function of ${\beta }_1$ with ${\beta }_2=0$ in Figure 1d. We can see that the nonlocal response function approaches saturation at small distance, and the response function in momentum space changes its signs at some interval, which is the key to the appearance of MI (discussed later). The nonlocality degree can be characterized by the parameter $\sigma \equiv R_b/R_{\perp }$, where $R_b$ is the radius of the blockade sphere, which is estimated by $R_b={|C_6/{\delta }_{\mathrm{EIT}}|}^{1/6}$ with ${\delta }_{\mathrm{EIT}}$ the linewidth of EIT transition spectrum, calculated by ${\delta }_{\mathrm{EIT}}={\left|{\Omega }_c\right|}^3/{\gamma }_{31}$ for ${\mathsf{\Delta }}_3=0$ and ${\delta }_{\mathrm{EIT}}={\left|{\Omega }_c\right|}^2/{\mathsf{\Delta }}_3$ for ${\mathsf{\Delta }}_3\gg {\gamma }_{31}$. Typically, $R_b\approx 8.34\mathsf{\mu }\mathrm{m}$ in our system with the parameters above.

3. Realization of Non-Hermitian Optical Potentials

3.1. Constant-Intensity Wave Solution

Our main goal is to investigate the MI and related optical pattern formations in the system with non-Hermitian optical potentials starting from the constant-intensity background. To this end, we first consider the physical realization of various non-Hermitian optical potentials and give a corresponding constant-intensity wave (CIW) solution. It is easy to show that the NNLS Equation (6) admits the 2D CIW solution

$$\begin{array}{c}\hfill u(\xi ,\eta ,s)=Aexp\left[i{\int }_l(d\xi W_{\mathrm{I}}+d\eta W_{\mathrm{II}})-i\mu s\right],\end{array}$$

(7)

if the optical potential $V_{\mathrm{op}}(\xi ,\eta )$ satisfies condition [10]

$$\begin{array}{cc}\hfill & V_{\mathrm{op}}(\xi ,\eta )=-|W_{\mathrm{I}}{|}^2-{\left|W_{\mathrm{II}}\right|}^2+i\left(\frac{\partial W_{\mathrm{I}}}{\partial \xi }+\frac{\partial W_{\mathrm{II}}}{\partial \eta }\right),\mathrm{and}\frac{\partial W_{\mathrm{II}}}{\partial \xi }=\frac{\partial W_{\mathrm{I}}}{\partial \eta }.\hfill \end{array}$$

(8)

where $\mu =\int \int d\xi d\eta \tilde{R}(\xi ,\eta )$, $W_{\mathrm{I}}$ and $W_{\mathrm{II}}$ are real functions of $\xi$ and $\eta$, A is an arbitrary real, positive constant, and l is any smooth open curve connecting an arbitrary point $({\xi }_0,{\eta }_0)$ to any different point $(\xi ,\eta )$. The detailed derivation of the constant intensity wave solution (7) is given in Appendix C. For simplicity and without loss of generality, we consider in the following only three special optical potentials: (i) the $\mathcal{PT}$-symmetric potential determined by $W_{\mathrm{I}}=acos\left(\xi \right)sin\left(\eta \right),W_{\mathrm{II}}=acos\left(\eta \right)sin\left(\xi \right)$; and (ii) the non-$\mathcal{PT}$-symmetric potential determined by $W_{\mathrm{I}}=b+a\eta +F_0\xi exp[-{\xi }^2/\left(2{\sigma }_1^2\right)]$ and $W_{\mathrm{II}}=b+a\xi +G_0\eta exp[-{\eta }^2/\left(2{\sigma }_2^2\right)]$; (iii) the ring-shaped optical potential (see below), where a, b, $F_0$, and $G_0$ are the arbitrary real constants.

3.2. Realization of $\mathcal{PT}$-Symmetric Optical Potential

In the special case where $W_{\mathrm{I}}$ and $W_{\mathrm{II}}$ are even functions, the actual optical potential $V_{\mathrm{op}}(\xi ,\eta )$ turns out to be $\mathcal{PT}$-symmetric. Thus we will consider how to realize a linear 2D $\mathcal{PT}$-symmetric optical potential in our Rydberg system (i.e., satisfying the condition of $V_r(-\xi ,\eta )=V_r(\xi ,\eta ),V_r(\xi ,-\eta )=V_r(\xi ,\eta )$ and $V_i(-\xi ,\eta )=-V_i(\xi ,\eta ),V_i(\xi ,-\eta )=-V_i(\xi ,\eta )$). Here the subscripts r and i represent the real and imaginary parts of the potential, respectively. To the aim, we assume that our target optical potential is $\mathcal{PT}$-symmetric with the form

$$\begin{array}{cc}\hfill V_{\mathrm{op},1}(\xi ,\eta )=& -V_0^2[{cos}^2\left(\xi \right){sin}^2\left(\eta \right)+{cos}^2\left(\eta \right){sin}^2\left(\xi \right)]-2V_{0}isin\left(\xi \right)sin\left(\eta \right),\hfill \end{array}$$

(9)

where $V_0^2$ and $2V_0$ are constants which characterize the amplitudes of the real and imaginary parts of the periodic optical potential, respectively. The potential (9) can be obtained by taking $W_{\mathrm{I}}=V_{0}cos\left(\xi \right)sin\left(\eta \right)$ and $W_{\mathrm{II}}=V_{0}cos\left(\eta \right)sin\left(\xi \right)$. The eigenvalue problem of the linearized Equation (6) reads $\widehat{L}v=\lambda v$, here $\widehat{L}\equiv {\partial }^2/\partial {\xi }^2+{\partial }^2/\partial {\eta }^2-V_{\mathrm{op},1}$ and the $\lambda$ is defined through the relation $u(\xi ,\eta ,s)=v(\xi ,\eta )exp\left(i\lambda s\right)$. The eigen spectrum is real (complex) for $V_0<1.25$ ($V_0>1.25$), working in unbroken (broken) $\mathcal{PT}$ phase with phase transition point $V_0=1.25$. Hereafter, we consider only the case of the unbroken $\mathcal{PT}$ phase and $V_0=1$.

To realize the such potential in the Rydberg-dressed EIT, we follow the standard method given in Refs. [68,69,70]. Firstly, we note that the optical potential with $\mathcal{PT}$ symmetry must satisfy the condition $\mathrm{Im}\left[V_{\mathrm{op}}\right](\xi =0,\eta =0)=0$ at the origin to balance gain and loss in the whole space. Then, we assume the half Rabi frequencies of the control and dressed fields be space-dependent, i.e., ${\Omega }_c={\Omega }_c(\xi ,\eta ),{\Omega }_d={\Omega }_d(\xi ,\eta )$. Hence the optical potential is a function of ${\Omega }_c$ and ${\Omega }_d$. Finally, by Taylor expanding $V_{\mathrm{op}}$ around ${\Omega }_c={\Omega }_{c0}(\xi =0,\eta =0)$ and ${\Omega }_d={\Omega }_{d0}(\xi =0,\eta =0)$ to the first-order and separating the real and imaginary parts, we obtain the following equations

$$\begin{array}{cc}\hfill & R_c{\Omega }_c+R_d{\Omega }_d=V_{\mathrm{tr}}(\xi ,\eta )-V_{\mathrm{tr}}({\Omega }_{c0},{\Omega }_{d0})+R_c{\Omega }_{c0}+R_a{\Omega }_{d0},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & I_c{\Omega }_c+I_d{\Omega }_d=V_{\mathrm{ti}}(\xi ,\eta )+I_c{\Omega }_{c0}+I_a{\Omega }_{d0},\hfill \end{array}$$

(11)

where $R_d\equiv \frac{\partial V_r}{\partial {\Omega }_d}$, $R_c\equiv \frac{\partial V_r}{\partial {\Omega }_c}$, $I_d\equiv \frac{\partial V_i}{\partial {\Omega }_d}$, and $I_c\equiv \frac{\partial V_i}{\partial {\Omega }_c}$ are the derivations at ${\Omega }_c={\Omega }_{c0},{\Omega }_d={\Omega }_{d0}$, and $V_{\mathrm{tr}},V_{\mathrm{ti}}$ are the real and imaginary of the target optical potential $V_{\mathrm{op},1}$, respectively. Finally, solving the coupled equations above, we can choose the control and dressed fields with spatial distribution, i.e.,

$$\begin{array}{cc}\hfill {\Omega }_c& (\xi ,\eta )\approx {\Omega }_{c0}\left\{0.993-0.022sin\left(\xi \right)sin\left(\eta \right)+0.0292cos\left(2\xi \right)cos\left(2\eta \right)\right\},\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\Omega }_d& (\xi ,\eta )\approx {\Omega }_{d0}\left\{0.989-0.0159cos\left(2\xi \right)sin\left(2\eta \right)+0.186sin\left(\xi \right)sin\left(\eta \right)\right\},\hfill \end{array}$$

(13)

where ${\Omega }_{c0}=2.0\times {10}^7{\mathrm{s}}^{-1}$ and ${\Omega }_{d0}=1.5\times {10}^7{\mathrm{s}}^{-1}$. The distributions of the obtained control and assistant fields are plotted in Figure 2a,b with $V_0=1$. Note that the periodic spatial modulation of the control and dressed fields can be realized experimentally by using a high-resolution space–light modulator with a pixel size smaller than the probe-beam radius $R_{\perp }$.

3.3. Realization of Non-$\mathcal{PT}$-Symmetric Optical Potential

The non-$\mathcal{PT}$-symmetric optical confining potential is of the form

$$\begin{array}{cc}\hfill V_{\mathrm{op},2}& (\xi ,\eta )=-W_{\mathrm{I}}^2-W_{\mathrm{II}}^2+i\left\{F_{0}exp\left(-\frac{{\xi }^2}{2{\sigma }_1^2}\right)\left(1-\frac{{\xi }^2}{{\sigma }_1^2}\right)+G_{0}exp\left(-\frac{{\eta }^2}{2{\sigma }_2^2}\right)\left(1-\frac{{\eta }^2}{{\sigma }_2^2}\right)\right\}.\hfill \end{array}$$

(14)

where $W_{\mathrm{I}}=a\eta +F_0\xi exp(-{\xi }^2/2{\sigma }_1^2)$, and $W_{\mathrm{II}}=a\xi +G_0\xi exp(-{\eta }^2/2{\sigma }_2^2)$. To achieve this potential, we design the laser fields with the form

$$\begin{array}{cc}\hfill {\Omega }_c(\xi ,\eta )& \approx {\Omega }_{c0}\left(1.004-0.022s_1+0.1s_2\right),\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill {\Omega }_d(\xi ,\eta )& \approx {\Omega }_{d0}\left(1.006-0.032s_1+0.65s_2\right),\hfill \end{array}$$

(16)

where $a=0.01,F_0=G_0=7,{\sigma }_1={\sigma }_2=5$, $s_1={[0.01\eta +7\xi exp(-{\xi }^2/50)]}^2+{[0.01\xi +7\eta exp(-{\eta }^2/50)]}^2$ and $s_2=(1-{\xi }^2/25)exp(-{\xi }^2/50)+(1-{\eta }^2/25)exp(-{\eta }^2/50)$. The distribution of laser fields are shown in Figure 2c,d.

4. MI and Pattern Formation for the $\mathcal{PT}$-Symmetric Periodic Potential

4.1. MI for the $\mathcal{PT}$-Symmetric Periodic Potential

The $\mathcal{PT}$-symmetric potential has the form $V_{\mathrm{op},1}(\xi ,\eta )=V_r(\xi ,\eta )+i{V}_i(\xi ,\eta )$ which is given in Equation (9). In Ref. [10], the authors consider the nonlinear Schrödinger equation with a general, non-Hermitian potential and a local Kerr nonlinearity and prove that the system allows constant-intensity wave solution. According to analogue analysis and simple mathematical deduction, a constant intensity probe field solution appears in our nonlocal Rydberg system described by Equation (6), i.e., $u(\xi ,\eta ,s)=Aexp\left[i{V}_{0}sin\left(\xi \right)sin\left(\eta \right)\right]exp\left(i\mu s\right)$, with $\mu =-A^2\int \int \tilde{R}(\xi ,\eta )d\xi d\eta$ the propagating constant and A a constant.

To study the stability of constant intensity solution in our Rydberg system, a simple linear stability analysis of constant intensity wave by the perturbation was implemented, i.e., $u=[A+ϵ\tilde{F}(\xi ,\eta )e^{i\overrightarrow{\zeta }·\overrightarrow{\beta }+i\lambda s}+ϵ{\tilde{G}}^{\ast }(\xi ,\eta )e^{-i\overrightarrow{\zeta }·\overrightarrow{\beta }-i{\lambda }^{\ast }s}]exp[i{V}_{0}sin\left(\xi \right)sin\left(\eta \right)]exp\left(i\mu s\right)$, where $\overrightarrow{\zeta }=(\xi ,\eta )$ is the dimensionless coordinates, and $\overrightarrow{\beta }=({\beta }_1,{\beta }_2)$ is the dimensionless wave vector corresponding to $(k_x,k_y)$. Inserting this equation into Equation (6) and linearizing with perturbation amplitudes $\tilde{F},\tilde{G}$, we have

$$\begin{array}{cc}\hfill ({\nabla }_{\perp }^2& -{\beta }^2+{\mathcal{L}}_{11})\tilde{F}+A^2\widehat{\mathcal{N}}=\lambda \tilde{F},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill ({\beta }^2& -{\nabla }_{\perp }^2+{\mathcal{L}}_{22})\tilde{G}-A^2\widehat{\mathcal{N}}=\lambda \tilde{G},\hfill \end{array}$$

(18)

where $\widehat{\mathcal{N}}=-\int \int d{\xi }^{\prime }d{\eta }^{\prime }\tilde{R}({\xi }^{\prime }-\xi ,{\eta }^{\prime }-\eta )[\tilde{F}({\xi }^{\prime },{\eta }^{\prime })+\tilde{G}({\xi }^{\prime },{\eta }^{\prime })]e^{i\overrightarrow{\beta }·(\overrightarrow{{\zeta }^{\prime }}-\overrightarrow{\zeta })}$, ${\beta }^2={\beta }_1^2+{\beta }_2^2$, ${\mathcal{L}}_{11}=2i{L}_1+2i{V}_0{L}_2-2V_0{L}_3$, ${\mathcal{L}}_{22}=-2i{L}_1+2i{V}_0{L}_2-2V_0{L}_3$, $L_1={\beta }_1\frac{\partial }{\partial \xi }+{\beta }_2\frac{\partial }{\partial \eta }$, $L_2=sin\left(\eta \right)cos\left(\xi \right)$$\frac{\partial }{\partial \xi }+sin\left(\xi \right)cos\left(\eta \right)\frac{\partial }{\partial \eta }$, and $L_3={\beta }_{1}sin\left(\eta \right)cos\left(\xi \right)+{\beta }_{2}sin\left(\xi \right)cos\left(\eta \right)$.

We use the Fourier collocation method [10,71] to study the MI of the CIW solution according to Equations (17) and (18). The eigenvalue equation tells us that if there exists an eigenvalue whose imaginary part $\mathrm{Im}\left[\lambda \right]$ is larger than zero, the solution is linearly unstable. Figure 3a–c show the growth rate $\lambda$ as functions of ${\beta }_1$ and ${\beta }_2$ for different nonlocal nonlinearity $\alpha =2,5$, and 10, respectively. From the gain spectrum, we know there exists an instability interval in which the system undergoes from CIW to some other states when we increase the nonlinearity. Figure 3d is the cross section of the eigenvalue spectrum at ${\beta }_2=-0.5$ as function of ${\beta }_1$ corresponding to $\alpha =2,3,5,10$, and 20 (from lines a to e), which further illustrates the instability interval increases with $\alpha$ increasing. The results show that nonlocal nonlinearity is the key to MI.

According to the MI analysis of combining with $\mathcal{PT}$ symmetric optical potential, we know it is the nonlocal Rydberg interaction that makes our system distinctive. It is generally known that there is always MI in self-focusing nonlinearity both in local and nonlocal nonlinearity [13,15,16], but in the nonlocal self-defocusing case, the MI is strongly dependent on the detailed form of the response function, which demands the response to have some sign change in the momentum space. Fortunately, the response function in our Rydberg system is the soft-core potential shown in Figure 1c, and the momentum space of the response shown in Figure 1d changes the sign.

4.2. Pattern Formation for the $\mathcal{PT}$-Symmetric Periodic Potential

Based on the MI, SSB appearing in our Rydberg system is obvious. We solve Equation (6) with different $\alpha$ and the nonlocality degree $\sigma$ of the nonlocal Kerr nonlinearity under the CIW $u_0(\xi ,\eta ,s=0)=Aexp[i{V}_{0}sin\xi sin\eta ]$ by imaginary time evolution method [71]. We study the MI with its consequences systematically, and eventually, we obtain the different solutions from CIW solution to tetragonal and hexagonal lattice solutions.

We first illustrate the phase diagram as functions of nonlocality degree $\sigma$ and nonlinear interaction strength $\alpha$, with the range of constant intensity probe field amplitude $\alpha$ from 3 to 12 and nonlocality degree $\sigma$ from $0.3$ to $1.2$ in Figure 4a, which represents the emergence of the different phases from the constant intensity wave (labeled by “CIW”) solution (green region) to the tetragonal lattice (purple region) and hexagonal lattice (red region), respectively. We then track the system energy by minimizing $E=\int |{\nabla }_{\perp }u\left(\overrightarrow{\zeta }\right){|}^2{d}^2\overrightarrow{\zeta }-1/2\int \int \tilde{R}\left({\overrightarrow{\zeta }}^{\prime }-\overrightarrow{\zeta }\right)|u\left(\overrightarrow{\zeta }\right){|}^2{\left|u\left({\overrightarrow{\zeta }}^{\prime }\right)\right|}^2{d}^2{\overrightarrow{\zeta }}^{\prime }{d}^2\overrightarrow{\zeta }$$+\int \int V_{\mathrm{op}}{\left|u\left(\overrightarrow{\zeta }\right)\right|}^2{d}^2\overrightarrow{\zeta }$. The inset in Figure 4b shows the energy change along $s=z/2L_{\mathrm{diff}}$, and its energy becomes lower, eventually remains unchanged, and comes to the ground state. To view the tendency of the energy change, we give the partially enlarged drawing inside the inset. We also show the difference in Hamiltonian energy $E_{\mathrm{gs}}-E_{\mathrm{CIW}}$ (blue line in left y-axis) and propagation constant ${\mu }_{\mathrm{gs}}-{\mu }_{\mathrm{CIW}}$ (red line in right y-axis) as functions of $\alpha$ from 2 to 12 with $\sigma =0.6$ in Figure 4c, which is corresponding to the vertical dotted line crossed from the point “P”, i.e., $\mathrm{P}(0.6,5.3)$ in Figure 4a with $\sigma =0.6$, and $\alpha$ from 3 to 12. We find that the difference in Hamiltonian energy or propagation constant has a phase transition at point “P”.

Furthermore, according to the phase diagram, we choose two groups of parameters to solve the equation and obtain the two types of pattern formations. One is the tetragonal lattice which is a function of dimensionless coordinate $\xi$ and $\eta$ with nonlocality degree $\sigma =1.1$, and $\alpha =4$ in Figure 4d; the other one is the hexagonal lattice which is in Figure 4e with $\alpha =6$ and $\sigma =0.7$. Different patterns can be understood as follows. On one hand, intuitively, by detailed analysis of the original NNLSE Equation (6), we find this equation includes diffraction, non-Hermitian potential, and nonlocal nonlinearity, which can interact with each other. While the diffraction effect tends to make the wave function extension much flatter, the non-Hermitian potential can modify the diffraction property of the system (i.e., in some points of the Brillouin region, the diffraction can change its sign). On the other hand, the optical potential and Van der Waals interaction simultaneously work on the system. When one interaction is stronger than another interaction, the former plays the dominant role.

5. MI and Pattern Formation for Non-$\mathcal{PT}$-Symmetric Confining Potential

5.1. MI for Non-$\mathcal{PT}$-Symmetric Optical Confining Potential

We then turn our attention to studying the MI of the CIW under non-$\mathcal{PT}$-symmetric confining (Hermite-Gaussian) optical potential, which is given in Equation (14). The solution is given by $u(\xi ,\eta ,s)=Aexp\left[i\left\{s\mu +a\xi \eta +F_0{\sigma }_1^2\left[1-exp(-{\xi }^2/2{\sigma }_1^2)\right]+G_0{\sigma }_2^2\left[1-exp(-{\eta }^2/2{\sigma }_2^2)\right]\right\}\right]$. Using a similar method to that used in the previous section and adding small amplitude perturbation $\tilde{F}$ and $\tilde{G}$ ($|\tilde{F}|\ll A$ and $|\tilde{G}|\ll A$) into CIW, we have

$$\begin{array}{cc}\hfill & [{\nabla }^2-{\beta }^2+{\mathcal{L}}_{11}-V_{\mathrm{op},2}-\mu ]\tilde{F}+A^2\widehat{\mathcal{N}}=\lambda \tilde{F},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & [-{\nabla }^2+{\beta }^2+{\mathcal{L}}_{22}+V_{\mathrm{op},2}+\mu ]\tilde{G}-A^2\widehat{\mathcal{N}}=\lambda \tilde{G},\hfill \end{array}$$

(20)

where ${\mathcal{L}}_{11}=2i{L}_1-2L_2+2i{L}_3$, ${\mathcal{L}}_{22}=-2i{L}_1-2L_2+2i{L}_3$, $L_1={\beta }_1\frac{\partial }{\partial \xi }+{\beta }_2\frac{\partial }{\partial \eta }$, $L_2={\beta }_1{W}_{\mathrm{I}}+{\beta }_2{W}_{\mathrm{II}}$, and $L_3=W_{\mathrm{I}}\frac{\partial }{\partial \xi }+W_{\mathrm{II}}\frac{\partial }{\partial \eta }$. Solving the eigenvalue Equations (19) and (20) numerically, we obtain the gain spectrum in Figure 5a, which shows the gain interval lying in the vicinity of $|{\beta }_1|=|{\beta }_2|$ in the first Brillouin region. We also give the cross-section curves of the eigenvalue spectrum in Figure 5b as a function of ${\beta }_1$ corresponding to ${\beta }_2=0$ (red line), ${\beta }_2=0.25$ (green dotted line), and ${\beta }_2=0.5$ (blue dashed-dotted line), respectively.

5.2. Pattern Formation for Non-$\mathcal{PT}$-Symmetric Confining Potential

We start with CIW $u(\xi ,\eta )=Aexp\left\{i\left\{0.01\xi \eta +175\left[2-exp(-{\xi }^2/50)-exp(-{\eta }^2/50)\right]\right\}\right\}$ perturbed by small noise and simulate Equation (6) numerically. First, the phase diagram is given in Figure 5c as functions of nonlocality degree $\sigma$ (range from 0.3 to 1.2) and nonlinear interaction strength $\alpha$ (from 2 to 10). It shows that the CIW solution can experience MI in other states by manipulating the nonlinearity or nonlocality degree. The white dashed line is the boundary of the stable pattern and constant intensity wave. We also show the pattern with $\sigma =0.9$ and $\alpha =10$ in Figure 5d, in which we also plot the contours of the optical potential with a dotted line.

The pattern can be explained by the co-interaction of diffraction effect, potential, and nonlocal self-defocused interaction. When both the nonlocality degree and nonlinearity are small, the nonlocal nonlinearity is reduced to local nonlinearity, and the diffraction effect tends to make the wave function extension much flatter; eventually, the CIW is still stable. While we increase nonlocality degree and nonlinearity, the nonlocal nonlinearity plays the leading role, which tends to make the wave function extensive, while the potential traps the wave function. Eventually, the wave will become a hexagonal pattern when reaching balance between potential and nonlinearity.

5.3. Pattern Formation for the Ring-Shaped Non-$\mathcal{PT}$-Symmetric Potential

We also consider another potential $V_r=V_0[exp(-{\zeta }^8/a^8)-exp(-{\zeta }^8/b^8)]$ and $V_i=-V_0^2{\left[a^2\Gamma \left(1/4,{\zeta }^8/a^8\right)-b^2\Gamma \left(1/4,{\zeta }^8/b^8\right)\right]}^2/\left(64{\zeta }^2\right)$. Here, $V_0$ can represent the amplitude of the potential, $\zeta =\sqrt{{\xi }^2+{\eta }^2}$, $a,b$ are the width of the super-Gaussian function, and $\Gamma (·,·)$ is the incomplete gamma function. The CIW solution gives the form $u=Aexp\left[i\theta \right(\xi ,\eta \left)\right]$ with A a constant, and the detailed form of $\theta (\xi ,\eta )$ is given in Appendix C.

The pattern formations obtained by numerical simulation in Equation (6) are shown in Figure 6. When the width of the super-Gaussian potential is $a=3,b=4$, the singer layer ring-shaped structure pattern with the smaller radius is obtained in Figure 6a. When we increase the width, i.e., $a=6,b=8$ and $a=4,b=7$, the larger ring-shaped pattern is obtained in Figure 6b,c. Furthermore, the two or more two layers of a ring-shaped pattern with a hexagonal structure have been obtained with $a=2$ and much larger $b=9$ which is shown in Figure 6d.

According to the simulation with the typical parameters, we find that by manipulating the parameters a or b, on the one hand, we can obtain the ring-like structure patterns with different radii. On the other hand, we can change the layers of the ring. That is, the larger a and b correspond to the pattern with the larger radius, but the larger $|b-a|$ results in more layers of the pattern.

6. Summary

The predictions of the optical patterns presented above may be observed experimentally by using a cold Rydberg atomic gas with the inverted Y-shaped level configuration. Since these lattice patterns form in the transverse x-y plane and the typical lattice separation is around 8 $\mathsf{\mu }$m, one should prepare a Rydberg gas with transverse size of several tens of micrometers (e.g., 30 µm) and with longitudinal size of several millimeters (e.g., 2 mm), which can be realized by current experimental techniques [72]. To realize the predicted optical structures, one can inject a continuous-wave probe field; under the condition of the EIT and by setting suitable system parameters, the probe field will undergo MI and then be transformed into the optical patterns when it propagates along the z-direction for several millimeters. In the present work, we have limited our considerations under the ultra-cold Rydberg atomic gas. It is possible that the system supports optical structures in thermal atomic gases. There are many inherent broadening effects, such as Doppler broadening, transition broadening, collision broadening, etc., which should be considered. We think it may have a richer self-organization structure, which is an interesting topic deserving to be explored further.

In conclusion, we have proposed a realistic scheme to realize non-Hermitian optical potentials with or without parity-time symmetry and hence find the pattern formations through the self-organization of laser light in a cold Rydberg atomic gas. We have shown that through nonlocal Kerr nonlinearity contributed from the atom–atom interaction, a CIW probe laser field can undergo an MI, and induce SSB and hence result in various self-organized optical patterns (including hexagonal and tetragonal lattices, etc.), which can be actively manipulated by changing the Kerr nonlinearity, its nonlocality degree, and the potential depth of the non-Hermitian optical potentials. The results reported here open a route for exploring the SSB and self-organization of laser light for non-Hermitian nonlinear optics based on cold Rydberg gases.