Higher-order exceptional points and asymmetric soliton scattering with triple-well \mathcal{P}\mathcal{T}-symmetry potential in an atomic gas

Abstract

We propose a physical scheme to realize a parity-time (\mathcal{P}\mathcal{T})-symmetric triple-well potential in a coherent atomic gas and investigate the emergence of higher-order exceptional points (EPs). Through careful parameter tuning, we successfully identify a third-order EP (EP3), observing its characteristic cube-root response to perturbations—a key signature distinguishing it from second-order EPs (EP2). Notably, EP3 exhibits significantly enhanced sensitivity, underscoring its potential for advancing quantum information technologies. Furthermore, we show that the imaginary part of the nonlinear \mathcal{P}\mathcal{T} potential critically governs \mathcal{P}\mathcal{T} phase transitions and dynamically reshapes the phase diagram. Our system also supports stable optical solitons, whose properties can be precisely controlled by tuning the interplay between linear and nonlinear \mathcal{P}\mathcal{T} potentials. When these potentials act as defects, soliton scattering displays strongly asymmetric behavior, directly adjustable via the potentials’ imaginary components. These findings open new avenues for applications in optical switch, optical sensing, and high-precision information transmission.

Graphical abstract

Keywords

non-Hermitian / \mathcal{P}\mathcal{T}-symmetry / phase diagram / exceptional point / soliton scattering

1 Introduction

Exceptional points (EPs), non-Hermitian degeneracies where eigenvalues and eigenvectors coalesce, have recently emerged as a powerful paradigm for controlling the response of open physical systems [1, 2]. Second-order EPs (EP2) are the most studied because of their abundance, which requires only the tuning of one parameter [3–6]. In contrast, higher-order EPs (greater than second order) generically require more fine-tuning parameters. Recent studies have demonstrated that third-order EPs (EP3) generically require only two real tuning parameters in the presence of a parity−time (\mathcal{P}\mathcal{T}) symmetry [7–13]. Furthermore, higher-order EPs can amplify the effect of perturbations, enhancing the sensitivity [14–20]. These advances may pave the way towards a new class of ultrasensitive sensing technologies. With the advance of research, it has been found that \mathcal{P}\mathcal{T} symmetric optics, may have many practical applications. These include unidirectional light propagation [21–27], coherent perfect absorbers [28–31], giant light amplification [32], novel lasers [33, 34], precision metrology [14, 17], quantum computation [35], among others. Beyond these phenomena, nonlinear waves, particularly optical solitons, have also been extensively studied in \mathcal{P}\mathcal{T}-symmetric systems [1, 36–39], expanding their utility in photonic technologies.

More recently, it has been shown that coherently atomic gases interacting with laser fields may provide a fertile platform for studying optical \mathcal{P}\mathcal{T} symmetry [38, 40–45]. This system provides exceptional tunability of multiple parameters, making it particularly suitable for exploring higher-order EPs. Notably, the optical refractive index in such systems can be actively manipulated, while the Kerr nonlinearity in such systems can be significantly enhanced through resonant nonlinear optical susceptibilities. Thus, the non-Hermitian optics with coherently atomic gases opens a new route for the investigation of \mathcal{P}\mathcal{T}-symmetric linear and nonlinear quantum mechanics.

Motivated by higher-order EPs, there is growing interest in developing \mathcal{P}\mathcal{T}-symmetric systems that not only have a gain-loss linear potential but also an additional neutral one [8, 14]. In this work, we propose a realistic scheme for physically realizing gain-neutral-loss optical potentials with \mathcal{P}\mathcal{T} symmetry in a coherent atomic gas, achieved through a spatial modulation of control laser field. We find that it is quite easy to achieve EP2 and EP3 by adjusting the imaginary part of the \mathcal{P}\mathcal{T} potential in this system. We shall demonstrate that the spatially dependent imaginary part of the nonlinear \mathcal{P}\mathcal{T} potential plays a crucial role in determining the \mathcal{P}\mathcal{T} phase transition and the actively modifying of \mathcal{P}\mathcal{T} phase diagram. Furthermore, we shall show that the system supports stable optical solitons, which can be managed via tuning the combined \mathcal{P}\mathcal{T} potentials. Additionally, the scattering of the optical solitons by the defect exhibits obvious asymmetric behavior by taking such combined \mathcal{P}\mathcal{T} potentials as a defect. It primarily controlled by the spatially dependent imaginary parts of the combined linear and nonlinear potentials.

Before proceeding, we emphasize our work differs significantly from previous studies [46–58] in several key aspects. First, we present a practical scheme for realizing the \mathcal{P}\mathcal{T} potentials and investigate the higher-order EPs, providing direct guidance for realistic implementations. Second, our analysis reveals that the property of the \mathcal{P}\mathcal{T} phase transition has a heavy dependence on the spatially dependent imaginary part of the nonlinear \mathcal{P}\mathcal{T} potential, offering new avenues for active control. Third, we conduct a detailed study of soliton scattering, uncovering striking nonreciprocal and asymmetric behavior, a phenomenon with profound implications for non-Hermitian photonics.

2 Physical model

We start with considering a cold, dilute N-type four-level {}^{85}Rb atomic gas, which consists of two hyperfine ground states |1⟩=|5^2{S}_{1/2},F=2⟩, |2⟩=|5^2{S}_{1/2},F=3⟩, and two excited states |3⟩=|5^2{P}_{1/2}⟩ and |4⟩=|5^2{P}_{3/2}⟩ [42]. A weak signal laser field with half Rabi frequency {\mathrm{\Omega }}_s couples the transition |1⟩\leftrightarrow |3⟩, a strong control and a strong pump laser fields with half Rabi frequencies {\mathrm{\Omega }}_c and {\mathrm{\Omega }}_p couple the transitions |2⟩\leftrightarrow |3⟩ and |1⟩\leftrightarrow |4⟩, respectively, as shown in Fig. 1(a) [41, 45]. The detuning {\mathrm{\Delta }}_{\alpha } (\alpha =2,3,4) represents the frequency difference between laser field and the corresponding atomic transition. {\mathrm{\Gamma }}_{\alpha \beta } denote spontaneous emission decay rates from state |\beta ⟩ to |\alpha ⟩ (). Furthermore, all laser fields are assumed to propagate along z direction through the atomic ensemble to suppress the Doppler effect. In fact, the atomic excitation scheme [Fig. 1(a)] can be taken as a system of three-level \mathrm{\Lambda }-type electromagnetically induced transparency (EIT) with an optical pumping that provides a gain to the signal field [59].

Under electric-dipole and rotating-wave approximations, the Hamiltonian of the system in the interaction picture is given by

\begin{array}{rl}{\hat{H}}_{\mathrm{i}\mathrm{n}\mathrm{t}}=& -\hslash \sum _{\alpha =2}^4{\mathrm{\Delta }}_{\alpha }|\alpha ⟩⟨\alpha |-\hslash [{\mathrm{\Omega }}_p|4⟩⟨1|+{\mathrm{\Omega }}_s|3⟩⟨1|\\ & +{\mathrm{\Omega }}_c|3⟩⟨2|+H.c.].\end{array}

Here {\mathrm{\Delta }}_2={\omega }_s-{\omega }_c-(E_2-E_1)/\hslash is two-photon detuning, {\mathrm{\Delta }}_3={\omega }_c-(E_3-E_1)/\hslash and {\mathrm{\Delta }}_4={\omega }_p-(E_4-E_1)/\hslash are one-photon detunings (with E_{\alpha } the eigenenergy of the atomic state |\alpha ⟩), {\mathrm{\Omega }}_s=({\mathit{e}}_s\cdot {\mathit{p}}_{13}){\mathcal{E}}_s/(2\hslash ), {\mathrm{\Omega }}_c=({\mathit{e}}_c\cdot {\mathit{p}}_{23}){\mathcal{E}}_c/(2\hslash ), and {\mathrm{\Omega }}_p=({\mathit{e}}_p\cdot {\mathit{p}}_{14}){\mathcal{E}}_p/(2\hslash ) are, respectively, half Rabi frequencies of the signal, control, and pump fields, where {\mathit{e}}_l and {\mathcal{E}}_l (l=s,c,p) are respectively polarization unit vectors and envelopes of the lasers (signal, control, and pump) fields, and {\mathit{p}}_{\alpha \beta } the electric-dipole matrix elements associated with the transition |\alpha ⟩\leftrightarrow |\beta ⟩.

The dynamics of atoms is governed by the optical Bloch equation

(1){\displaystyle \frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}}=-{\displaystyle \frac{\mathrm{i}}{\hslash }}\left[{\hat{H}}_{\mathrm{i}\mathrm{n}\mathrm{t}},\rho \right]+\mathcal{L}[\rho ],

where \rho is the density matrix (with density matrix elements {\rho }_{\alpha \beta }; \alpha ,\beta =1,2,3,4), describing the atomic population and coherence; and \mathcal{L} is the relaxation matrix, describing the spontaneous emission and dephasing. The explicit expression of Eq. (1) is presented in Appendix A.

Propagation of the probe field is governed by the Maxwell equation:

