Shock wave generation and propagation in dissipative and nonlocal nonlinear Rydberg media

We investigate the generation of optical shock waves in strongly interacting Rydberg atomic gases with a spatially homogeneous dissipative potential. The Rydberg atom interaction induces an optical nonlocal nonlinarity. We focus on local nonlinear ($R_b\ll R_0$) and nonlocal nonlinear ($R_b\sim R_0$) regimes, where $R_b$ and $R_0$ are the characteristic length of the Rydberg nonlinearity and beam width, respectively. In the local regime, we show spatial width and contrast of the shock wave change monotonically when increasing strength of the dissipative potential and optical intensity. In the nonlocal regime, the characteristic quantity of the shock wave depend on $R_b/R_0$ and dissipative potential nontrivially and on the intensity monotonically. We find that formation of shock waves dominantly takes place when $R_b$ is smaller than $R_0$, while the propagation dynamics is largely linear when $R_b$ is comparable to or larger than $R_0$. Our results reveal nontrivial roles played by dissipation and nonlocality in the generation of shock waves, and provide a route to manipulate their profiles and stability. Our study furthermore opens new avenues to explore non-Hermitian physics, and nonlinear wave generation and propagation by controlling dissipation and nonlocality in the Rydberg media.

Recently, it has been shown that cold atomic gases interacting with laser fields provide a fertile ground for studying shock waves [10][11][12][13][14][15].When additionally coupling the light to highly excited Rydberg states [54,55], strong and long-range interactions between Rydberg atoms can be mapped to light fields through electromagnetically induced transparency (EIT) [56], generating strong nonlocal nonlinearities [57,58].The characteristic length of the nonlocal nonlinearity, determined by the blockade radius of the Rydberg gas, is in the or-der of micrometers, which is comparable to typical beam width.Using the strong Rydberg nonlocal nonlinearity (NNL), it has been shown that DSWs can be generated and manipulated in Rydberg atom gases [38].Dissipation plays an important roles in the study of Rydberg systems [59][60][61].In cold atom gases, dissipation, on the other hand, can be induced and controlled [28,30].This opens new opportunities for exploring shock waves in the interplay between the nonlocal nonlinearity and controllable dissipation that is otherwise difficult to achieve in other systems.
In this work, we study the generation and propagation of shock waves within a cold Rydberg atomic gas setting, incorporating an engineered, homogeneous dissipative potential that can be changed from loss to gain.This change is controlled by employing an incoherent pumping [30] and controlling the laser detuning [28].A nonlocal optical nonlinear interaction is induced by coupling low-lying electronic states to Rydberg S state via EIT [56,57].Depending on the blockade radius R b of the NNL and beam width R 0 , the system is in a local regime when R b ≪ R 0 , and nonlocal regime when R b ∼ R 0 .In the local regime, the nonlinear Schrödinger (NLS) equation governing the light propagation is cast into coupled Riemann equations.Formation of shock waves is signified by wave breaking.Excluding dissipation, wave breaking points are obtained analytically through the Riemann equations.Oscillation contrast [39,62] and spatial width of the shock wave depends on the strength of the dissipation and optical intensity monotonically in the local regime.In the nonlocal regime, the nonlocal degree of the optical nonlinearity modifies amplitudes and width of shock waves.Importantly properties of the shock wave exhibit complicated dependence on the nonlocality.We show that shock wave generation and propagation is important when R b < R 0 .When R b and R 0 are comparable, the medium is effectively linear, where the NNL becomes a homogeneous dispersive potential approximately.
The paper is arranged as follows.In Sec.II, we present the physical model that can lead to the dissipative and nonlocal nonlinear potential.The NLS equation that describes the propagation of the probe laser field is derived.In Sec.III, light propagation in the local regime is discussed.The impact of the dissipative potential on the oscillation contrast, width, and shock width is investigated.In Sec.IV, we explore the influence of the nonlocality on the generation and propagation of shock waves.Finally, conclusions are given in Sec.V.