(2)\mathrm{i}\left({\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }z}}+{\displaystyle \frac{1}{c}}{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}}\right){\mathrm{\Omega }}_s+{\displaystyle \frac{c}{2{\omega }_s}}{\mathrm{\nabla }}_{\perp }^2{\mathrm{\Omega }}_s+{\displaystyle \frac{{\omega }_s}{2c}}{\chi }_s{\mathrm{\Omega }}_s=0,

where {\mathrm{\nabla }}_{\perp }^2=\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{y}^2}, {\chi }_s={\mathcal{N}}_a|{\mathit{e}}_s\cdot {\mathit{p}}_{13}{|}^2{\rho }_{31}/({\epsilon }_0\hslash {\mathrm{\Omega }}_s) is the optical susceptibility of the signal field, with {\mathcal{N}}_a the atomic density and {\epsilon }_0 the vacuum dielectric constant. In deriving Eq. (2), the paraxial and slowly varying envelope approximations have been applied.

For weak probe field, the Maxwell-Bloch (MB) Eqs. (1) and (2) can be solved perturbatively. Furthermore, we are interested in the stationary state of the system, i.e., the time derivatives in the MB equations can be neglected, which is valid for the signal field with a longer time duration. The solution of Eq. (1) up to third-order has been presented in Appendix B. With the solution, we can obtain the expression of the total optical susceptibility of the signal field, which reads {\chi }_s={\chi }_s^{(1)}+{\chi }_s^{(3)}|{\mathcal{E}}_s{|}^2. Detailed expressions of linear and nonlinear susceptibilities can be found in Appendix B. For simplicity, we assume that the signal, control, and pump fields have a wide distribution in y direction, so that their y-dependence plays no significant role during the propagation of the signal field. Thus, the term {\mathrm{\partial }}^2{\mathrm{\Omega }}_s/\mathrm{\partial }{y}^2 in Eq. (2) can be disregarded.

To get space-dependent linear and nonlinear optical potentials, we assume that the control field is modulated along x direction, i.e., {\mathrm{\Omega }}_c={\mathrm{\Omega }}_c(x). Then one has {\chi }_s^{(1)}={\chi }_s^{(1)}(x) and {\chi }_s^{(3)}={\chi }_s^{(3)}(x). Thus Eq. (2) reduces into

\mathrm{i}{\displaystyle \frac{\mathrm{\partial }{\mathrm{\Omega }}_s}{\mathrm{\partial }z}}+{\displaystyle \frac{1}{2k_s}}{\displaystyle \frac{{\mathrm{\partial }}^2{\mathrm{\Omega }}_s}{\mathrm{\partial }{x}^2}}+{\displaystyle \frac{k_s}{2}}{\chi }_s^{(1)}(x){\mathrm{\Omega }}_s+{\displaystyle \frac{b{k}_s}{2}}{\chi }_s^{(3)}(x)|{\mathrm{\Omega }}_s{|}^2{\mathrm{\Omega }}_s=0,

here b=(2\hslash /|{\mathit{p}}_{13}|{)}^2. For the convenience of the following discussion, we convert it into the dimensionless form

(3)\mathrm{i}{\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }\zeta }}=-{\displaystyle \frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{\xi }^2}}+V(\xi )u+W(\xi )|u{|}^{2}u,

where \xi =x/l_s, \zeta =z/(2k_s{l}_s^2), and u={\mathrm{\Omega }}_s/U_0, with U_0 and l_s respectively the typical half Rabi frequency and typical length in x direction (comparable to the wavelength {\lambda }_s=2\pi c/{\omega }_s) of the signal field. In Eq. (3), V(\xi )=-k_s^2{l}_s^2{\chi }_s^{(1)} and W(\xi )=-b{k}_s^2{l}_s^2|U_0{|}^2{\chi }_s^{(3)}, which are called respectively linear and nonlinear potentials, or combined linear and nonlinear potentials [1, 45].

3 Physical realization of combined \mathcal{P}\mathcal{T}-symmetric potentials

We now propose a realistic scheme to realize \mathcal{P}\mathcal{T}-symmetric potentials with combined linear and nonlinear potentials. By appropriately selecting system parameters and implementing a spatially modulated control field, we will demonstrate that such potentials can indeed be achieved in our system. Before presenting the concrete form of the spatially modulated {\mathrm{\Omega }}_c for realizing the combined linear and nonlinear optical potentials obeying \mathcal{P}\mathcal{T}-symmetry, we give the relation between the linear and nonlinear optical susceptibilities of the signal field and the control field frequency detuning {\mathrm{\Delta }}_2 for a fixed {\mathrm{\Omega }}_c={\mathrm{\Omega }}_{c0}. Other system parameters {\mathrm{\Delta }}_3={\mathrm{\Delta }}_4=0, {\mathrm{\Gamma }}_3=2\pi \times 5.75 MHz, {\mathrm{\Gamma }}_4=2\pi \times 6.07 MHz, {\mathrm{\Omega }}_s=2\pi \times 0.37 MHz, {\mathrm{\Omega }}_p=2\pi \times 4 MHz, {\mathrm{\Omega }}_{c0}=2\pi \times 2 MHz, and {\mathcal{N}}_a=1.8\times {10}^{13}{\mathrm{c}\mathrm{m}}^{-3}. Because |1⟩ and |2⟩ are two hyperfine levels of the same ground-state, the electric-dipole transition between them is forbidden, and the spontaneous population decay rate of |2⟩ can be taken as {\mathrm{\Gamma }}_2\simeq 0.

According to solving Eq. (1), we obtain the expressions of {\chi }_s^{(1)} and {\chi }_s^{(3)}, which are calculated numerically and then plotted in Figs. 2(a) and (d), respectively. The blue solid and red dashed lines in Fig. 2(a) represent the real part Re[{\chi }_s^{(1)}] and imaginary part Im[{\chi }_s^{(1)}] of the linear susceptibility, respectively. Figure 2(d) shows the nonlinear susceptibility. We observe that the imaginary parts of both the linear and nonlinear susceptibilities Im[{\chi }_s^{(1)}] and Im[{\chi }_s^{(3)}] simultaneously vanish at {\mathrm{\Delta }}_2\approx -30.07 MHz and {\mathrm{\Delta }}_2\approx 30.07 MHz, as marked by the points “P” and “Q” in Figs. 2(a) and (d), respectively. In the vicinity of “P” point, {\chi }_s^{(1)} exhibits absorption on the left side and gain on the right side, while {\chi }_s^{(3)} exhibits gain on the left side and absorption on the right side. Conversely, in the vicinity of “Q” point, {\chi }_s^{(1)} exhibits gain on the left side and absorption on the right side, whereas {\chi }_s^{(3)} exhibits absorption on the left side and gain on the right side.

Obviously, a symmetric real parts and an anti-symmetric imaginary parts of {\chi }_s^{(1)} and {\chi }_s^{(3)} can be acquired by choosing appropriate {\mathrm{\Delta }}_2 in the vicinity of “P” and “Q” points. Shown in Figs. 2(b) and (e) are Re[{\chi }_s^{(1)}] and Re[{\chi }_s^{(3)}] as functions of {\mathrm{\Omega }}_{c0}/(2\pi ). The red solid, black dashed, and blue solid lines correspond to {\mathrm{\Delta }}_2=\pm 29.73 MHz, \pm 30.07 MHz, and \pm 30.45 MHz, respectively. We see that \mathrm{R}\mathrm{e}[{\chi }_s^{(j)}]{|}_{{\mathrm{\Delta }}_2=\pm 29.73\mathrm{M}\mathrm{H}\mathrm{z}}\approx\mathrm{R}\mathrm{e}[{\chi }_s^{(j)}]{|}_{{\mathrm{\Delta }}_2=\pm 30.07\mathrm{M}\mathrm{H}\mathrm{z}}\approx\mathrm{R}\mathrm{e}[{\chi }_s^{(j)}]{|}_{{\mathrm{\Delta }}_2=\pm 30.45\mathrm{M}\mathrm{H}\mathrm{z}} (j=1,3), which means both the real parts of the linear and nonlinear susceptibilities are matched well (i.e., they are symmetric) [41, 45]. Figures 2(c) and (f) show results of Im[{\chi }_s^{(1)}] and Im[{\chi }_s^{(3)}] as functions of {\mathrm{\Omega }}_{c0}/(2\pi ). One finds that \mathrm{I}\mathrm{m}[{\chi }_s^{(j)}]{|}_{{\mathrm{\Delta }}_2=\pm 30.07\mathrm{M}\mathrm{H}\mathrm{z}}=0 and \mathrm{I}\mathrm{m}[{\chi }_s^{(j)}]{|}_{{\mathrm{\Delta }}_2=\pm 29.73\mathrm{M}\mathrm{H}\mathrm{z}}\approx -\mathrm{I}\mathrm{m}[{\chi }_s^{(j)}]{|}_{{\mathrm{\Delta }}_2=\pm 30.45\mathrm{M}\mathrm{H}\mathrm{z}} (j=1,3), which means the imaginary parts of the linear and nonlinear susceptibilities are asymmetric [45].

The above analysis demonstrates that the \mathcal{P}\mathcal{T} symmetry of the combined linear and nonlinear potentials, V(\xi ) and W(\xi ), can be achieved by employing a control field composed of three spatially separated laser beams with distinct frequency detunings. These beams generate three independent waveguides, each contributing differently to the signal field: one providing gain, another absorption, and the third remaining neutral. To fulfill the condition of \mathcal{P}\mathcal{T} symmetry, we need V^{\ast }(-\xi )=V(\xi ) and W^{\ast }(-\xi )=W(\xi ) [equivalent to Re[{\chi }_s^{(j)}(x)] =Re[{\chi }_s^{(j)}(-x)] and Im[{\chi }_s^{(j)}(x)]=-Im[{\chi }_s^{(j)}(-x)] (j=1,3)].

Based on the above analysis, we assume that the control field is composed of three Gaussian beams with identical spatial profiles. The transverse intensity distribution of the control field can therefore be expressed as

(4){\mathrm{\Omega }}_c(x)={\mathrm{\Omega }}_{c0}\left[{\mathrm{e}}^{-\frac{(x-x_0{)}^2}{2{\sigma }^2}}+{\mathrm{e}}^{-\frac{x^2}{2{\sigma }^2}}+{\mathrm{e}}^{\frac{-(x+x_0{)}^2}{2{\sigma }^2}}\right],

where three beams are located at distinct positions (\pm x_0,0) and (0,0), with \sigma denoting the width of each beam. Three beams can be generated by splitting a single coupling laser into three optical paths. Each beam passes through an independently driven acousto-optic modulator (AOM), allowing its optical frequency, and hence its coupling-field detuning, to be adjusted independently. Identical beam-shaping optics are employed to ensure that the three beams have the same peak intensity and Gaussian beam width. Here, “identical” refers to their spatial profiles and does not require their optical frequencies to be exactly the same; the small frequency shifts introduced by the AOMs enable independent control of the three detunings. Although the beams originate from the same laser, mutual interference can be neglected when their spatial overlap is sufficiently weak.

Substituting Eq. (4) into the expressions of {\chi }_s^{(j)}(x), we obtain the following \mathcal{P}\mathcal{T}-symmetric (dimensionless) linear and nonlinear potentials

(5a)\begin{array}{c}\begin{array}{rl}V(\xi )=& V_0+v_0[({\mathrm{e}}^{\frac{-(\xi -d{)}^2}{2a^2}}+{\mathrm{e}}^{\frac{-{\xi }^2}{2a^2}}+{\mathrm{e}}^{\frac{-(\xi +d{)}^2}{2a^2}})\\ & +\mathrm{i}{v}_1({\mathrm{e}}^{\frac{-(\xi -d{)}^2}{2a^2}}-{\mathrm{e}}^{\frac{-(\xi +d{)}^2}{2a^2}})],\end{array}\end{array}

(5b)\begin{array}{c}\begin{array}{rl}W(\xi )=& W_0+w_0[({\mathrm{e}}^{\frac{-(\xi -d{)}^2}{2a^2}}+{\mathrm{e}}^{\frac{-{\xi }^2}{2a^2}}+{\mathrm{e}}^{\frac{-(\xi +d{)}^2}{2a^2}})\\ & +\mathrm{i}{w}_1({\mathrm{e}}^{\frac{-(\xi -d{)}^2}{2a^2}}-{\mathrm{e}}^{\frac{-(\xi +d{)}^2}{2a^2}})],\end{array}\end{array}

where d=x_0/l_s, a=\sigma /l_s, V_0 and W_0 are constants, v_0 (w_0) is the amplitude of the space-dependent part for the linear potential V (nonlinear potential W), v_1 (w_1) is the relative amplitude between the imaginary part to the real part of V (W). The effective coefficients v_0, v_1, w_0, and w_1 introduced in Eq. (5) are not arbitrary or mutually independent, see Appendix. B. Rather, they are jointly determined and constrained by the the atomic and laser parameters, including the two-photon detuning {\mathrm{\Delta }}_2, the control and probe field Rabi frequencies {\mathrm{\Omega }}_c and {\mathrm{\Omega }}_p, the atomic density {\mathcal{N}}_a, and others.