A. Physical model
We consider a gas of cold atoms with an inverted Ytype four-level configuration [see Fig. 1(a)], where a weak probe laser field with half Rabi frequency Ω p couples the transitions |1⟩ ↔ |3⟩.A strong control and a dressed laser fields with half Rabi frequencies Ω c and Ω d couple the transition |2⟩ ↔ |3⟩ and |3⟩ ↔ |4⟩, correspondingly.Detuning ∆ α (α = 2, 3,4) gives difference between laser frequency and atomic transition.And Γ αβ are spontaneous emission decay rates from |β⟩ to |α⟩ (α < β).Here the incoherent decay from state |3⟩ causes loss of the probe field.To generate gain, an incoherent pumping (with pumping rate Γ 21 ) is used to pump atoms from |1⟩ to |2⟩.Driven by the control laser Ω c a small number of atoms are populated in state |3⟩ which provides a gain effect [30,63].
In this setting, state |4⟩ is a high-lying Rydberg state.The interaction between two Rydberg atoms respectively at positions r and r ′ is described by van der Waals potential V vdW ≡ ℏV (r ′ − r) = −ℏC 6 /|r ′ − r| 6 [64].When the light propagates in the medium [see Fig. 1(b)], Rydberg excitation in the vicinity of a Rydberg atom is strongly suppressed, due to the long-range Rydberg-Rydberg interaction.Such spatial dependent Rydberg blockade leads to nonlocal nonlinear optical interactions [65].Note that in the inverted Y-shaped excitation scheme shown in Fig. 1(a), the transition |1⟩ → |2⟩ → |3⟩ forms a Λ-type EIT, while |1⟩ → |3⟩ → |4⟩ forms a ladder type Rydberg-EIT.The interplay of the two paths can be controlled by the external lasers, giving rise to dissipative, nonlocal nonlinear interactions.
Under electric-dipole and rotating-wave approximations, Hamiltonian of the system is Ĥ = N a d 3 r Ĥ with Hamiltonian density Ĥ (ℏ ≡ 1), where N a is atomic density, and is the atomic transition operator between states |α⟩ and |β⟩.For weak excitation, Ŝαβ (r) are approximated by bosonic operators [66].Taking into account of decay, dynamics of the density matrix (matrix elements where Γ is the relaxation matrix describing the spontaneous emission and dephasing (see Appendix A).Propagation of the semiclassical probe field is governed by the Maxwell equation, where κ 13 = N a ω p |p 13 | 2 /(2ε 0 cℏ), with ω p the weak probe laser field of center frequency, p 13 the electric dipole matrix element associated with the transition |1⟩ ↔ |3⟩, ε 0 the vacuum dielectric constant, and c the speed of light in vacuum.In deriving the propagation equation, the paraxial and slowly varying envelope approximations have been applied.

B. Nonlinear envelope equation
For weak probe fields, the Maxwell and Bloch equations can be solved perturbatively.
The Rydberg-Rydberg interaction, on the other hand, is treated beyond the simple mean-field approximation [28,67].We then solve the Bloch equation up to the third-order of Ω p .This allows to derive a (2+1)D nonlocal nonlinear Schrödinger (NNLS) equation of the probe field [30,63,68], with r ⊥ = (x, y).Here V 1 is the dissipative potential that gives homogeneous gain or loss in the medium controlled by the laser parameters [30].Details on the control of the dissipative potential can be found in Appendix B. Nonlinear coefficient W characterizes the local Kerr nonlinearity (contributed by the weak, short-range interactions between photons and atoms [29]).G(r) is a nonlocal nonlinear response function characterizing respectively NNL [28,67,69].
Compared to the Rydberg induced nonlinearity, the conventional Kerr nonlinearity W is negligible.Using 87 Rb as an example with electronic states 0 W . Considering the probe beam radius R 0 = 5 µm and U 0 = 10 MHz, we obtain R b ≈ 1.94 µm and L diff ≈ 0.2 mm.The dimensionless strength of the local nonlinearity is W ≈ 0.01 ≪ 1, which is much smaller than the strength (in the order of 1) of the NNL, and will be neglected in the following.With this approximation and focusing on the diffraction along the x-axis, we convert the propagation equation in a dimensionless form, where we have defined dimensionless quantities The optical field and spatial coordinate have been scaled with respect to the maximal Rabi frequency U 0 and beam radius R 0 .The imaginary potential V = iV I , where V I is gain (loss) when V I > 0 (V I < 0).
Due to the Rydberg blockade, the response function has a soft-core shape [see Fig. 1(c)].However the expression of g(ξ ′ , ξ) is typically lengthy and complicated (see Appendix B for discussions).It can be approximated by an analytical form [30,38] 3) has a soft-core shape, which also agrees with the numerical data.Note that this system has a defocusing nonlinearity as g(ξ ′ , ξ) < 0, which is crucial for the generation of shock waves.
When R b is comparable to R 0 , the nonlinear interaction is nonlocal (i.e.σ is finite).In the opposite regime when R b ≪ R 0 , we have a local regime as σ ∼ 0. The nonlocality parameter σ can be varied by varying R b or R 0 .The blockade radius can be tuned by changing detuning, laser intensities, or choosing different Rydberg states as C 6 ∝ n 11 with n to be the principal quantum number.When σ ∼ 0, we can make a local field approximation, i.e.
Carrying out the spatial integration, we arrive at a propagation equation with local nonlinear interactions, where ḡ0 = −g 0 πσ (b 1 /b 2 ) 1/6 /(3b 1 ) is the effective interaction strength.The local interaction is similar to the conventional Kerr nonlinearity, though the nonlinearity is much stronger.In the following, we will discuss the local and nonlocal regime separately.