Spatial profiles of V(\xi ) and W(\xi ) are shown in Figs. 3(a) and (b), respectively. When plotting the figures, system parameters are chosen as {\mathrm{\Delta }}_2=-30.07 MHz [“P” point in Fig. 2(a)], -29.73 MHz [right side of “P” point in Fig. 2(a)], -30.45 MHz [left side of “P” point in Fig. 2(a)], x_0=15 μm, \sigma =2.2 μm, and U_0=2\pi \times 0.68 MHz. These correspond to the dimensionless parameters v_0=2, w_0=2, v_1=0.25, w_1=0.06, d=15, and a=2.2. Other parameters are same as Fig. 2. As illustrated in Fig. 3, both the linear potential V(\xi ) and the nonlinear potential W(\xi ) exhibit the expected \mathcal{P}\mathcal{T} symmetry. For the specific form of the control field in Eq. (4), the real parts of the potentials Re(V) and Re(W) have shape of triple barriers (blue solid lines) as shown in Figs. 3(a) and (b); for the linear potential, there is an absorption on the left side [Im(], a gain on the right side [Im(V)>0], and a neutral [Im(V)=0] region in the middle. In contrast, the nonlinear potential exhibits the opposite behavior: gain on the left side [Im(W)>0], absorption on the right side [Im(], and a neutral (zero) region in the middle. These features are depicted by the red dashed lines in the figure.

Furthermore, when choosing {\mathrm{\Delta }}_2=30.07 MHz [“Q” point in Fig. 2(a)], 29.73 MHz [left side of “Q” point], and 30.45 MHz [right side of “Q” point], v_0=w_0=-2, and other parameters identical to those in Figs. 3(a) and (b). From Figs. 3(c) and (d) we see that both the linear potential V(\xi ) and the nonlinear potential W(\xi ) retain their \mathcal{P}\mathcal{T} symmetry. However, in this case, Re(V) and Re(W) adopt triple potential wells (blue solid lines). For the linear potential, there is a gain on the left side [Im(V)>0] and an absorption on the right side [Im(], and a neutral [Im(V)=0] region in the middle. However, the nonlinear potential reverses this behavior. There is an absorption on the left side [Im(] and a gain on the right side [Im(W)>0], illustrated by the red dashed lines in the figure. Therefore, it can be easily realized \mathcal{P}\mathcal{T} potential wells and barriers by selecting appropriate system parameters in atomic gases.

4 Higher-order exceptional points and controllable \mathcal{P}\mathcal{T} phase transition

4.1 Higher-order exceptional points

Now, we turn to consider the EPs in the linear \mathcal{P}\mathcal{T} symmetry system. Among these, EP2 have been the most extensively studied due to their abundance, requiring only the tuning of one real parameter. In contrast, it is argued that a close approach of an EP3 cannot be achieved with adjusting only one real parameter [7]. The ultracold atomic system is popular platform in the literature [40–43], since they allow for experimental, but are flexible enough to provide insight into characteristic phenomena of the more complex physical situations. Our system is expected to search the higher-order EPs. In fact, the finding point is to high sensitivity in the parameters near to the EP3 [7].

For v_1\ne 0, the model is underlying \mathcal{P}\mathcal{T} symmetry we except real eigenvalues until reach a coalescence of at last two eigenvalues where two eigenvalues become complex at an EP2. Such value depends on the other parameters a and d. Here we seek these parameters in such a way that the fourfold of the for parameters (v_0,v_1,a,d) leads to the coalescence to an EP3 [7]. To achieve the coalescence of three real eigenvalues at least two parameters have to be judiciously chosen. It is at this point where we would expect the need for a careful fine-tuning the physical parameters. To find the EP3, the second parameter v_0 (definite value) must be judiciously varied, with the fixed parameters a and d.

Figure 4 shows the linear eigenvalues as functions of the amplitude v_1 of the imaginary part of the linear potential, with v_0 held fixed. This effective-parameter scan determines the critical strength of the non-Hermitian component at which the three or two eigenvalues coalesce, thereby locating the EP3 or EP2. Here, “fixed v_0” denotes a constrained trajectory in the multidimensional physical-parameter space, rather than implying that v_0 and v_1 are fundamentally independent or arbitrarily tunable.

From Fig. 4(a), we observe that increasing v_1 induces the coalescence of the three real eigenvalues to form EP3. The spectrum remains entirely real for , indicating that the system preserves unbroken \mathcal{P}\mathcal{T} symmetry in this regime. However, at the critical value v_1=v_{1\mathrm{E}\mathrm{P}3}=0.022, the \mathcal{P}\mathcal{T} symmetry is broken, marked by the emergence of a pair of complex conjugate eigenvalues in Fig. 4(b).

When varying v_0, deviates from its optimal value, EP3 becomes inaccessible in the physical system. To further illustrate this, we numerically analyze the real and imaginary parts of the eigenvalues as functions of v_1, as shown in Figs. 4(c) and (d). Here, only two real eigenvalues coalesce, corresponding to the typical EP2, i.e., v_1=0.06. This confirms that achieving an EP3, where three eigenvalues coalesce, requires precise tuning of at least two independent parameters, highlighting the stringent conditions necessary for realizing higher-order EPs in realistic systems.

4.2 Bifurcations of complex eigenvalues around the EP3

We now turn to focus on sensitivity change of specture near higher EPs, which is essential for the implement of sensor devices [17, 18]. To analyze how a small perturbation affects the system, we impose the small perturbation ϵ on the linear operator \hat{L}=-{\mathrm{\partial }}_{\xi \xi }+V(\xi ), varying the perturbation ϵ from 0 to 0.4.

At EP3, we numerically compute the eigenvalues and plot their real and imaginary parts as functions of ϵ in Fig. 5. In this case, the three originally degenerate eigenvalues undergo splitting. From Fig. 5(a), we observe that the perturbation lifts the degeneracy in the real parts of the eigenvalues. Meanwhile, the imaginary parts exhibit splitting [Fig. 5(b)], confirming that the originally coalesced eigenvalues (both real and imaginary components) no longer remain degenerate. Furthermore, the difference between two eigenvalues {\mu }_1 and {\mu }_2 is also plotted as a function of ϵ in Fig. 5(c). As a result, the eigenvalues split with a characteristic scaling of {ϵ}^{1/3}. By examining the logarithmic scaling of this curve, we determine that the slope of the response is 1/3, as shown in Fig. 5(d). This confirms that perturbations near an EP3 follow an enhanced sensitivity scaling of the form {ϵ}^{1/3}.

However at the EP2, the eigenvalues exhibit a square-root dependence on perturbations, scaling as {ϵ}^{1/2} [17]. Hence near EP3, the splitting would be more abrupt, which displays a sharper cubic-root dependence, scaling as {ϵ}^{1/3}. This enhanced sensitivity at EP3 suggests significant advantages for sensing applications requiring exceptional resolution.

4.3 \mathcal{P}\mathcal{T} phase transitions and their active control

We now turn to consider the properties of the \mathcal{P}\mathcal{T} phase transition with the combined linear and nonlinear \mathcal{P}\mathcal{T} potentials. A key advantage of this configuration is its active controllability, the \mathcal{P}\mathcal{T} phase transition may be precisely manipulated by tuning system parameters [i.e., barrier and potential well parameters v_1 and w_1 in Eq. (5)]. The effective dimensionless coefficients can be varied within experimentally accessible ranges through coordinated adjustments of the atomic and laser parameters. After the critical values have been identified by the spectral analysis, the coefficients used in the subsequent calculations are selected below, near, or above these critical values to investigate the propagation and scattering dynamics in different parameter regimes.

To demonstrate behavior, we consider solutions of the form u(\zeta ,\xi )=u_0(\xi ){\mathrm{e}}^{\mathrm{i}\mu \zeta }, substituted into Eq. (3), yield the nonlinear eigenvalue problem {\mathrm{\partial }}^2{u}_0/\mathrm{\partial }{\xi }^2-V(\xi )u_0-W(\xi )|u_0{|}^2{u}_0=\mu u_0, where \mu represents the eigenvalue (also called propagation constant) and u_0(\xi ) is corresponding eigenfunction. We solve this problem numerically using the Newton iteration method [60]. One of key characters of \mathcal{P}\mathcal{T} symmetry is that the eigenvalue displays a transition from purely real to complex as the relative amplitudes of the linear and nonlinear potentials (i.e., v_1 and w_1) are varied. This transition marks the spontaneous breaking of \mathcal{P}\mathcal{T} symmetry and can be systematically controlled by adjusting these parameters.

Figure 6(a) displays the result for the phase diagram of the \mathcal{P}\mathcal{T} phase transition in the parameter plane of v_1 and w_1 for fixed a=2.2, d=15, v_0=2, and w_0=2, with other system parameters the same as used in Fig. 3(a) and (b). The purple solid line in the figure represents the boundary line of the \mathcal{P}\mathcal{T} phase transition, with the left-bottom (light blue) domain being the phase with \mathcal{P}\mathcal{T} symmetry and the top-right (orangey) domain being the phase with broken \mathcal{P}\mathcal{T} symmetry. The diagram reveals that the \mathcal{P}\mathcal{T} phase transition critically depends not only on v_1 but also on w_1. Notably, the strength of the imaginary nonlinear potential W in Eq. (3) significantly influences the symmetry breaking threshold, demonstrating the crucial interplay between linear and nonlinear effects in governing the system’s \mathcal{P}\mathcal{T} behavior. To be specific, the phase transition will be decreased when w_1 increases in Fig. 6(a). The \mathcal{P}\mathcal{T} phase transition is determined by the ratio between the imaginary and real parts of the combined linear and nonlinear potentials, which is determined by the quantity (v_0{v}_1+w_0{w}_1|u{|}^2)/(v_0+w_0|u{|}^2), with |u{|}^2 the intensity of the signal field. Therefore, the \mathcal{P}\mathcal{T} phase transition depends on the parameters v_0, v_1, w_0, w_1, and |u{|}^2.

Except for v_1 and w_1, the \mathcal{P}\mathcal{T} phase diagram is also tunable by the beam width a of the control field. Shown in Fig. 6(b) are phase boundary lines of \mathcal{P}\mathcal{T} phase transition for different values of a, with other system parameters the same as those used in Fig. 6(a), where the purple solid and purple dash-dotted lines are for a=2.2 and a=2.1, respectively. The comparison reveals a significant expansion of the \mathcal{P}\mathcal{T} symmetry phase as a decreases. This trend is explicitly indicated by the purple arrow, which shows the systematic shift of the phase boundary with varying beam width. The results demonstrate that narrower beam profiles (smaller a) enhance the stability of the \mathcal{P}\mathcal{T}-symmetric phase in the system.

The threshold of the \mathcal{P}\mathcal{T} symmetry breaking is further influenced by the beam separation (characterizing by the parameter d) of the control field. Figure 6(c) illustrates this dependence, comparing the \mathcal{P}\mathcal{T} phase diagrams for d=15 (solid purple line) and d=16 (dash-dotted purple line), while keeping all other system parameters identical to those in Fig. 6(a). A clear expansion of the \mathcal{P}\mathcal{T}-symmetric phase (unbroken phase) is observed as d increases. This demonstrates that larger beam separations enhance the stability of the \mathcal{P}\mathcal{T}-symmetric regime, effectively raising the threshold for symmetry breaking.

Another key characters can change \mathcal{P}\mathcal{T} symmetry domain is that the \mathcal{P}\mathcal{T} potential displays a barriers and potential wells. Shown in Fig. 6(d) is the phase diagram of the \mathcal{P}\mathcal{T} phase transition in the plane of v_1 and w_1 for v_0=-2, and w_0=-2 (a potenial wells), with other system parameters the same as used in Fig. 6(a). The purple solid line in the figure represents the boundary line of the \mathcal{P}\mathcal{T} phase transition, with the left-bottom domain being the phase with \mathcal{P}\mathcal{T} symmetry and the top-right domain being the phase with broken \mathcal{P}\mathcal{T} symmetry. From the figures, we see that the sign of the linear and nonlinear potentials in Eq. (3) can extend the region of unbroken \mathcal{P}\mathcal{T} phase. It plays an important role for the application of \mathcal{P}\mathcal{T} phase transition.

Furthermore, the \mathcal{P}\mathcal{T} phase diagram can also be changed by adjusting the width a and the separation d of the beams of the control field. Shown in Fig. 6(e) [(f)] are phase boundary lines of \mathcal{P}\mathcal{T} phase transition for different values of a [d] [with other system parameters are same as those used in Fig. 6(d)], where the purple solid and purple dash-dotted lines are for a=2.2 and a=2.0 [d=15 and d=16], respectively. From the figure, we see that the domain of the \mathcal{P}\mathcal{T} symmetry phase is increased greatly as a is decreased. The domain of the \mathcal{P}\mathcal{T} symmetry phase is increased when d increases.

From the results presented in Fig. 6, we conclude that the \mathcal{P}\mathcal{T}-phase transition can be effectively controlled by tuning the relative amplitudes between the imaginary and real parts of the combined linear and nonlinear potentials, as well as by modifying the structure of the \mathcal{P}\mathcal{T}-symmetric potential (e.g., barriers versus potential wells). These findings provide a concrete example of active control in non-Hermitian nonlinear optical systems.

4.4 Stable solitons and their active control

One of the most important applications of controlling the \mathcal{P}\mathcal{T}-phase transition is the propagation and active manipulation of solitons in systems with combined linear and nonlinear \mathcal{P}\mathcal{T}-symmetric potentials. By numerically solving Eq. (3), we demonstrate that the system supports stable optical solitons when operating in the \mathcal{P}\mathcal{T}-symmetric regime. However, these solitons become unstable when the system enters the broken \mathcal{P}\mathcal{T}-symmetry phase. These results highlight the critical role of \mathcal{P}\mathcal{T} symmetry in governing soliton dynamics and their controllability in non-Hermitian optical systems.

We first study the optical soliton with three barriers \mathcal{P}\mathcal{T} symmetry potential. To investigate the system’s dynamics, we perform numerical simulations of Eq. (3) using the ground state solution of the linear eigenequation as the initial condition. In this case, the regions between the three potential barriers effectively form two potential wells; consequently, the ground state solution is distributed over the two wells. Shown in Fig. 7(a) is the numerical result on the propagation of optical soliton, with (w_1,v_1)=(5,0) to make the system work in the \mathcal{P}\mathcal{T}-symmetric phase [i.e., the point “A” in Fig. 6(a)]. We see that the optical soliton is fairly stable during propagation, with their intensity arrested in two spatial region between the adjacent Gaussian peaks of the control field. In contrast, Fig. 7(b) shows the propagation of another optical soliton for (w_1,v_1)=(35,0.085), corresponding to point “B” in Fig. 6(a), where the system is in the broken \mathcal{P}\mathcal{T}-symmetric phase. Here, the potential on the left side of the three barriers structure is of the gain channel. The ground state solution is thus predominantly localized on the left side, while the associated gain channel gives rise to pronounced instability of the soliton.

Further simulations in Fig. 7(c) examine soliton dynamics in a triple-well potential with (w_1,v_1)=(10,0.05), maintaining \mathcal{P}\mathcal{T} symmetry. The soliton shows remarkable stability, with intensity confinement at the center of the potential well. However, when the parameters are adjusted to (w_1,v_1)=(10,0.2) (broken \mathcal{P}\mathcal{T} phase), Fig. 7(d) reveals pronounced instability in the soliton propagation.

5 Study on the soliton scattering

5.1 Symmetric soliton scattering

The scattering properties of solitons are of fundamental interest for deep understanding of \mathcal{P}\mathcal{T}-symmetric systems and also for possible practical applications. It is natural to ask the question how about the soliton scattering in our optically-pumped EIT system if the combined \mathcal{P}\mathcal{T}-symmetric optical potentials act as a defect. To address this, we initialize the soliton at a position far from the defect, ensuring no initial interaction at \zeta =0. The scattering dynamics may lead to diverse outcomes, including full reflection, transmission, trapping, or a combination.

These scattering behaviors can be described by the coefficients of reflection (R), transmission (T), and trapping (G), defined respectively by

\begin{array}{rl}& R=\frac{1}{Q}{\int }_{-\mathrm{\infty }}^{-X}|\psi (\xi ){|}^2\mathrm{d}\xi ,\\ & T=\frac{1}{Q}{\int }_X^{+\mathrm{\infty }}|\psi (\xi ){|}^2\mathrm{d}\xi ,\\ & G=\frac{1}{Q}{\int }_{-X}^X|\psi (\xi ){|}^2\mathrm{d}\xi ,\end{array}

where Q={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}|\psi (\xi ){|}^2\mathrm{d}\xi is the total power of the optical soliton and X denotes position on the \xi-axis at which the influence of the defect to the soliton is negligible.

We first consider the situation with (w_1,v_1)=(0,0), corresponding to soliton scattering by the pure real combined linear and nonlinear potentials. Other parameters a=2.2, d=15, v_0=2, and w_0=2. That is, we consider the scattering problem associated with the potential barriers. Shown in Figs. 8(a)−(c) present the results of the soliton scattering when the soliton for different incident velocities (\nu = 0.4, 0.93, and 1.5, respectively) when the soliton is incident from the left side of the defect. These results are obtained through numerically solving Eq. (3) by using the split-step Fourier method [60] and taking the initial condition u(\xi ,\zeta =0)=0.5\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{h}(2\xi )e^{i\nu \xi } [45, 58]. In each panel, the defect region is indicated by the area between three vertical white dashed lines. The results demonstrate distinct scattering behaviors across the different velocity regimes, revealing the critical dependence of soliton-defect interactions on the incident velocity. We see that, for smaller (larger) incident velocity, the soliton is completely reflected (transmitted), and a small interval of incident velocity, i.e., , the soliton is partially trapped.

Figure 8(d) illustrates the result of the reflection coefficient R (blue solid line), trapping coefficient G (red dash-dotted line) and transmission coefficient T (green dashed line) as functions of incident velocity \nu. The purple dots “a”, “b”, and “c” in the figure indicate the values of R, G, and T, which correspond to the cases shown in panels (a), (b), and (c), respectively. We find that when \nu \ge {\nu }_{cr} (where {\nu }_{cr}\approx 0.92 is a critical value), the scattering of the soliton changes sharply from a full reflection to trapping, and then the scattering of the soliton changes sharply from trapping to full transmission (where \nu >0.94 is a critical value).

For comparison, in Figs. 8(e)−(g), we show the result of the soliton scattering when the soliton is incident from right side of the defect with incident velocity \nu =0.4, 0.93, and 1.5, respectively. One sees that for smaller (larger) incident velocity the soliton is also completely reflected (transmitted). Figure 8(h) plots the reflection coefficient R (blue solid line), trapping coefficient G (red dash-dotted line), and transmission coefficient T (green dashed line) as functions of incident velocity \nu. The purple dots “e”, “f”, and “g” indicate the values of R, G, and T for the cases shown in panels (e), (f), and (g), respectively. Similar to the case of the scattering from left side, a small interval of incident velocity, i.e., , the soliton is partially trapped. These results collectively demonstrate that the combined real linear and nonlinear \mathcal{P}\mathcal{T}-symmetric defect potentials exhibit left-right symmetric scattering behavior. The identical response to solitons incident from either direction confirms the fundamental symmetry of the system.

5.2 Asymmetric soliton scattering

We now investigate what will happen if the imaginary parts of the combined linear and nonlinear \mathcal{P}\mathcal{T} defect potentials are not zero by taking (w_1,v_1)=(0,0.01). Panels (a), (b), and (c) of Fig. 9 show results of soliton scattering when the soliton is incident from the left side of the defect with \nu =0.5,0.93, and 1.5, respectively. We see that the soliton gets a complete reflection when it collides with the defect for a small incidence velocity [\nu =0.5; panel (a)], and a complete transmission for a large incidence velocity [\nu =1.5; panel (c)]. However, for an intermediate incidence velocity, the soliton experiences a state with a combination of reflection and transmission [\nu =0.93; panel (b)], in which a radiation is generated.

Plotted in Fig. 9(d) are reflection coefficient R (blue solid line), trapping coefficient G (red dash-dotted line), and transmission coefficient T (green dashed line) as functions of \nu. Purple dots indicate the value of R, G, and T corresponding to the panels (a), (b) and (c), respectively. We see that there exists an interval of the incident velocity, i.e., , in which the trapping coefficient G is nonzero.

Shown in panels (e), (f), and (g) of Fig. 9 are results of the soliton scattering when the soliton is incident from the right side of the defect with the incident velocity \nu =0.5, 1.05, and 1.5, respectively. One sees that for small and large incident velocity [panels (e) and (g)] the soliton scattering displays similar behaviors like those shown in panels (a) and (c), but there exists an large interval, i.e., , in which the trapping coefficient G is nonzero. Thus by tuning \nu one can control the ratio between the reflected and transmitted parts of the soliton, which might be useful to design an optical soliton beam splitter. Figure 9(h) gives the reflection coefficient R (blue solid line), trapping coefficient G (red dash-dotted line) and transmission coefficient T (green dashed line) as functions of v. Purple dots indicate the values of R and T corresponding to the panels (e), (f) and (g), respectively.

A comparison of the upper and lower panels of Fig. 9 reveals that the left–right symmetry of soliton scattering is broken. Specifically, the scattering behavior depends on whether the input pulse is launched from the left or the right side of the defect. Since the imaginary parts of both the linear and nonlinear \mathcal{P}\mathcal{T}-symmetric potentials are odd functions, a soliton incident from the left sequentially encounters the loss, neutral, and gain channels. In contrast, a soliton incident from the right experiences these regions in the reverse order, namely, gain, neutral, and loss. This reversal of the gain-loss sequence gives rise to distinct scattering outcomes for left- and right-incident solitons.

Lastly, we investigate the soliton scattering by a complex \mathcal{P}\mathcal{T}-symmetric linear and nonlinear defect potentials with (w_1,v_1)=(0.1,0.01). Figures 10(a)−(c) present the scattering dynamics for left-incident solitons with velocities \nu =0.5, 1.0, and 1.5, respectively. The corresponding scattering coefficients are quantified in Fig. 10(d), which plots the reflection (R, blue solid line), trapping (G, red dash-dotted line), and transmission (T, green dashed line) coefficients as functions of incident velocity. The purple markers identify the specific cases shown in panels (a)−(c). The results for right-incident solitons are shown in Figs. 10(e)−(g) for velocities \nu =0.5, 1.5, and 2.0, respectively. The velocity-dependent scattering coefficients for this configuration are presented in Fig. 10(h), maintaining the same color and line-style conventions, with purple markers again indicating the cases shown in panels (e)−(g).

The numerical results reveal distinct velocity-dependent scattering behaviors in Fig. 10. For both left and right incidence, we observe similar scattering characteristics at extremely low and high velocities. However, significant asymmetry emerges in the intermediate velocity regime (), where the scattering dynamics differ markedly depending on the incidence direction. Figures 10(b) and (f) show the soliton scattering for \nu =1.0 and \nu =1.5, where the soliton is completely collapse. In fact, no complete self-trapping is found for the soliton when it collides the defect from right side; but there exists an interval of \nu, i.e., . These findings demonstrate that the introduction of non-zero imaginary components in both linear and nonlinear potentials breaks the left-right symmetry of soliton scattering. The observed asymmetry fundamentally distinguishes \mathcal{P}\mathcal{T}-symmetric systems from their purely real counterparts, highlighting the unique role of balanced gain and loss in soliton-defect interactions.

6 Summary

In this work, we have proposed a realistic physical scheme for realizing triple-well potential with \mathcal{P}\mathcal{T} symmetry and investigated the higher-order EPs. Our system generates combined linear and nonlinear \mathcal{P}\mathcal{T}-symmetric potentials through controlled spatial modulation of laser fields. The nonlinear potential’s imaginary component proves crucial for controlling \mathcal{P}\mathcal{T} phase transitions and modifying phase diagrams. We have shown that the combined linear and nonlinear \mathcal{P}\mathcal{T}-symmetric potentials can be produced through the design of the spatial modulation of the control laser field. We have demonstrated that the imaginary part of the nonlinear \mathcal{P}\mathcal{T} potential plays a very important role for the occurrence of \mathcal{P}\mathcal{T} phase transition and the change of \mathcal{P}\mathcal{T} phase diagram, which can be actively manipulated in our system. The system supports stable, tunable optical solitons through \mathcal{P}\mathcal{T} potential adjustments. When these potentials act as defects, they induce asymmetric soliton scattering behaviors that can be controlled via the imaginary parts of both linear and nonlinear potentials. These findings offer potential applications in optical switch, beam splitter, optical sensing and information transmission technologies.

7 Appendix A: Optical Bloch equation

The optical Bloch equation describing the time evolution of the density-matrix elements {\rho }_{jl} reads [61]

\begin{array}{rl}& \mathrm{i}\frac{\mathrm{\partial }}{\mathrm{\partial }t}{\rho }_{22}-i{\mathrm{\Gamma }}_{24}{\rho }_{44}-i{\mathrm{\Gamma }}_{23}{\rho }_{33}+{\mathrm{\Omega }}_c^{\ast }{\rho }_{32}-{\mathrm{\Omega }}_c{\rho }_{32}^{\ast }=0,\\ & \mathrm{i}(\frac{\mathrm{\partial }}{\mathrm{\partial }t}+{\mathrm{\Gamma }}_3){\rho }_{33}+{\mathrm{\Omega }}_s{\rho }_{31}^{\ast }-{\mathrm{\Omega }}_s^{\ast }{\rho }_{31}+{\mathrm{\Omega }}_c{\rho }_{32}^{\ast }-{\mathrm{\Omega }}_c^{\ast }{\rho }_{32}=0,\\ & \mathrm{i}\frac{\mathrm{\partial }}{\mathrm{\partial }t}{\rho }_{44}+i{\mathrm{\Gamma }}_4{\rho }_{44}+{\mathrm{\Omega }}_p{\rho }_{41}^{\ast }-{\mathrm{\Omega }}_p^{\ast }{\rho }_{41}=0,\\ & \left(\mathrm{i}\frac{\mathrm{\partial }}{\mathrm{\partial }t}+d_{21}\right){\rho }_{21}+{\mathrm{\Omega }}_c^{\ast }{\rho }_{31}-{\mathrm{\Omega }}_s{\rho }_{32}^{\ast }-{\mathrm{\Omega }}_p{\rho }_{42}^{\ast }=0,\\ & \left(\mathrm{i}\frac{\mathrm{\partial }}{\mathrm{\partial }t}+d_{31}\right){\rho }_{31}+{\mathrm{\Omega }}_s({\rho }_{11}-{\rho }_{33})+{\mathrm{\Omega }}_c{\rho }_{21}-{\mathrm{\Omega }}_p{\rho }_{43}^{\ast }=0,\\ & \left(\mathrm{i}\frac{\mathrm{\partial }}{\mathrm{\partial }t}+d_{41}\right){\rho }_{41}+{\mathrm{\Omega }}_p({\rho }_{11}-{\rho }_{44})-{\mathrm{\Omega }}_s{\rho }_{43}=0,\\ & \left(\mathrm{i}\frac{\mathrm{\partial }}{\mathrm{\partial }t}+d_{32}\right){\rho }_{32}+{\mathrm{\Omega }}_c({\rho }_{22}-{\rho }_{33})+{\mathrm{\Omega }}_s{\rho }_{21}^{\ast }=0,\\ & \left(\mathrm{i}\frac{\mathrm{\partial }}{\mathrm{\partial }t}+d_{42}\right){\rho }_{42}-{\mathrm{\Omega }}_c{\rho }_{43}+{\mathrm{\Omega }}_p{\rho }_{21}^{\ast }=0,\\ & \left(\mathrm{i}\frac{\mathrm{\partial }}{\mathrm{\partial }t}+d_{43}\right){\rho }_{43}+{\mathrm{\Omega }}_p{\rho }_{31}^{\ast }-{\mathrm{\Omega }}_s^{\ast }{\rho }_{41}-{\mathrm{\Omega }}_c^{\ast }{\rho }_{42}=0,\end{array}

where the symbol “*” denotes complex conjugate, d_{jl}={\mathrm{\Delta }}_j-{\mathrm{\Delta }}_l+i{\gamma }_{jl}, {\gamma }_{jl}=({\mathrm{\Gamma }}_j+{\mathrm{\Gamma }}_l)/2+{\gamma }_{jl}^{\mathrm{c}\mathrm{o}\mathrm{l}} (j\ne l), and , with {\mathrm{\Gamma }}_{jl} the spontaneous emission decay rate and {\gamma }_{jl}^{\mathrm{c}\mathrm{o}\mathrm{l}} the dephasing rate from |l⟩ to |j⟩. The atomic population in the ground state {\rho }_{11} can be obtained by using the condition \sum _{j=1}^4{\rho }_{jj}=1 [61].

8 Appendix B: Solutions of the Bloch equation

By taking {\mathrm{\Omega }}_s\sim ϵ as an expansion parameter, assuming the expansion {\rho }_{jl}={\rho }_{jl}^{(0)}+ϵ{\rho }_{jl}^{(1)}+{ϵ}^2{\rho }_{jl}^{(2)}+\cdots, and substituting the expansion into the Bloch Eq. (1), we obtain a chain of linear but inhomogeneous equations, which can be solved order by order.

At the zero-order, the solution of {\rho }_{jj}^{(0)} (j=1,2,3,4), {\rho }_{32}^{(0)}, and {\rho }_{41}^{(0)} are given by

(B1a)\begin{array}{c}{\rho }_{11}^{(0)}={\displaystyle \frac{{\mathrm{\Gamma }}_{13}{X}_{32}|{\mathrm{\Omega }}_c{|}^2(\mathrm{i}{\mathrm{\Gamma }}_4+X_{41}|{\mathrm{\Omega }}_p{|}^2)}{D}},\end{array}

(B1b)\begin{array}{c}{\rho }_{22}^{(0)}={\displaystyle \frac{{\mathrm{\Gamma }}_{24}{X}_{41}|{\mathrm{\Omega }}_p{|}^2(\mathrm{i}{\mathrm{\Gamma }}_3+X_{32}|{\mathrm{\Omega }}_c{|}^2)}{D}},\end{array}

(B1c)\begin{array}{c}{\rho }_{33}^{(0)}={\displaystyle \frac{{\mathrm{\Gamma }}_{24}{X}_{41}{X}_{32}|{\mathrm{\Omega }}_p{|}^2|{\mathrm{\Omega }}_c{|}^2}{D}},\end{array}

(B1d)\begin{array}{c}{\rho }_{44}^{(0)}={\displaystyle \frac{{\mathrm{\Gamma }}_{13}{X}_{41}{X}_{32}|{\mathrm{\Omega }}_p{|}^2|{\mathrm{\Omega }}_c{|}^2}{D}},\end{array}

(B1e)\begin{array}{c}{\rho }_{32}^{(0)}=-{\displaystyle \frac{\mathrm{i}{\mathrm{\Gamma }}_3{\mathrm{\Gamma }}_{24}{X}_{41}|{\mathrm{\Omega }}_p{|}^2{\mathrm{\Omega }}_c}{d_{32}D}},\end{array}

(B1f)\begin{array}{c}{\rho }_{41}^{(0)}=-{\displaystyle \frac{\mathrm{i}{\mathrm{\Gamma }}_{13}{\mathrm{\Gamma }}_4{X}_{32}|{\mathrm{\Omega }}_c{|}^2{\mathrm{\Omega }}_p}{d_{41}D}},\end{array}

where D=2({\mathrm{\Gamma }}_{13}+{\mathrm{\Gamma }}_{24})X_{32}{X}_{41}|{\mathrm{\Omega }}_p{|}^2|{\mathrm{\Omega }}_c{|}^2+\mathrm{i}{\mathrm{\Gamma }}_4{\mathrm{\Gamma }}_{31}{X}_{32}|{\mathrm{\Omega }}_c{|}^2+\mathrm{i}{\mathrm{\Gamma }}_3{\mathrm{\Gamma }}_{24}{X}_{41}|{\mathrm{\Omega }}_p{|}^2, with X_{41}=1/d_{41}^{\ast }-1/d_{41} and X_{32}=1/d_{32}^{\ast }-1/d_{32}. Other {\rho }_{jl}^{(0)} are zero.

The first-order solution {\rho }_{31}^{(1)} is given by

(B2)\begin{array}{c}\begin{array}{rl}{\rho }_{31}^{(1)}=& \frac{a_1({\rho }_{33}^{(0)}-{\rho }_{11}^{(0)})+a_2{\mathrm{\Omega }}_p{\rho }_{41}^{\ast (0)}-a_3{\mathrm{\Omega }}_c{\rho }_{32}^{\ast (0)}}{a_1{d}_{31}+a_2|{\mathrm{\Omega }}_p{|}^2-a_3|{\mathrm{\Omega }}_c{|}^2}{\mathrm{\Omega }}_s\\ \equiv & a_{31}^{(1)}{\mathrm{\Omega }}_s,\end{array}\end{array}

with a_1=d_{21}{d}_{42}^{\ast }{d}_{43}^{\ast }-d_{21}|{\mathrm{\Omega }}_c{|}^2+d_{43}^{\ast }|{\mathrm{\Omega }}_p{|}^2, a_2=|{\mathrm{\Omega }}_p{|}^2-|{\mathrm{\Omega }}_c{|}^2+d_{21}{d}_{42}^{\ast }, and a_3=d_{42}^{\ast }{d}_{43}^{\ast }+|{\mathrm{\Omega }}_p{|}^2-|{\mathrm{\Omega }}_c{|}^2. Explicit expressions of {\rho }_{21}^{(1)}, {\rho }_{42}^{(1)}, and {\rho }_{43}^{(1)} are very lengthy and are not explicitly written here [45]. Other {\rho }_{jl}^{(1)} are zero.

Explicit expressions of {\rho }_{jl}^{(2)} are not explicitly written here [45]. The third-order solution for {\rho }_{31}^{(3)} is given by

(B3)\begin{array}{c}\begin{array}{rl}{\rho }_{31}^{(3)}=& \frac{a_1({\rho }_{33}^{(2)}-{\rho }_{11}^{(2)})+a_2{\mathrm{\Omega }}_p{\rho }_{41}^{\ast (2)}-a_3{\mathrm{\Omega }}_c{\rho }_{32}^{\ast (2)}}{a_1{d}_{31}+a_2|{\mathrm{\Omega }}_p{|}^2-a_3|{\mathrm{\Omega }}_c{|}^2}|{\mathrm{\Omega }}_s{|}^2{\mathrm{\Omega }}_s\\ \equiv & a_{31}^{(3)}|{\mathrm{\Omega }}_s{|}^2{\mathrm{\Omega }}_s.\end{array}\end{array}

With the solution, we can obtain the expression of the total optical susceptibility of the signal field, which reads {\chi }_s={\chi }_s^{(1)}+{\chi }_s^{(3)}|{\mathcal{E}}_s{|}^2. Here {\chi }_s^{(1)}={\mathcal{N}}_a|{\mathit{e}}_s\cdot {\mathit{p}}_{13}{|}^2{a}_{31}^{(1)}/({\epsilon }_0\hslash ) and {\chi }_s^{(3)}={\mathcal{N}}_a|{\mathit{e}}_s\cdot {\mathit{p}}_{13}{|}^4{a}_{31}^{(3)}/({\epsilon }_0{\hslash }^3) are respectively first-order linear and third-order nonlinear optical susceptibilities, with explicit expressions of a_{31}^{(1)} and a_{31}^{(3)} are respectively given by Eqs. (B2) and (B3).

Substituting Eq. (4) into the expressions of {\chi }_s^{(j)}(x), we obtain

(B4a)\begin{array}{c}{\chi }_s^{(1)}=\mathrm{R}\mathrm{e}({\chi }_s^{(1)})+\mathrm{i}\mathrm{I}\mathrm{m}({\chi }_s^{(1)}),\end{array}

(B4b)\begin{array}{c}{\chi }_s^{(3)}=\mathrm{R}\mathrm{e}({\chi }_s^{(3)})+\mathrm{i}\mathrm{I}\mathrm{m}({\chi }_s^{(3)}).\end{array}

The two terms on the right-hand side of Eq. (B4) represent the spatial distributions of the real and imaginary parts of the susceptibility, respectively. Then Substituting Eq. (B4) into linear and nonlinear potentials

V(\xi )=-k_s^2{l}_s^2{\chi }_s^{(1)},W(\xi )=-b{k}_s^2{l}_s^2|U_0{|}^2{\chi }_s^{(3)},

which appear in the dimensionless nonlinear Schrödinger equation (3). To relate these expressions to Eq. (5), we define the characteristic amplitudes of their spatially varying parts as

{\chi }_{s,R}^{(j)}=max\left|\text{Re}\left[{\chi }_s^{(j)}(\xi )\right]\right|,{\chi }_{s,I}^{(j)}=max\left|\text{Im}\left[{\chi }_s^{(j)}(\xi )\right]\right|.

Comparison with Eq. (5) then gives

\begin{array}{l}{v}_0=-k_s^2{l}_s^2{\chi }_{s,R}^{(1)},w_0=-b{k}_s^2{l}_s^2|U_0{|}^2{\chi }_{s,R}^{(3)},\\ v_1={\chi }_{s,I}^{(1)}/{\chi }_{s,R}^{(1)},w_1={\chi }_{s,I}^{(3)}/{\chi }_{s,R}^{(3)}.\end{array}

Thus, v_0 and w_0 determine the overall amplitudes of the spatially varying linear and nonlinear potentials, respectively, whereas v_1 and w_1 characterize the relative strengths of their imaginary and real components. For the configuration considered here, the susceptibilities contain no spatially uniform background components. Therefore the constant terms in both the linear and nonlinear potentials follow that V_0=W_0=0.

References

[1] V. V. Konotop , J. Yang , and D. A. Zezyulin , Nonlinear waves in PT-symmetric system, Rev. Mod. Phys. 88(3), 035002 (2016)

[2] C. M. Bender and D. W. Hook , PT-symmetric quantum mechanics, Rev. Mod. Phys. 96(4), 045002 (2024)

[3] H. Cartarius , J. Main , and G. Wunner , Exceptional points in atomic spectra, Phys. Rev. Lett. 99(17), 173003 (2007)

[4] S. Klaiman , U. Günther , and N. Moiseyev , Visualization of branch points in PT-symmetric waveguides, Phys. Rev. Lett. 101(8), 080402 (2008)

[5] M. A. Miri and A. Alù , Exceptional points in optics and photonics, Science 363(6422), eaar7709 (2019)

[6] W. Liu , Y. Wu , C. K. Duan , X. Rong , and J. Du , Dynamically encircling an exceptional point in a real quantum system, Phys. Rev. Lett. 126(17), 170506 (2021)

[7] W. D. Heiss and G. Wunner , A model of three coupled wave guides and third order exceptional points, J. Phys. A Math. Theor. 49(49), 495303 (2016)

[8] J. Schnabel , H. Cartarius , J. Main , G. Wunner , and W. D. Heiss , PT-symmetric waveguide system with evidence of a third-order exceptional point, Phys. Rev. A 95(5), 053868 (2017)

[9] K. Ding , G. Ma , M. Xiao , Z. Q. Zhang , C. T. Chan , Emergence , and coalescence , and topological properties of multiple exceptional points and their experimental realization, Phys. Rev. X 6(2), 021007 (2016)

[10] Q. Zhong , D. N. Christodoulides , M. Khajavikhan , K. G. Makris , and R. El-Ganainy , Power-law scaling of extreme dynamics near higher-order exceptional points, Phys. Rev. A 97, 020105(R) (2018)

[11] A. Laha , D. Beniwal , S. Dey , A. Biswas , and S. Ghosh , Third-order exceptional point and successive switching among three states in an optical microcavity, Phys. Rev. A 101(6), 063829 (2020)

[12] I. Mandal and E. J. Bergholtz , Symmetry and higherorder exceptional points, Phys. Rev. Lett. 127(18), 186601 (2021)

[13] Y. Zhang , S. Xia , X. Zhao , L. Qin , X. Feng , W. Qi , Y. Jiang , H. Lu , D. Song , L. Tang , Z. Zhu , W. Liu , and Y. Liu , Symmetry-protected third-order exceptional points in staggered flatband rhombic lattices, Photon. Res. 11(2), 225 (2023)

[14] H. Hodaei , A. U. Hassan , S. Wittek , H. Garcia-Gracia , R. El-Ganainy , D. N. Christodoulides , and M. Khajavikhan , Enhanced sensitivity at higher-order exceptional points, Nature 548(7666), 187 (2017)

[15] Y. Wu , P. Zhou , T. Li , W. Wan , and Y. Zou , High-order exceptional point based optical sensor, Opt. Express 29(4), 6080 (2021)

[16] X. Wang , G. Guo , and J. Berakdar , Enhanced sensitivity at magnetic high-order exceptional points and topological energy transfer in magnonic planar waveguides, Phys. Rev. Appl. 15(3), 034050 (2021)

[17] W. Chen , Ş. Kaya Özdemir , G. Zhao , J. Wiersig , and L. Yang , Exceptional points enhance sensing in an optical microcavity, Nature 548(7666), 192 (2017)

[18] S. Yu , Y. Meng , J. S. Tang , X. Y. Xu , Y. T. Wang , P. Yin , Z. J. Ke , W. Liu , Z. P. Li , Y. Z. Yang , G. Chen , Y. J. Han , C. F. Li , and G. C. Guo , Experimental investigation of quantum PT-enhanced sensor, Phys. Rev. Lett. 125(24), 240506 (2020)

[19] Q. Zhong , J. Kou , Ş. K. Özdemir , and R. El-Ganainy , Hierarchical construction of higher-order exceptional points, Phys. Rev. Lett. 125(20), 203602 (2020)

[20] C. Hang and G. Huang , Parity−time symmetry with coherent atomic gases: A concise review on recent progress, Adv. Phys. 74(1−4), 97 (2025)

[21] Z. Lin , H. Ramezani , T. Eichelkraut , T. Kottos , H. Cao , and D. Christodoulides , Unidirectional invisibility induced by PT-symmetric periodic structures, Phys. Rev. Lett. 106(21), 213901 (2011)

[22] L. Feng , Y.L. Xu , W. S. Fegadolli , M.H. Lu , J. E. B. Oliveira , V. R. Almeida , Y.F. Chen , and A. Scherer , Experimental demonstration of a unidirectional reflectionless parity−time metamaterial at optical frequencies, Nat. Phys. 12, 108 (2012)

[23] C. Bender , S. Factor , J. D. Bodyfelt , H. Ramezani , D. N. Christodoulides , F. M. Ellis , and T. Kottos , Observation of asymmetric transport in structures with active nonlinearities, Phys. Rev. Lett. 110(23), 234101 (2013)

[24] B. Peng , S. K. Özdemir , F. Lei , F. Monifi , M. Gianfreda , G. L. Long , S. Fan , F. Nori , C. M. Bender , and L. Yang , Parity−time-symmetric whispering-gallery microcavities, Nat. Phys. 10(5), 394 (2014)

[25] L. Chang , X. Jiang , S. Hua , C. Yang , J. Wen , L. Jiang , G. Li , G. Wang , and M. Xiao , Parity−time symmetry and variable optical isolation in active-passive-coupled microresonators, Nat. Photonics 8(7), 524 (2014)

[26] X. Zhang , T. Sun-Huo , and X. Zhang , Symmetry-guided nonreciprocal transport in non-hermitian scattering systems: Parity-flux reversal effects, Front. Phys. (Beijing) 21(4), 042203 (2026)

[27] Y. P. Wu , G. Q. Zhang , C. X. Zhang , J. Xu , and D. W. Zhang , Interplay of nonreciprocity and nonlinearity on mean-field energy and dynamics of a Bose–Einstein condensate in a double-well potential, Front. Phys. (Beijing) 17(4), 42503 (2022)

[28] S. Longhi , PT-symmetric laser absorber, Phys. Rev. A 82, 031801(R) (2010)

[29] Y. D. Chong , L. Ge , and A. D. Stone , PT-symmetry breaking and laser-absorber modes in optical scattering systems, Phys. Rev. Lett. 106(9), 093902 (2011)

[30] Y. Sun , W. Tan , H. Li , J. Li , and H. Chen , Experimental demonstration of a coherent perfect absorber with PT phase transition, Phys. Rev. Lett. 112(14), 143903 (2014)

[31] C. Wang , W. R. Sweeney , A. D. Stone , and L. Yang , Coherent perfect absorption at an exceptional point, Science 373(6560), 1261 (2021)

[32] V. V. Konotop , V. S. Shchesnovich , and D. A. Zezyulin , Giant amplification of modes in parity−time symmetric waveguides, Phys. Lett. A 376(42−43), 2750 (2012)

[33] L. Feng , Z. J. Wong , R. Ma , Y. Wang , and X. Zhang , Single mode laser by parity−time symmetry breaking, Science 346(6212), 972 (2014)

[34] H. Hodaei , M. A. Miri , M. Heinrich , D. N. Christodoulides , and M. Khajavikhan , Parity−time symmetric microring lasers, Science 346(6212), 975 (2014)

[35] L. Xiao , X. Zhan , Z. H. Bian , K. K. Wang , X. Zhang , X. P. Wang , J. Li , K. Mochizuki , D. Kim , N. Kawakami , W. Yi , H. Obuse , B. C. Sanders , and P. Xue , Observation of topological edge states in parity−time-symmetric quantum walks, Nat. Phys. 13(11), 1117 (2017)

[36] S. V. Suchkov , A. A. Sukhorukov , J. Huang , S. V. Dmitriev , C. Lee , and Y. S. Kivshar , Nonlinear switching and solitons in PT-symmetric photonic systems, Laser Photonics Rev. 10(2), 177 (2016)

[37] S. Shen , M. R. Belić , Y. Zhang , Y. Li , T. Wang , Z. N. Tian , and Q. D. Chen , Edge solitons in non-Hermitian photonic lattices with type-II Dirac cones, Front. Phys. (Beijing) 20(4), 042203 (2025)

[38] S. Xu , J. Li , Y. Hou , J. R. He , Z. Fan , Y. Zhao , and L. Dong , Vortex light bullets in Rydberg atoms trapped in twisted PT-symmetric waveguide arrays, Phys. Rev. A 110(2), 023508 (2024)

[39] S. Gao , W. Peng , Q. Huang , Y. Zhao , H. Gong , and S. Xu , Vortex solitons in fractional Schrödinger equation with cold Rydberg atoms and symmetric potentials, Front. Phys. 21, 043201 (2026)

[40] C. Hang , G. Huang , and V. V. Konotop , PT symmetry with a system of three-level atoms, Phys. Rev. Lett. 110(8), 083604 (2013)

[41] J. Sheng , M.-A. Miri , D. N. Christodoulides , and M. Xiao , PT-symmetric optical potentials in a coherent atomic medium, Phys. Rev. A 88, 041803(R) (2013)

[42] Z. Zhang , Y. Zhang , J. Sheng , L. Yang , M. A. Miri , D. N. Christodoulides , B. He , Y. Zhang , and M. Xiao , Observation of parity−time symmetry in optically induced atomic lattices, Phys. Rev. Lett. 117(12), 123601 (2016)

[43] C. Hang and G. Huang , Parity−time symmetry with coherent atomic gases, Adv. Phys. X 2(3), 737 (2017)

[44] Z. Zhang , D. Ma , J. Sheng , Y. Zhang , and M. Xiao , Non-Hermitian optics in atomic systems, J. Phys. At. Mol. Opt. Phys. 51(7), 072001 (2018)

[45] L. Qin , C. Hang , and G. Huang , Controllable PT phase transition and asymmetric soliton scattering in atomic gases with linear and nonlinear potentials, Phys. Rev. A 99(4), 043832 (2019)

[46] F. K. Abdullaev , V. V. Konotop , M. Salerno , and A. V. Yulin , Dissipative periodic waves, solitons, and breathers of the nonlinear Schrödinger equation with complex potentials, Phys. Rev. E 82(5), 056606 (2010)

[47] F. K. Abdullaev , Y. V. Kartashov , V. V. Konotop , and D. A. Zezyulin , Solitons in PT-symmetric nonlinear lattices, Phys. Rev. A 83, 041805(R) (2011)

[48] D. A. Zezyulin , Y. V. Kartashov , and V. V. Konotop , Stability of solitons in PT-symmetric nonlinear potentials, Europhys. Lett. 96(6), 64003 (2011)

[49] Y. He , X. Zhu , D. Mihalache , J. Liu , and Z. Chen , Lattice solitons in PT-symmetric mixed linear-nonlinear optical lattices, Phys. Rev. A 85(1), 013831 (2012)

[50] C. Huang , C. Li , and L. Dong , Stabilization of multipolemode solitons in mixed linear−nonlinear lattices with a PT symmetry, Opt. Express 21(3), 3917 (2013)

[51] Y. Meng and Y. Liu , Power-dependent shaping of solitons in parity-time symmetric potentials with spatially modulated nonlinearity, J. Opt. Soc. Am. B 30(5), 1148 (2013)

[52] J. Xie , W. Chen , J. Lv , Z. Su , C. Yin , and Y. He , Nonlocal defect solitons in parity−time-symmetric photonic lattices with spatially modulated nonlinearity, J. Opt. Soc. Am. B 30(5), 1216 (2013)

[53] X. Zhu , Z. Shi , and H. Li , Gap solitons in parity-time symmetric mixed linear−nonlinear optical lattices with fourth-order diffraction, Opt. Commun. 382, 455 (2017)

[54] F. K. Abdullaev , V. A. Brazhnyi , and M. Salerno , Scattering of gap solitons by PT-symmetric defects, Phys. Rev. A 88(4), 043829 (2013)

[55] U. Al Khawaja , S. M. Al-Marzoug , H. Bahlouli , and Y. S. Kivshar , Unidirectional soliton flows in PT-symmetric potentials, Phys. Rev. A 88(2), 023830 (2013)

[56] S. M. Al-Marzoug , Scattering of solitons by complex PT symmetric gaussian potentials, Opt. Express 22(18), 22080 (2014)

[57] N. Karjanto , W. Hanif , B. A. Malomed , and H. Susanto , Interactions of bright and dark solitons with localizedsymmetric potentials, Chaos 25(2), 023112 (2015)

[58] L. Qin , C. Hang , Z. Shi , J. Qian , X. Feng , Y. Zhang , S. Xia , Z. Zhu , W. Liu , and X. Zhao , Soliton molecules and their scattering by a localized PT-symmetric potential in atomic gases, Opt. Express 31(7), 11116 (2023)

[59] H. Kang , L. Wen , and Y. Zhu , Normal or anomalous dispersion and gain in a resonant coherent medium, Phys. Rev. A 68(6), 063806 (2003)

[60] J. Yang , Nonlinear Waves in Integrable and Nonintegrable System, SIAM, Philadelphia, (2010)

[61] R. W. Boyd , Nonlinear Optics, 3rd Ed., Academic, Elsevier, (2008)

Rights & permissions

Higher Education Press