III. LOCAL RYDBERG NONLINEARITY REGIME
A. Euler-like fluid equation In the local regime, we start to investigate the generation of shock waves with the hydrodynamic approach.By treating the light field as a classical fluid, the hydrodynamic equation can be obtained by using the Madelung transformation u(ξ, ζ) = ρ(ξ, ζ)e iϕ(ξ,ζ) , Eq. ( 4) can be transformed into two Euler-like fluid equations, where ∂ξ 2 is quantum pressure, and v = ∂ϕ/∂ξ is the flow velocity of the light fluid.Neglecting Q for the moment, Eq. (5b) becomes Equations ( 5a) and ( 6) can be cast into the diagonal Riemann form where the Riemann invariant and hyperbolic speeds are defined by The light fluid intensity and the flow velocity are given by ρ = (r 1 − r 2 ) 2 /(4ḡ 0 ) and v = r 1 + r 2 .
We numerically solve the Riemann equation without dissipative potential with the initial condition Here ρ b and ρ h are the background and hump intensity, and ξ 0 is the width of the hump.In other words the left and right-moving components will form shock waves after propagating equal distances.

B. Breaking point
In the local regime and without dissipation, the breaking point can be found analytically.We first linearize Eq. ( 7) by means of the hodograph transform (see, Refs.[33,39,70]), which treats ξ and ζ as functions of r 1 and r 2 .The transformation yields, We then introduce two functions w 1 (r 1 , r 2 ) and w 2 (r 1 , r 2 ) such that Using the initial condition, w 1 and w 2 can be obtained, Wave breaking corresponds to the occurrence of a gradient catastrophe for which |∂r 1,2 /∂ξ| → ∞.As the right and left-moving wave are symmetric with respect to the ξ, we determine the breaking point of the right-moving branch explicitly.From Eq. ( 10), wave breaking occurs at a distance ζ such that [39] ζ = 2 3 One can evaluate ζ b approximately when the point of largest gradient in r 1 (ξ) lies in a region where r 2 ≈ √ ḡ0 ρ b .In this case, Eq. ( 12) becomes The shortest of distance ζ is reached close to the point ξ * for which |∂ρ/∂ξ| is maximal.We denote ξ * the coordinate of this point in the ξ direction and . This leads to the approximate relation, The breaking point as a function of ρ h is shown in the Fig. 2(d), which matches the numerical calculation well.The breaking point ζ b is reduced when increasing the hump intensity.Such relation is useful in controlling the generation of shock waves.For example, increasing the intensity of the hump peak allows for a shorter distance and faster visibility of the shock wave.Moreover, the breaking point ξ b along the ξ axis when the wave breaks along the ζ direction can be obtained [40], Here c s = √ ḡ0 ρ is the local sound speed.The results are shown in Fig. 2(e), which agrees with the numerical calculation well.
Including dissipation in the Riemann function, analytical solutions are in general not possible.Instead, we find the breaking point numerically.Breaking points for V I = 0.1 and −0.1 are shown in Fig. 2(d).It is found that the breaking point ζ b decrease as the dissipative potential changes from loss to gain.These results highlight the importance of dissipative potential on the generation of shock waves.For example, the gain potential (i.e.V I > 0) accelerates the generation of shock waves, as the wave breaks earlier.The loss potential (V I < 0) slows down their generation.
On the other hand, the nonlocality can also affect the breaking point.To be specific, the breaking point ζ b increases with σ.Therefore, a strong nonlocality will postpone the occurrence of shock waves [20].

C. Wave propagation
We now turn to investigate the propagation of shock waves by numerically solving Eq. ( 4), where the quantum pressure is taken into account explicitly.In Fig. 3 propagation of shock waves without external potential [Fig.3(a)], in loss potential [Fig.3(b)], and in gain potential [Fig.3(c)-(d)] are shown.Without dissipation, the initial hump splits into two density peaks firstly.When the shock wave forms, the wave front oscillates rapidly in the ξ direction, as depicted in the upper and middle panel of Fig. 3(a).The end of oscillations edge corresponds to small-amplitude edge of the shock wave [38,71].Before the shock wave reaches the boundary, the background field is not perturbed.
In a loss potential (V = −0.2i), the intensity of the background wave decays exponentially, as shown in Fig. 3(b).Shock waves form in the central region, characterized by the rapid oscillation at the shock edge.Amplitudes of the oscillation edge become smaller than that of the V = 0 case.Results of shock waves in a gain potential with V I = 0.1 are shown in Fig. 3(c).From the upper panel of Fig. 3(c), the intensity of the shock wave and background field, and the small-amplitude edge are all lager than the cases V I = 0 and V I = −0.2.This is a direct manifestation of the gain effect.On the other hand, further increasing the strength of the gain potential, the shock wave and background field quickly become unstable, causing catastrophic collapse [38].
In the presence of the dissipative potential, power of the field will decay or grow exponentially with the propagation distance ζ.V I = 0, the total, background as well as the exchange power is conserved as a function of ζ [lower panel of Fig. 3(a)].Their values are determined by the initial values.On the other hand, the power will grow (decay) exponentially when V I > 0 (V I < 0) when propagating in the medium, as shown in the lower panel of Fig. 3(b) and (c).

D. Contrast and shock width
As shown in Fig. 3, profiles of the shock wave exhibit a non-trivial dependence on the dissipative potential and initial state.Once the shock wave forms, rapid oscillations are found along the ξ axis.Including the quantum pressure, the gradient divergence of u is not available, which makes it impossible to calculate the breaking point.As shock waves oscillate rapidly, the maximal and minimal values of the oscillations provide a way to characterize the amplitude of the shock wave.Therefore we calculate the visibility of the oscillations near the soliton edge (i.e., the start of oscillation edge) of the shock wave by measuring the contrast [39,62] where ρ max and ρ min are the maximum and minimum values of ρ, as depicted in the left lower insert of Fig. 4(a).At a fixed propagation distance ζ, C as a function V I is shown in Fig. 4(a).As V I increases (from loss to gain), the contrast of the formed shock wave decreases.In Fig. 4(b), we show contrast C when varying ρ h .For a given V I , the contrast of the shock wave increases with increasing ρ h .Changing V I , such trend remains the same.These results show that we could enhance the contrast by using larger ρ h .We also calculate the shock wave width L 1 (measured from the center to the end of the intensity oscillation) [18,72], and the oscillation width L 2 (measured from the start to the end of the oscillation), as indicated in the insert of Fig. 4(a).We find L 1 increases with V I , as shown in Fig. 4(c).The reason is that the small-amplitude edge has a slight increase when the potential changes from loss to gain.However, the width of oscillation L 2 (orange solid line) decreases due to the soliton edge increases.
In Fig. 4(d) width L 1 and L 2 as function of the hump peak intensity ρ h are shown.The larger the hump density ρ h , the wider the width L 1 and L 2 .We can understand these distances by examining the local sound speed c s ∝ √ ρ. c s is a function of local density ρ that consists of both the background and hump density.Increasing the hump density will increase the local sound speed.This means L 1 will be larger with higher ρ h , after propagating certain ζ.The inner region defined by L 2 , on the other hand, describes propagation of solitons [38].Its front travels at the sound speed approximately.Hence L 2 increases when c s (ρ h ) is larger.

IV. NONLOCAL RYDBERG NONLINEARITY REGIME
We will consider the full soft-core potential using Eq.(3).To understand the role played by the nonlocality, we will solve Eq. ( 2) numerically by taking into the full soft-core potential Eq. (3).There are two different ways to change the soft-core potential Eq. (3).One can vary σ by fixing B 1 and B 2 , which requires to change parameters b 1 and b 2 correspondingly.The depth of the potential [see Fig. 1 , the depth of the potential increases rapidly as σ decreases, which eventually makes the numerical calculation unstable.Hence σ can not be too small in practice.To avoid this divergence, one can alternatively change σ while keeping the potential depth (hence b 1 ) constant.

A. Fix B1 and B2
To avoid numerical instabilities, we have chosen the smallest σ = 0.2 in the numerical simulation.When σ is small, the formation of shock waves are featured by conspicuous oscillations, as illustrated in Fig. 5(a).The behavior is similar to that of the local nonlinearity regime.Increasing σ, the oscillation frequency decreases.As a result, the contrast increases rapidly with increasing σ, and reaches its maximum value C ≈ 0.47 around σ ≈ 0.33, as shown in Fig 5(b).It subsequently decreases and saturates to a constant value.The width L 1 and L 2 depend on σ sensitively when σ is small.They decrease rapidly with increasing σ, as shown in Fig. 5(c)  and (d).This dependence can be understood by analyzing the sound speed.When B 1 and B 2 are fixed, we obtain c s = πg 0 ρB 1 /(3B 5/6 2 σ 5 ).The sound speed decreases rapidly (c s ∝ 1/ √ σ 5 ) when σ increases.Therefore both L 1 and L 2 reduce when σ is large.
When increasing σ, we in fact drive the response from highly nonlinear to a linear regime.In other words, when σ is small, the wave dynamics is strongly nonlinear, which promotes the generation of shock waves.By increasing σ, however, the response of the Rydberg medium becomes effectively linear.When σ (R b ) is large, the softcore potential is nearly a constant compared to the typical wavelength of the excitation.Assuming the nonlocal potential is a constant, one can carry out the integration in Eq. ( 2) and obtain a linear potential [72], i.e.
When increasing ρ h , the contrast increases monotonically, as shown in Fig. 6(a).Both L 1 and L 2 become larger for higher ρ h , as L 1 , L 2 ∝ c s ∝ √ ρ h .Similar dependence are also found in the local nonlinear case [see Fig. 4(b) and (d)].The difference is that both L 1 and L 2 are slightly larger in the nonlocal Rydberg medium than that of the local medium (for given ρ h ), mainly due to that the strength of the nonlocal interaction is different in the two figures.
We then examine the characteristic quantities as a function of V I numerically.The results are shown in Fig. 6(c) and (d).When varying V I , the contrast, L 1 and L 2 exhibit similar dependence on V I as found in the local nonlinearity case.Changing V I from −0.3 to 0.15, the contrast decreases slowly, as shown in Fig. 6(c), akin to the finding in the local regime, as shown in Fig. 4(a).Compared to the local case, shock width L 1 is barely change when increasing V I , as shown in Fig. 6(d).The weak dependence comes from the fact that the smallamplitude edge only has a slight change.Oscillation width L 2 decreases apparently, due to the soliton edge increases, similar to the local case shown in Fig. 4(c).The numerical data show that the contrast, L 1 and L 2 are all smaller than that of the local nonlinear case.

B. Fix the potential depth
In this case, b 1 is a constant.The sound speed c s ∝ √ σρ, which means the larger σ, the more separation between the left and right moving shock wave, as shown in Fig. 7(a).When σ = 0.5 and 1, we find high density peaks at in the inner region, which result from that the system is effectively linear.They also affect the contrast, shown in Fig. 7(b).It decreases gradually, arrives at a minimum, and then increases again with increasing σ.For large σ, the rising contrast is purely caused by the inner peaks and the background density, where the amplitude of the shock wave is marginal.Moreover, both L 1 and L 2 increase monotonically as we increase σ, due to c s ∝ √ σρ.Thus after propagating distance ζ, the leftand right-moving shock waves separately significantly.Dissipation, on the other hand, leads to a nearly global shifts to the contrast.As shown in Fig. 7(b), the contrast becomes larger when V I is negative (loss potential).When V I is positive, the contrast is only shifted lower slightly when σ < 0.3.The dissipation merely modifies L 1 as we increase σ, which shows the robustness of the shock wave propagation.On the contrary, non-negligible changes to L 2 are found when V I is nonzero.This means that the speed of the soliton varies apparently.Compared to results shown in Fig. 5(b)-(d), one notes that modifications of these quantities due to V I are identical in the two cases.When increasing ρ h , the contrast, L 1 and L 2 all increase [see Fig. 8(a) and (b)].This is because not only the sound speed, but also amplitudes of the oscillation increases with larger ρ h .Similar trends are also found in the previous case shown in Fig. 6(a) and (b).However, the contrast depends on V I non-trivially in the current case.By increasing V I , we find contrast has a minimal value C min [see Fig. 8(c)], which is different from the situation shown in Fig. 6(c).In the latter case, the contrast declines monotonically with increasing V I in the given parameter range.C min depends on not only V I , but also σ.We numerically obtain C min and the respective parameter V I and σ.In the inset of Fig. 8(c), the corresponding V I and σ are plotted.It shows that when σ increases, one has to decreases V I in order to find C min .Finally, L 1 and L 2 as a function of V I are shown in Fig. 8(d).The trend is similar to that of the previous case [Fig.6(d)].Their values are, however, larger in general.This results from the fact that the sound speed has different dependence on parameters in the two cases.

V. CONCLUSION
In this paper, we have elaborated a scheme that enables the generation and propagation of shock waves within an atomic gas involving a homogeneous dissipative potential and long-range Rydberg interaction under the condition of EIT.We have demonstrated that the homogeneous gain or loss potential significantly alters the power in the local nonlinearity regime.Both the oscillation contrast of shock waves and the oscillation width exhibit monotonic decreases as the potential changes from loss to gain.Furthermore, we have shown that in the NNL regime, the nonlocal degree of the nonlinearity mod-ifies the oscillation amplitude and width of shock waves.Additionally, we have illustrated the hump intensity of the initial state can enhance the visibility of shock waves in both local and NNL regimes.Our results reveal the nontrivial roles of dissipation and nonlocality in the generation of shock waves, providing new routes to manipulate their profiles and stability.Our study opens new avenues for exploring non-Hermitian dynamics [73][74][75][76], and nonlinear wave dynamics [77][78][79] modulated by the interplay between the NNL, and local and nonlocal [59,80] dissipation in highly controllable Rydberg gases.44 where From the zeroth order solution, we find that the incoherent population pumping rate Γ 21 is a key parameter in the zeroth order solution.If Γ 21 = 0, all populations are in the ground state, i.e., ρ (0) 11 = 1, and other state population and coherence are both zero, ρ (0) αβ = 0.However, when Γ 21 ̸ = 0, we have ρ (0) 33 ̸ = 0, and hence a gain to the probe field will be realized when the probe field is coupled to the states |1⟩ and |3⟩.
First-order solutions: At the first order, the solutions of ρ Second-order solutions: At the second order, the matrix elements can be solved by the equation Solving the matrix equations above yields ρ where 31 at the third order is obtained from Eq. (B4) where a where should be solved firstly in order to solve the nonlocal nonlinear Schrödinger equation above.
In dimensionless NNLS equation Eq. ( 2), the linear potential reads (B7) Here, V R and V I are the real and imaginary parts of linear potential, respectively.V R is a constant and independent on the intensity of the shock wave, and is indeed negligible.V I < 0 show a loss to the probe field due to the dissipation of system when ρ (0) 33 = 0 in the zeroth order for Γ 21 = 0.With increasing of incoherent pumping Γ 21 , ρ (0) 33 ̸ = 0 and a gain to the probe field will be realized.As results show in Fig. 9, a loss/gain potential V I can be realized by adjusting the Γ 21 .When Γ 21 < 2.45 MHz, the potential is a loss.Otherwise, the potential is a gain.

Solutions of two-body density-matrix elements
We solve the second order solution of the correlators of ρ (2) 41,41 which can be obtained from  (B8) The solution for ρ
and |4⟩ = |nS 1/2 ⟩, and typical parameters ∆ 2 = 20 MHz, ∆ 3 = −200 MHz, ∆ 4 = 10 MHz, Γ 3 = 2π × 6.1 MHz, Γ 4 = 2π × 16.7 kHz, Ω c = 40 MHz, Ω d = 10 MHz, and N a = 2.3 × 10 10 cm −3 , we evaluate the local nonlinearity, characterized by a dimensionless strength W = 2L diff U 2 where σ = R b /R 0 characterizes the nonlocal degree of the nonlinearity.Here R b = |C 6 /δ EIT | 1/6 is the radius of the blockade sphere with δ EIT = |Ω c | 2 /∆ 3 the linewidth of EIT transition spectrum.B 1 and B 2 are the coefficients determined by laser parameters.The relation between B 1,2 and b 1,2 are B 1 = σ 6 /b 2 and B 2 = b 1 /b 2 .As shown in Fig. 1(c), Eq. ( Figure 2(a) and (b) are the propagation of the right and left-moving Riemann waves, respectively.The wave is stable before the shock onsets.After a critical distance (marked by stars on the figure) the shock forms.A distinctive feature is that the wave becomes steepening at the critical distance.These points are the shock wave breaking points (see more discussions in Sec.III B).These points are symmetric for the left and right-moving solutions, as depicted in Fig. 2(c).

FIG. 2 .
FIG. 2. (a) Right-moving and (b) left-moving Riemann waves.The stars represent the breaking position of the wave, after which the shock wave forms.(c) Riemann waves r1 and r2 as function as ξ for ζ = 0, 0.47, 0.94.The stars represent the breaking position.Upper panel: The right-moving part r1.Lower panel: The left-moving part r2.(d) Breaking point in the ζ direction (ζ b ) as function as the hump peak intensity ρ h .The star marks the breaking point corresponding to panel (a) and (b).(e) The same as (d), but for breaking point in the ξ direction (ξ b ).Other parameters are ρ b = 1, ρ h = 2, ξ0 = 1, ḡ0 = −1.44,and VI = 0.In panel (d), we additionally show data with dissipative potentials.

FIG. 4 .
FIG. 4. (a) Oscillation contrast versus VI .The insert illustrate the maximum and minimum values of ρ, i.e., ρmax and ρmin, shock width L1, and oscillation width L2.The parameters ρ b = 1, ρ h = 2, ξ0 = 1, ḡ0 = −1.44,and ζ = 3.(b) Contrast versus ρ h with VI = 0, 0.1, and −0.1, respectively.(c) Shock width L1, measured from the center to the end of oscillations, with respect to the intensity of imaginary potential VI (blue solid line).The oscillation width L2, measured from the start to the end of oscillations (orange dotted-dashed line).(d) L1 and L2 by varying the hump peak intensity ρ h with VI = 0.The other parameters same as panel (a).

FIG. 9 .
FIG. 9. (a) The imaginary potential as function as incoherent pumping Γ21.When Γ21 < 2.45 MHz, the potential is a loss.Otherwise, the potential is a gain.The red star represents VI = 0 when Γ21 = 2.45 MHz.The lower green and upper orange regions mark the loss and gain potential, respectively.

Ĥ1H 1 H 1 H 1 ξH 1 H 1 H
FIG. 10.(a) Wave propagation when σ = 0.5.The dashed line marks the trajectory of the sound wave.It is clear that the wave propagates at the speed of sound cs, indicating the system is in the linear regime.Other parameters are same as in Fig. 5(a), corresponding to the linear propagation discussed in Sec.IV A. (b) The same as (a).Other parameters are same as in Fig. 7(a), corresponding to the linear propagation discussed in Sec.IV B.