Formation of dispersive shock waves in a saturable nonlinear medium
Sergey K. Ivanov, Jules-Elemir Suchorski, Anatoly M. Kamchatnov, Mathieu Isoard, Nicolas Pavloff
FFormation of dispersive shock waves in a saturable nonlinear medium
Sergey K. Ivanov,
1, 2
Jules-El´emir Suchorski, Anatoly M. Kamchatnov,
1, 2
Mathieu Isoard, and Nicolas Pavloff Moscow Institute of Physics and Technology, Institutsky lane 9, Dolgoprudny, Moscow region, 141700, Russia Institute of Spectroscopy, Russian Academy of Sciences, Troitsk, Moscow, 108840, Russia Universit´e Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France
We use the Gurevich-Pitaevskii approach based on the Whitham averaging method for studyingthe formation of dispersive shock waves in an intense light pulse propagating through a saturablenonlinear medium. Although the Whitham modulation equations cannot be diagonalized in thiscase, the main characteristics of the dispersive shock can be derived by means of an analysis of theproperties of these equations at the boundaries of the shock. Our approach generalizes a previ-ous analysis of step-like initial intensity distributions to a more realistic type of initial light pulseand makes it possible to determine, in a setting of experimental interest, the value of measurablequantities such as the wave-breaking time or the position and light intensity of the shock edges.
1. INTRODUCTION
The propagation of nonlinear waves in dispersive mediahas attracted much attention in various fields of researchsuch as water waves, plasma physics, nonlinear optics,Bose-Einstein condensates and others. In particular, theexpansion of an initial state with a fairly smooth andlarge profile is accompanied by gradual steepening fol-lowed by wave breaking resulting in the formation of adispersive shock wave (DSW), and this phenomenon wasexperimentally observed, for example, in Bose-Einsteincondensates [1–3] and nonlinear optics [4–9]. In these sys-tems, the dynamics of the pulse is well described by thenonlinear Schr¨odinger (NLS) or Gross-Pitaevskii equa-tion and the initial stage of evolution admits a purelyhydrodynamic description in terms of the classical Rie-mann method, see, e.g., Refs. [10–14]. After the wavebreaking moment, when DSWs are formed, the evolutionof the wave structure can be described by the Whithammethod [15–18]. In the case of cubic Kerr-like nonlin-earity the NLS equation is completely integrable, hencethe Whitham modulation equations can be put in a di-agonal Riemann form [19, 20], and, with the use of theGurevich-Pitaevskii approach [21], it has been possibleto develop a detailed analytic theory of the evolution ofDSWs [12, 22–24]. This produces an excellent descriptionof experiments [6, 7] on evolution of initial discontinuitiesin the intensity distribution of optical pulses propagatingin optical fibers.As first understood by Sagdeev [25] in the context ofviscous stationary shocks whose width is much greaterthat a typical soliton width, the formation of DSWs isa universal phenomenon which occurs in a number ofsystems demonstrating nonlinear dispersive waves. Inparticular, the formation of DSWs has been observed inthe propagation of light beams through a photorefrac-tive medium [26, 27]. However, in such a system thenonlinearity is not of the Kerr type, and the theory ofRefs. [12, 22–24] thus needs to be modified. Of course,the validity of the Whitham averaging approach for thedescription of the modulations of a nonlinear wave doesnot depend on the complete integrability of the wave equation. Nevertheless, it is difficult to put it into prac-tice for non-completely integrable equations because ofthe lack of diagonalized form of the Whitham modula-tion equations.The first general statement about the properties ofDSWs applicable to a non-integrable situation was madeby Gurevich and Meshcherkin in Ref. [28]. In this workthe authors claimed that, when a DSW is formed afterwave breaking of a “simple wave”, that is of a wave forwhich one of the non-dispersive Riemann invariants isconstant, then the value of this constant Riemann invari-ant remains the same at both extremities of the DSW. Inother words, this means that such a wave breaking leadsto the formation of a single DSW, at variance with the sit-uation with viscous shocks where more complicated wavestructures are generated (see, e.g., Ref. [29]). The nextimportant step was made by El [30] who showed that, insituations of the Gurevich-Meshcherkin type, when theWhitham equations at the small-amplitude edge includethe linear “number of waves” conservation law, this equa-tion can be reduced to an ordinary differential equationwhose solution provides the wave number of a linear waveat this edge. Under some reservations, a similar approachcan be developed for the soliton edge of the DSW. As aresult, one can calculate the speeds of both edges (soli-tonic, large amplitude one, and linear, small amplitudeone) in the important case of a step-like initial condition.This approach was applied to many concrete physical sit-uations [30–38], and was recently extended in Ref. [39] toa general case of a simple-wave breaking and then appliedto shallow water waves described by the Serre equation[40]. In the present paper, we use the same approachto study the propagation of optical pulses and beamsthrough a saturable nonlinear medium. This makes itpossible to considerably extends the theory of Ref. [32]and provides a more realistic explanation of the results ofRef. [27] and of future experimental studies in nonlinearoptics [41]. a r X i v : . [ n li n . PS ] J u l
2. FORMULATION OF THE PROBLEM
The formation of DSWs has been observed in the spa-tial evolution of light beams propagating through self-defocusing photorefractive crystals in Refs. [26, 27]. Ini-tial non-uniformities of the beam give rise to breaking sin-gularities resulting in the formation of dispersive shocks.As is well known [42], in the paraxial approximation, thepropagation of the complex amplitude A = A ( X, Y, Z ) ofthe electric field of a monochromatic beam is describedby the equation i ∂A∂Z + 12 k n ∆ ⊥ A + k δn ( I ) A = 0 , (1)where k = 2 π/λ is the carrier wave number, Z is thecoordinate along the beam, X , Y are transverse coordi-nates, ∆ ⊥ = ∂ /∂ X + ∂ /∂ Y is the transverse Lapla-cian, n is the linear refractive index, and δn is a nonlin-ear index which depends on the light intensity I = | A | .It is often the case that the nonlinearity is not of pureKerr type (i.e., not exactly proportional to I ) but satu-rates at large intensity. We consider here the case of adefocusing saturable medium where δn is of the form δn = − n II + I sat , (2)where n is a constant positive coefficient. This situa-tion is encountered for instance in semiconductor dopedglasses [43] and in photorefractive media [44]. In this lat-ter case n linearly depends on the applied electric fieldand on the electro-optical index. A near-resonant laserfield propagating inside a hot atomic vapor is also de-scribed by a saturable mixed absorptive and dispersivesusceptibility [45, 46]. At negative detuning the nonlin-earity is defocusing, and if the detuning is large enough,absorption effects are small compared to nonlinear ones.As a result, the propagation of the beam is described byan envelope equation which can be cast in the form ofEq. (1), with a nonlinearity of the type (2), see e.g., Ref.[47, 48].We define dimensionless units by choosing a refer-ence intensity I ref which can be chosen for instance asthe background intensity in the situation considered inRefs. [26, 27] where the Z = 0 light distribution has theform of a region of decreased [26] or increased [27] lightintensity perturbating a uniform background. Anothernatural choice would be I ref = I sat . We define dimen-sionless variables t = k n I ref I sat Z, x = k (cid:114) n n I ref I sat X,y = k (cid:114) n n I ref I sat Y, ψ = A √ I ref . (3)We consider a geometry where the transverse profile istranslationally invariant and depends on a single coordi- nate x . Then, the dimensionless generalized NLS equa-tion (1) takes the form iψ t + 12 ψ xx − | ψ | γ | ψ | ψ = 0 , (4)where γ = I ref /I sat . When the light intensity is smallcompared to the saturation intensity, i.e., when γ | ψ | (cid:28)
1, Eq. (4) reduces to the usual defocusing NLS equation;but at large intensity the nonlinearity saturates.The transition from the function ψ to the dimension-less intensity ρ and chirp u , is performed by means of theMadelung transform ψ ( x, t ) = (cid:112) ρ ( x, t ) exp (cid:18) i (cid:90) x u ( x (cid:48) , t ) dx (cid:48) (cid:19) . (5)After substitution of this expression into Eq. (4), separa-tion of the real and imaginary parts and differentiation ofone of the equations with respect to x , we get the system ρ t + ( ρu ) x = 0 ,u t + uu x + ρ x (1 + γρ ) + (cid:18) ρ x ρ − ρ xx ρ (cid:19) x = 0 . (6)In the problem of light beam propagation, the last termof the left-hand side of the second equation accounts fordiffractive effects. It is sometimes referred to as a “disper-sive term”, or “quantum pressure” or “quantum poten-tial” depending on the context. The Madelung transformreveals that light propagating in a nonlinear medium be-haves as a fluid, amenable to a hydrodynamic treatment,see Ref. [49]. The light intensity ρ in the hydrodynamicformulation (6) is interpreted as the density of this ef-fective fluid and the chirp u is the effective flow velocity.The coordinate t along the beam plays the role of time,and in this case the last term in (6) should be referred toas “dispersive”.The present paper in organized in two main parts. Inthe first one (Sec. 3) we describe the dispersionless evo-lution of a light pulse in the presence of a background.After a certain distance, the wave breaks and a DSW isformed. The main characteristics of this shock wave arestudied in the second part of the work (Sec. 4).We consider the case where the initial profile has theform of an increased parabolic intensity bump over a con-stant and stationary background (with u ( x, t = 0) = 0).If this profile is smooth enough (in a way which will bequantified in Sec. 3.2) the separation of the initial bumpin two counter propagating pulses can be described by adispersionless approach: this is performed in Sec. 3. Thedifficulty lies in the fact that the initial profile is not asimple wave in a hydrodynamic sense of this terminol-ogy. This means that the regions where the initial non-dispersive Riemann invariants depend on position signif-icantly overlap. Therefore, to study this problem, weshould resort to the Riemann hodograph method [10–14]which we specialize for the current photorefractive sys-tem in Sec. 3.1. Since the initial intensity profile has adiscontinuity in its first derivative, then in the process ofevolution, two simple rarefaction waves are formed alongthe edges (they are described in Sec. 3.3), and the cen-tral part is the region where the hodograph transformshould be employed. Due to nonlinearity, the profiles atboth extremities of the pulse gradually steepen and thisleads to wave breaking after some finite sample length,or equivalently some finite “time” which we denote asthe wave breaking time, t WB . Typically this occurs in aregion where only one Riemann invariant varies, and thecorresponding DSW results from the breaking of a simplewave, which is the case considered in the second part ofthe present paper.After the wave breaking, the dispersive effects can-not longer be neglected. For their description, we re-sort in Sec. 4 to Whitham modulation theory [15, 16]which is based on the large difference in spatial andtemporal scales between the rapid nonlinear oscillationsand their slow envelope. However, Eq. (4) being non-integrable, the Whitham modulation equations cannotbe transformed to a diagonal Riemann form: the lackof Riemann invariants hinders the application of the fullWhitham modulation theory to our system after wavebreaking. Nonetheless, the method of Refs. [30, 39] isbased on the universal applicability of Whitham’s “num-ber of waves conservation law” as well as on the conjec-ture of the applicability of its soliton counterpart to theabove mentioned class of initial conditions. Such an ap-proach is substantiated by comparison with similar situa-tions in the case of completely integrable wave equations.It makes it possible to calculate the limiting character-istic velocities of the Whitham modulation equations atthe boundary with the smooth part of the pulse whoseevolution obeys the dispersionless approximation equa-tions, even after the wave breaking time t WB . We willtreat in two separate subsections the case of a positive(Sec. 4.1) and of a negative (Sec. 4.2) initial intensitypulse.Having formulated the problem, we now proceed to thestudy of the dispersionless stage of evolution.
3. DISPERSIONLESS STAGE
For smooth enough wave patterns we can neglect thedispersive term in the second equation of the system (6)and the initial evolution of the system is described by theso-called dispersionless equations which can be put in aform equivalent to the equations of inviscid gas dynamics ρ t + ( ρu ) x = 0 , u t + uu x + c ( ρ ) ρ ρ x = 0 , (7)where c ( ρ ) = √ ρ γρ (8)is the local sound velocity in the medium. These equa-tions can be cast into a diagonal form by introducing new variables, known as the Riemann invariants r ± ( x, t ) = u ( x, t )2 ± (cid:90) ρ ( x,t )0 c ( ρ (cid:48) ) ρ (cid:48) dρ (cid:48) = u ( x, t )2 ± √ γ arctan (cid:112) γρ ( x, t ) , (9)whose evolution is described by the following equations ∂r + ∂t + v + ( r + , r − ) ∂r + ∂x = 0 , (10a) ∂r − ∂t + v − ( r + , r − ) ∂r − ∂x = 0 , (10b)with Riemann velocities v ± = u ± c. (11)In this last equation, the dependence of u and c on r + and r − is obtained from Eqs. (9). This yields u = r + + r − , (12)whereas c is computed as follows: one inverts the relation r + − r − = (cid:90) ρ c ( ρ (cid:48) ) ρ (cid:48) dρ (cid:48) , (13)which yields ρ = 1 γ tan (cid:20) √ γ r + − r − ) (cid:21) , (14)and then one evaluates c ( ρ ( r + − r − )) from (8). Once thesolution of Eqs. (10) are found and the Riemann invari-ants are known, then the physical variables are obtainedfrom relations (12) and (14). The Riemann method which we presenting now, en-ables one to find the solutions of Eqs. (10) in the genericregion where both Riemann invariants are changing. Rie-mann noticed that Eqs. (10) become linear with respectto the independent variables x and t if these are con-sidered as functions of the Riemann invariants: x = x ( r + , r − ), t = t ( r + , r − ). After this “hodograph trans-form” from the real space ( x, t ) to the hodograph plane( r + , r − ) we arrive at the system ∂x∂r − − v + ( r + , r − ) ∂t∂r − = 0 ,∂x∂r + − v − ( r + , r − ) ∂t∂r + = 0 . (15)It should be noticed that the Jacobian of this transfor-mation is equal to J = (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( x, t ) ∂ ( r + , r − ) (cid:12)(cid:12)(cid:12)(cid:12) = ∂t∂r + ∂t∂r − ( v − − v + ) , (16)thus the hodograph transform breaks down whenever ∂t/∂r + = 0 (or ∂t/∂r − = 0), which, by virtue ofEqs. (15), implies ∂x/∂r + = 0 (or ∂x/∂r − = 0). Thismeans that this approach makes sense only in a regionwhere both Riemann invariants are position-dependent.We look for the solution of the system (15) in the form x − v + ( r + , r − ) t = w + ( r + , r − ) , (17a) x − v − ( r + , r − ) t = w − ( r + , r − ) . (17b)Inserting the above expressions in the system (15) yields,with account of Eqs. (11),1 w + − w − ∂w + ∂r − = 1 v + − v − ∂v + ∂r − , w + − w − ∂w − ∂r + = 1 v + − v − ∂v − ∂r + . (18)From expressions (11), (12) and (14) one sees that ∂v + /∂r − = ∂v − /∂r + . This implies that ∂w + /∂r − = ∂w − /∂r + ; we can thus represent w ± in term of a type of“potential” function W ( r + , r − ): w + = ∂W∂r + , w − = ∂W∂r − . (19)Inserting the expressions (19) into the system (18) shows,with account of Eqs. (11), (12) and (14), that W is solu-tion of the Euler-Poisson equation ∂ W∂r + ∂r − + a ( r + , r − ) ∂W∂r + + b ( r + , r − ) ∂W∂r − = 0 , (20)with a ( r + , r − ) = − b ( r + , r − ) == − v + − v − ∂v + ∂r − = − − c (cid:48) ( r + − r − )2 c ( r + − r − ) , (21)where c (cid:48) ( r ) = dc/dr and c is computed as a function of r = r + − r − from expressions (8) and (14).The characteristics of the second order partial differ-ential equation (20) are the straight lines r + = ξ = constand r − = η = const, parallel to the coordinates axisin the hodograph plane. The Riemann method is basedon the idea that one can find the solution of the Euler-Poisson equation in a form similar to d’Alembert solutionof the wave equation, with explicit account of the initialconditions which fix the value of W on some curve C inthe hodograph plane. These initial data are transferredalong the characteristics into the domain of interest, sothat the function W can be found at any point P = ( ξ, η ).Riemann showed (see, e.g., Refs. [50, 51]) that W ( P )can be represented in the form W ( P ) = 12 ( RW ) A + 12 ( RW ) B − (cid:90) BA ( V dr + + U dr − ) , (22) where the points A and B are projections of the “obser-vation” point P onto C along the r + and r − axis respec-tively. The integral in (22) is performed along C , with U = 12 (cid:18) R ∂W∂r − − W ∂R∂r − (cid:19) + aW R,V = 12 (cid:18) W ∂R∂r + − R ∂W∂r + (cid:19) − bW R. (23) R ( r + , r − ; ξ, η ) is the “Riemann function” which satisfiesthe equation ∂ R∂r + ∂r − − a ∂R∂r + − b ∂R∂r − − (cid:18) ∂a∂r + + ∂b∂r − (cid:19) R = 0 , (24)with the additional conditions: ∂R∂r + − bR =0 along the characteristic r − = η,∂R∂r − − aR =0 along the characteristic r + = ξ, (25)and R ( ξ, η ; ξ, η ) = 1.One might think, looking at Eq. (24), that we didnot progress towards the determination of the solutionto Eq. (20). But we have replaced the initial conditionsfor Eq. (20) by standard boundary conditions (25) forEq. (24), independent of the initial values of ρ and u . Theknowledge of the initial properties of the flow is encap-sulated in the known value of W along C , which appearsin the right hand side of Eq. (22). We now consider a specific type of initial condition, forwhich u ( x,
0) = 0, with an initial parabolic bump profile ρ ( x,
0) = ρ ( x ) with ρ ( x ) = ρ + ( ρ m − ρ ) (cid:18) − x l (cid:19) , | x | ≤ l,ρ , | x | > l. (26)Here ρ is the background intensity, l and ρ m are thewidth and maximal intensity of the initial bump. Thisinitial density profile is represented in the left panel onthe top row of Fig. 1, it is an even function of x ; gen-eralization to non-symmetric distributions is straightfor-ward. If l (cid:29) ρ m − ρ , then deviations of the exact solu-tion from its hydrodynamic approximation are negligiblysmall almost everywhere, except in small regions at theboundaries of the bump.We denote the characteristic values of the Riemanninvariant as r and r m with r = r + ( | x | ≥ l, t = 0) = 1 √ γ arctan √ γρ , (27a) r m = r + ( x = 0 , t = 0) = 1 √ γ arctan √ γρ m , (27b) − −
20 0 20 400 . . ρ ρ m t = 01 3 x ρ − −
20 0 20 400 . . t = 10 ρ I IIIx ρ − −
20 0 20 400 . . t = 40 ρ I II l II r IIIx ρ − −
50 0 50 1000 . . t = 100 ρ I II l II r IIIx ρ . . r r m − −
20 0 20 40 − . − . − − r − r m x r ± . . r r m − −
20 0 20 40 − . − . − − r − r m x r ± . . r r m − −
20 0 20 40 − . − . − − r − r m x r ± . . r r m − −
50 0 50 100 − . − . − − r − r m x r ± . . . . − − . − . − . − . r − r − r m r m r + r − . . . . − − . − . − . − . r − r − r m r m I IIIr + r − . . . . − − . − . − . − . r − r − r m r m III l II r IIIr + r − . . . . − − . − . − . − . r − r − r m r m I , II l III , II r r + r − Fig. 1 . Behavior of characteristic quantities of the system. Each column corresponds to a given value of t , each row to adifferent type of quantity. The top row displays the light intensity distribution ρ ( x, t ) plotted as a function of x . Both thenumerical solution of Eq. (4) (red thick curves) and the analytic dispersionless solution (blue curves) are shown. The initialstate is represented in the left panel and corresponds to Eq. (26) with ρ = 0 . ρ m = 2, l = 20. The dynamics of the system isgoverned by Eq. (4) with γ = 1. The middle row displays the (analytic result for) the Riemann invariants r + ( x, t ) ( r + ∈ [ r , r m ])and r − ( x, t ) ( r − ∈ [ − r m , − r ]) plotted as functions of x . The set of initial condition corresponds to r = 0 .
651 and r m = 0 . r , r m ] × [ − r m , − r ]. Arabicand Roman numbers correspond to the notations of the regions from the top rows of the figure. Here the red curves correspondto t = 0, the blue curves to t = 10, the orange curves to t = 40 and the green curves to t = 100. The red dotted line representsthe diagonal coinciding with the initial curve C . see the left panel in the middle row of Fig. 1. It is con-venient to denote by x ( r ) the positive branch of the re-ciprocal function of r + ( x, t = 0): r + ( x,
0) = (cid:90) ρ c ( ρ (cid:48) ) ρ (cid:48) dρ (cid:48) . (28)One gets ρ ( r + ) = 1 γ tan ( √ γr + ) , (29)and x ( r ) = l (cid:115) ρ m − ρ ( r ) ρ m − ρ . (30) For the zero velocity initial profile we consider, Eq. (12)shows that the curve C which represents the initial condi-tion ( r + ( x, , r − ( x, r − = − r + = − r . Along itEqs. (17) with t = 0 give ∂W∂r − (cid:12)(cid:12)(cid:12)(cid:12) C = ∂W∂r + (cid:12)(cid:12)(cid:12)(cid:12) C = x. (31)hence W is constant along C . The value of this constantcan be arbitrarily fixed to zero: W ( r, − r ) = 0, and ex-pressions (23) reduce to U = x R ( r, − r ; ξ, η ) , V = − x R ( r, − r ; ξ, η ) . (32)The top row of Fig. 1 shows the light intensity ρ ( x, t )plotted as a function of position for different t . The mid-dle row represents the corresponding distributions of theRiemann invariants. At a given time, the x axis can beconsidered as divided into several domains, each requir-ing a specific treatment. Each domain is characterized bythe behavior of the Riemann invariants and is identifiedin the upper row of Fig. 1. The domains in which bothRiemann invariants depend on x are labeled by Arabicnumbers and the ones in which only one Riemann invari-ant depends on x are labeled by Roman numbers. As wasnoted above, at the initial moment r + = − r − , as can beseen, e.g., in the left panel of the middle row of Fig. 1.Then, one of the Riemann invariants begins to move inthe positive direction of the x axis, and the other movesin the opposite direction. This behavior initially leadsto the configuration represented in the second column ofFig. 1, where two simple-wave regions ( I and III ) and anew region (labeled “2”) have appeared. At later times(i.e., for longer sample lengths) region 2 persists whileregions 1 and 3 vanish and new simple-wave regions II l and II r appear; this is illustrated in the third columnsof Fig. 1. At even larger lengths (right column of Fig. 1),region 2 also disappears and only simple-wave regionspersist: The initial pulse has completely split into twopulses propagating in opposite directions.It is worth noticing that the wave breaking correspondsto an overlap between different regions which results ina multi-valued solution. If we consider for instance theright propagating pulse, at the wave breaking time t WB ,region III starts overlapping with the (unnamed) quies-cent region where both Riemann invariants are constants,at the right of the plot. Then these two regions bothoverlap with region II r . These two configurations are re-spectively illustrated in the columns t = 40 and t = 100of Fig. 1. We will not dwell on this aspect in detail, sincea multi-valued solution is nonphysical and, when disper-sion is taken into account, it is replaced by a dispersiveshock wave (considered in Sec. 4).The bottom row in Fig. 1 represents the behavior of r + and r − in the hodograph plane. Here, as before, thesimple-wave regions are indicated by Roman numeralswhile the Arabic numerals correspond to regions whereboth Riemann invariants depend on x . In each of thethree domains 1, 2, and 3, the solution W of the Euler-Poisson equation has a different expression. In order todescribe these three branches, following Ludford [52], weintroduce several sheets in the characteristic plane byunfolding the domain [ r , r m ] × [ − r m , − r ] into a fourtimes larger region as illustrated in Fig. 2. The poten-tial W ( r + , r − ) can now take a different form in each ofthe regions labeled 1, 2, and 3 in Fig. 2 and still be con-sidered as a single valued. From the relation (22) wecan obtain W ( ξ, η ) in regions 1 and 3, by computingthe right hand side integral along the anti-diagonal C ,between the points A of coordinates ( − η, η ) and B of r − r − r m r m r − r I III II l II r II l II r −−−−−−−−→ r + increases −−−−−−−−→ r + decreases r − i n c r e a s e s ← −−−−−−−− r − d ec r e a s e s ← −−−−−−−− t = 0 t = 10 t = 40 t = 100 Fig. 2 . Distributions of r + and r − in the four-sheeted hodo-graph plane (see the text). The curves are the unfolded ver-sions of the ones of the bottom row of Fig. 1. The red linecorresponds to the input distribution ( t = 0), while othercolors correspond to other values of t , with the same colorcode as in the bottom row of Fig. 1. Curved arrows indicatethe direction of unfolding of the domain [ r , r m ] × [ − r m , − r ].The whole region above red line is unreachable for the initialdistribution which we consider. coordinates ( ξ, − ξ ): W (3) ( ξ, η ) = − W (1) ( ξ, η ) = (cid:90) ξ − η x ( r ) R ( r, − r ; ξ, η ) dr (33)The difference in signs in the above expressions of W comes from the fact that x = ∓ x ( r ) depending on if one isin region 1 or 3. For the region 2, using unfolded surfaces(see Fig. 2 and also Ref. [14]) and upon integrating byparts one obtains W (2) ( ξ, η ) = (cid:16) RW (1) (cid:17) B + (cid:16) RW (3) (cid:17) A + (cid:90) CA (cid:18) ∂R∂r − − aR (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) r + = r m W (3) dr − − (cid:90) BC (cid:18) ∂R∂r + − bR (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) r − = − r m W (1) dr + , (34)where the coordinates of the relevant points are: A =( r m , η ), B = ( ξ, − r m ) and C = ( r m , − r m ).From the conditions (25) one can find that R ( r + , η ; ξ, η ) = (cid:115) c ( ξ − η ) c ( r + − η ) exp (cid:18) (cid:90) r + − ηξ − η drc ( r ) (cid:19) ,R ( ξ, r − ; ξ, η ) = (cid:115) c ( ξ − η ) c ( ξ − r − ) exp (cid:32) (cid:90) ξ − r − ξ − η drc ( r ) (cid:33) . (35)These expressions suggest that R can be looked for in theform R ( r + , r − ; ξ, η ) = (cid:115) c ( ξ − η ) c ( r + − r − ) exp (cid:18) (cid:90) r + − r − ξ − η drc ( r ) (cid:19) F ( r + , r − ; ξ, η ) , (36)where F ( r + , η ; ξ, η ) = F ( ξ, r − ; ξ, η ) = 1. For smallenough t , when ξ is close to r m and η is close to − r m ,the integrand functions in Eq. (34) are small by virtueof Eqs. (25); accordingly F (cid:39) γ = 0. Thus, using this approximation, we obtainan approximate expression for the Riemann function (seealso Ref. [14]): R ( r + , r − ; ξ, η ) (cid:39) R ( r + − r − , ξ − η ) , (37)where R ( r , r ) = (cid:115) c ( r ) c ( r ) exp (cid:18) (cid:90) r r drc ( r ) (cid:19) = (cid:115) c ( r ) ρ ( r ) c ( r ) ρ ( r ) , (38)where ρ is expressed as a function of r (or r ) throughEq. (14). A simple calculation yields R ( r , r ) = (cid:118)(cid:117)(cid:117)(cid:117)(cid:116) sin (cid:16) √ γ r (cid:17) cos (cid:16) √ γ r (cid:17) sin (cid:16) √ γ r (cid:17) cos (cid:16) √ γ r (cid:17) . (39)We can thus adapt expressions (33) for W in regions 1and 3 W (1) ( r + , r − ) = − W (3) ( r + , r − ) = − (cid:90) r + − r − x ( r ) R (2 r, r + − r − ) dr, (40)whereas expression (34) in region 2 now reads W (2) ( r + , r − ) = R ( r m − r − , r + − r − ) W (3) ( r m , r − )+ R ( r + + r m , r + − r − ) W (1) ( r + , − r m ) . (41)In expressions (40) and (41) we have made, with respectto Eqs. (33) and (34), the replacements ξ → r + , η → r − so that these expressions can be used in Eqs. (19) andthen (17). The knowledge of the expression of W is theregions where both Riemann invariants vary (regions 1,2and 3) makes it possible to compute r ± ( x, t ) in theseregions, as detailed in Ref. [14]. The density and velocityprofiles are then obtained from Eqs. (12) and (14). In the simple-wave regions, where one of the Riemanninvariants is constant, the hodograph transform does notapply. In such a region however, the solution of the prob-lem is relatively easy. For instance, during the initialstages of evolution, when the regions I and III still ex-ist, we can solve the problem by means of the method ofcharacteristics. We look for a solution in the form x − v + ( r + , − r ) t = h III ( r + ) , for region III,x − v − ( r , r − ) t = h I ( r − ) , for region I. (42)where the functions h III and h I are determined byboundary conditions: the simple wave solution shouldmatch with the solution in regions where both Riemanninvariants vary at the interface between the two regions( III and 3 or I and 1). This yields from Eqs. (17) x − v + ( r + , − r ) t = ∂W (3) ( r + , − r ) ∂r + , (43)for the simple-wave region III , and x − v − ( r , r − ) t = ∂W (1) ( r , r − ) ∂r − , (44)for region I .After a certain time, two new simple-wave regions ap-pear, which are denoted as II l and II r in Fig. 1). Simi-larly, we get x − v + ( r + , − r ) t = ∂W (2) ( r + , − r ) ∂r + , (45)for region II r , and for region II l : x − v − ( r , r − ) t = ∂W (2) ( r , r − ) ∂r − . (46)From the results presented in Secs. 3.2 and 3.3 we ob-tain a complete description of the dispersionless stageof evolution of the system. The top row of Fig. 1 dis-plays a comparison of the results of this approach withthe numerical solution of Eq. (4). There is a very goodagreement up to the wave breaking moment. For sub-sequent time ( t > t WB ) the dispersionless evolution be-comes multi-valued in some regions of space (indicatinga breakdown of the approach and the occurrence of aDSW), but remains accurate in others which we denoteas “dispersionless regions” below. Note that these re-gions can be precisely defined only after the extension ofthe DSW has been properly determined. This questionis addressed in the following section.
4. DISPERSIVE SHOCK WAVES
The basic mathematical approach for the descriptionof DSWs is based on the Whitham modulation theory[15, 16]. This treatment relies on the large differencebetween the fast oscillations within the wave and theslow evolution of its envelope. It results in the so-calledWhitham modulation equations which constitute a verycomplex system of nonlinear first order differential equa-tions. Whitham’s great achievement was to be able totransform such a system — in the very important anduniversal case of the Korteweg-de Vries (KdV) equa-tion — into a diagonal (Riemann) form analogous to thesystem (10). This achievement enabled Gurevich andPitaevskii [21] to successfully apply the Whitham the-ory to the description of DSWs dynamics for the KdVequation in a system experiencing simple wave breaking.It became clear later that the possibility to diagonal-ize the system of Whitham equations is closely relatedwith the complete integrability of the KdV equation, andWhitham theory was extended to many completely inte-grable equations, as described, for example, in the re-views [53] and [54].However, a large number of equations are not com-pletely integrable — among which Eq. (4) — and requirethe development of a more general theory for describingDSWs. An important success in this direction was ob-tained by El who developed in Ref. [30] a method thatmade it possible to find the main parameters of a DSWarising from the evolution of an initial discontinuity. Re-cently, it was shown in Ref. [39] that El’s method canbe generalized to a substantial class of initial conditions,such as simple waves, for which the limiting expressionsfor the characteristic velocities of the Whitham systemat the edges of the DSW are known from general con-siderations as being equal to either the group velocity ofthe wave at the boundary with a smooth solution, or thesoliton velocity at the same boundary. This information,together with the knowledge of the smooth solution in thedispersionless regions, is enough to find the law of motionof the corresponding edge of the DSW for an arbitraryinitial profile. We will use the methods of Ref. [39] tofind the dynamics of the edges of the DSW in our case.The specific dispersive properties of the system underconsideration enter into the general theory in the form ofWhitham’s “number of waves” conservation law [15, 16] ∂k∂t + ∂ω ( k ) ∂x = 0 , (47)where k = 2 π/L and ω = kV are the wave vector andthe angular frequency of a single phase nonlinear wave, L being its wavelength and V its phase velocity. Because ofthe large difference between the scales characterizing theenvelope and those characterizing the oscillations withinthe wave, Eq. (47) is valid along the whole DSW. Atthe small amplitude edge of the DSW the function ω ( k )becomes the linear dispersion relation ω lin ( k ) in the sys-tem. In our case, this dispersion law is obtained by lin-earizing equations (6) around the uniform state ρ = ρ , u = u (we keep here a nonzero value of u for futureconvenience); that is, we write ρ ( x, t ) = ρ + ρ (cid:48) ( x, t ) and u ( x, t ) = u + u (cid:48) ( x, t ), where | ρ (cid:48) | (cid:28) ρ and | u (cid:48) | (cid:28) u . A plane wave expansion of ρ (cid:48) and u (cid:48) immediately yields ω lin ( k ) ≡ ku ± k (cid:115) ρ (1 + γρ ) + k . (48)Thus, if one can calculate the wave number at the small-amplitude edge of the DSW, one obtains the speed ofpropagation of this edge as being equal to the correspond-ing group velocity v g ( k ) = dω lin /dk .Unfortunately, this approach cannot be appliedstraightforwardly at the soliton edge of a DSW. In spiteof that, El showed [30] that under some additional as-sumptions one can obtain from Eq. (47) an ordinary dif-ferential equation relating the physical variables alongthe characteristic of Whitham equations at the solitonedge of the DSW. The idea is based on the following re-mark: At the tails of a soliton, the density profile hasan exponential form ∝ exp( − (cid:101) k | x − V s t | ), where V s is thespeed of the soliton. Of course a soliton propagates atthe same velocity as its tails, which, being a small pertur-bations of the background, obey the linear dispersion law ω = ω lin ( k ). Hence, as was noticed by Stokes [55] (seealso Sec. 252 in [56]), we arrive at the statement that thesoliton velocity V s is related to its “inverse half-width” (cid:101) k by the formula V s = (cid:101) ω ( (cid:101) k ) / (cid:101) k, where (cid:101) ω ( (cid:101) k ) ≡ − iω lin ( i (cid:101) k ) . (49)However, the “solitonic counterpart” ∂ (cid:101) k∂t + ∂ (cid:101) ω ( (cid:101) k ) ∂x = 0 (50)of Eq. (47) does not apply in all possible DSW config-urations, even at the soliton edge. Interestingly, it doesapply for a step-like type of initial conditions, and Elused this property to determine the inverse half-widthand velocity of a soliton edge, using a procedure similarto the one employed at the small amplitude edge. As aresult, a number of interesting problems were successfullyconsidered by this method, see, e.g., Refs. [31–38].To go beyond the initial step-like type of problems, onehas to use some additional information about the prop-erties of the Whitham modulation equations at the edgesof DSWs (see Refs. [39, 40]). For integrable equations, itis known that the “soliton number of waves conservationlaw” (50) is valid in the case where the DSW is triggeredby a simple wave breaking [39]. Therefore, it is natural toassume that Eq. (50) also applies for non-integrable equa-tions in situations where the pulse considered propagatesinto a uniform and stationary medium and experiencesa simple wave breaking.We can consider that the smooth solution of the dis-persionless equations is known from the approach pre-sented in Sec. 3. At the boundary between the disper-sionless simple wave region and the DSW, Eqs. (47) and More precisely: propagates into a medium for which r + and r − are constant. (50) reduce to ordinary differential equations which canbe extrapolated to the whole DSW. Their solution, withknown boundary condition at both edges, yields the wavenumber k and the inverse half-width of solitons (cid:101) k at theboundary with the smooth part of the pulse. Conse-quently, the corresponding group velocity or the solitonvelocity at a DSW edge can be expressed in terms ofthe parameters of the smooth solution at its boundarywith a DSW. These velocities can be considered as thecharacteristic velocities of the limiting Whitham equa-tions at this edge, which makes it possible to representthese equations in the form of first order partial differ-ential equations after a hodograph transformation. Thecompatibility condition for this partial differential equa-tion with a smooth dispersionless solution gives the lawof motion of this DSW edge. If the soliton solution of thenonlinear equation at hand is known, then the knowledgeof its velocity makes it possible to also determine its am-plitude at the boundary of the DSW. In this paper, weshall apply this method to study the evolution of initialsimple-wave pulses in the generalized NLS equation (4).We suppose that in the whole region of the DSW, onehas r − = u − √ γ arctan √ γρ = − √ γ arctan √ γρ = − r = cst . (51)This occurs for instance in the DSW formed at the rightof the right propagating bump issued from the initial con-dition (26) — for which the behavior of the dispersionlessRiemann invariants is sketched in Fig. 1 — and also inthe model cases studied in Secs. 4.1.1 and 4.2 below.Equation (10b) is satisfied identically by virtue of ourassumption (51). Also, since r − is constant, r + and u can be considered as functions of ρ only [see Eqs. (12)and (13)]: u ( ρ ) = 2 √ γ (arctan √ γρ − arctan √ γρ ) ,r + ( ρ ) = 1 √ γ (2 arctan √ γρ − arctan √ γρ ) , (52)and Eq. (10a) for r + thus takes the form ∂ρ∂t + ( u ( ρ ) + c ( ρ )) ∂ρ∂x = 0 , (53)where the expression of the local sound velocity is givenin (8). The general solution (17a) of Eq. (10a) here spe-cializes to x − ( u ( ρ ) + c ( ρ )) t = w + ( r + ( ρ ) , − r ) . (54)All the relevant characteristic quantities of the DSWformed after the wave breaking of such a pulse can beconsidered to be functions of ρ only. We can thus write the right propagating soliton dispersion law (49) in theform (cid:101) ω ( (cid:101) k ) = (cid:101) k (cid:20) u ( ρ ) + (cid:113) c ( ρ ) − (cid:101) k (cid:21) = (cid:101) k [ u ( ρ ) + c ( ρ ) (cid:101) α ( ρ )] , (55)where (cid:101) α ( ρ ) = (cid:115) − (cid:101) k c ( ρ ) . (56)We then have (cid:101) k ( ρ ) = 2 c ( ρ ) (cid:112) − (cid:101) α ( ρ ) , (57)where (cid:101) α ( ρ ) is yet an unknown function.As implied by El’s method, along the soliton edge ofthe DSW one considers that r − is also a constant, andthat, at all time, the quantities u , c , r + , (cid:101) k , (cid:101) α and (cid:101) ω canbe considered as functions of ρ only. Then the solitonic“number of wave conservation” (50) yields, with accountof (53): d (cid:101) ωdρ = [ u ( ρ ) + c ( ρ )] d (cid:101) kdρ . (58)Substitution of (55) and (57) into this equation yields anordinary differential equation for (cid:101) α ( ρ ), d (cid:101) αdρ = − (1 + (cid:101) α )(1 + 3 γρ + 2 (cid:101) α (1 − γρ ))2 ρ (1 + γρ )(1 + 2 (cid:101) α ) , (59)which should be solved with the boundary condition ap-propriate to the situation considered (see below).The same method is employed for studying the motionof the small amplitude edge of the DSW. This edge prop-agates at the linear group velocity, and for determiningthe relevant wave-vector, we rewrite equation (48) in theform ω lin ( k ) = k [ u ( ρ ) + c ( ρ ) α ( ρ )] , (60)that is α ( ρ ) = (cid:115) k c ( ρ ) , (61)and k ( ρ ) = 2 c ( ρ ) (cid:112) α ( ρ ) − . (62)A simple calculation leads to the following differentialequation dαdρ = − (1 + α )(1 + 3 γρ + 2 α (1 − γρ ))2 ρ (1 + γρ )(1 + 2 α ) , (63)which actually coincides with the equation (59). Nowwe “extrapolate” this equation across the entire DSWan determine the wave-vector as the value of k ( ρ ) at thedensity of the small amplitude edge. The appropriateboundary condition for solving (63) is fixed at the soli-tonic edge, and depends on the configuration under study(see below).0 xρ ( x, ρ (a) ρ x ( ρ ) ρ (b) Fig. 3 . (a) Initial profile ρ ( x ) of a monotonous positive pulse.(b) Inverse function x ( ρ ). In all this section we consider an initial condition witha region of increased light intensity over a stationary uni-form background. r − As a first illustration of the method we consider aninitial condition for which the Riemann invariant r − hasthe constant value − r for all x . We shall start with apulse with the following initial intensity distribution: ρ ( x ) = (cid:40) ρ + (cid:101) ρ ( x ) if x < ,ρ if x ≥ , (64)where (cid:101) ρ ( x ) is a decreasing function of x with a verticaltangent at x → − , see Fig. 3(a). From (9) one thenobtains u ( x,
0) = u ( x ) with u ( x ) = 2 √ γ arctan (cid:112) γρ ( x ) − r . (65)This non-standard type of initial condition is of interestbecause it enables to test the theory just presented in aparticularly simple setting in which (i) the evolution ofthe non-dispersive part of the system is that of a simplewave, i.e., one can replace w + in the right hand side ofEq. (54) by x ( ρ ), where x ( ρ ) is the function inverse to theinitial distribution of the intensity ρ ( x ) (see Fig. 3(b)),(ii) the wave breaking occurs instantaneously, at the frontedge of the pulse. This wave breaking is followed by theformation of a DSW with a small-amplitude wave at itsfront edge and with solitons at its trailing edge, at theboundary with the smooth part of the pulse [describedby the equation (54)]. Therefore, we can find the lawof motion of the rear solitons edge of the DSW by themethod of Ref. [39].Eq. (59) should be solved with the boundary condition (cid:101) α ( ρ ) = 1, that is, the soliton inverse half-width vanishestogether with its amplitude at the small-amplitude edge. When the function (cid:101) α ( ρ ) is known, the velocity of thefront edge can be represented as the soliton velocity V s = (cid:101) ω (cid:101) k = u ( ρ s ) + c ( ρ s ) (cid:101) α ( ρ s ) , (66)where ρ s is the density at the solitonic edge of the DSW.We notice again that V s is the characteristic velocityof the Whitham modulation equations at the solitonedge. The system being nonintegrable, the correspond-ing Whitham equation is not known, but its limiting format the solitonic edge can be written explicitly: along thisboundary one has dx s − V s dt = 0, where x s ( t ) is the posi-tion of the solitonic edge at time t and ρ s = ρ ( x s ( t ) , t ). Itis convenient to re-parametrize all the quantities in termof ρ s : x ( ρ s ) and t ( ρ s ). This leads to ∂x s ∂ρ s − V s ( ρ s ) ∂t∂ρ s = 0 . (67)This equation should be compatible with Eq. (53) at thesolitonic boundary between the DSW and the dispersion-less region gives, with account of (54). Here, the specificinitial condition we have chosen (with r − ( x,
0) = cst)simplifies the problem because one has w + ( r + ( ρ ) , − r ) =¯ x ( ρ ). The compatibility equation then reads √ ρ s γρ s (1 − (cid:101) α ) dtdρ s + 3 + γρ s √ ρ s (1 + γρ s ) t = − dx ( ρ s ) dρ s . (68)One sees from Fig. (3) that for the initial condition weconsider, the wave breaking occurs instantly at t = 0,with ρ s = ρ . Hence Eq. (68) should be integrated withthe initial condition t ( ρ ) = 0. The corresponding solu-tion reads t ( ρ s ) = (cid:90) ρ s ρ (1 + γρ ) x (cid:48) ( ρ ) G ( ρ s , ρ ) √ ρ ( (cid:101) α ( ρ ) − dρ, (69)where x (cid:48) = dx/dρ and G ( ρ s , ρ ) = exp (cid:18)(cid:90) ρ s ρ γρ (cid:48) ρ (cid:48) (1 + γρ (cid:48) )( (cid:101) α ( ρ (cid:48) ) − dρ (cid:48) (cid:19) . (70)Consequently, x s ( ρ s ) = [ u ( ρ s ) + c ( ρ s )] t ( ρ s ) + x ( ρ s ) , (71)where u ( ρ ) and c ( ρ ) are given by Eqs. (52) and (8), re-spectively.As a generalisation of the model initial profile repre-sented in Fig. 3, we now consider the case where theinitial distribution ρ ( x ) is non monotonous, with a singlemaximum. This amounts to assume that (cid:101) ρ in (64) hasa maximum ρ m at x = x m <
0, tends to 0 fast enoughas x → −∞ , with still a vertical tangent at x → − ,see Fig. 4(a). This last condition means that we sup-pose again for convenience that the wave breaks at themoment t = 0. The reciprocal function of ρ ( x ) becomesnow two-valued and we denote its two branches as x ( ρ )1 xρ ( x, ρ ρ m x m (a) ρx ( ρ ) ρ ρ m x m x x (b) Fig. 4 . (a) Initial profile ρ ( x ) of a non-monotonous posi-tive pulse. (b) The inverse function x ( ρ ) is represented bytwo branches x ( ρ ) and x ( ρ ). Here x m is the point where ρ reaches its maximal value ρ m . and x ( ρ ), see Fig. 4(b). In an initial stage of develop-ment of the DSW, the formulas obtained previously fora monotonous initial condition straightforwardly apply.This occurs when the maximum of the smooth part ofthe profile issued from the initial distribution has not yetpenetrated the DSW. In this case the soliton edge propa-gates along the branch x ( ρ ) and one should just use thisfunction instead of x ( ρ ) in Eqs. (69) and (71). A secondstage begins when the maximum of the smooth part ofthe profile reaches the DSW. In this case one can show that Eqs. (69) and (71) are modified to t ( ρ s ) = (cid:90) ρ m ρ (1 + γρ ) x (cid:48) ( ρ ) G ( ρ s , ρ ) √ ρ ( (cid:101) α ( ρ ) − dρ + (cid:90) ρ s ρ m (1 + γρ ) x (cid:48) ( ρ ) G ( ρ s , ρ ) √ ρ ( (cid:101) α ( ρ ) − dρ, (72a) x s ( ρ s ) = [ u ( ρ s ) + c ( ρ s )] t ( ρ s ) + x ( ρ s ) . (72b)To determine the position of the small-amplitude edgewe solve Eq. (63) with the boundary condition α ( ρ m ) = 1 , (73)which corresponds to the moment at which the solitonedge reaches the maximal point of the initial distribution.The corresponding solution α ( ρ ) yields the spectrum (62)of all possible wave numbers k ( ρ ), with ρ ranging from ρ to ρ m . The maximal value of the group velocity dω lin /dk is reached at ρ = ρ and it provides the asymptotic valueof the velocity of the small amplitude edge: dx r dt = dω lin dk (cid:12)(cid:12)(cid:12)(cid:12) ρ = ρ = √ ρ γρ α ( ρ ) − α ( ρ ) . (74)As an illustration of the accuracy of the method, wenow compare the theoretical results with those of a nu-merical solution of the generalized NLS equation (4) for See the detailed study of a similar situation in Sec. 4.1.2.
50 0 50 100 150 x ρ ( x , t ) ρ m ρ t = 0 t = 100 Fig. 5 . Blue solid curve: simple-wave initial condition (75)with ρ = 0 . ρ m = 1 and x = 50. Red solid curve: resultingdispersive shock wave computed by numerically solving Eq.(4) with γ = 1. the non-monotonous initial intensity distribution ρ ( x ) = ρ + ( ρ m − ρ ) (cid:16) xx + 1 (cid:17) if x ∈ [ − x , ,ρ elsewhere , (75)see the blue curve in Fig. 5. The typical form of theDSW generated by such a pulse is represented in the samefigure by a red curve. For the particular initial profile(75) the branch x ( ρ ) shrinks to zero, the first term inthe right hand side of (72a) cancels and we need takeinto account only the contribution of the second branch,given by the expression x ( ρ ≥ ρ ) = x (cid:18) ρ − ρ ρ m − ρ (cid:19) / − x . (76)For the numerical simulation, we took instead of the ide-alized profile (75) the numerical initial condition: ρ ( x ) = ρ + ( ρ m − ρ ) (1 + x/x ) − tanh( x/w )2 , (77)with w = 1, x = 50, ρ = 0 . ρ m = 1. Equation(77) yields a maximum value of ρ which is not exactlyequal to ρ m = 1, but to ρ max = 0 . x s ( t ) and density ρ s ( t ), which are represented in Fig. 6by continuous colored curves. For the analytic determi-nation of ρ s ( t ) and x s ( t ) we used in Eqs. (72) two possiblevalues of ρ m : 1 (green curve in Fig. 6) and 0.9362 (redcurve). As one can see, there is no much difference, andthe theoretical results both agree well with the numer-ical simulations (black dots). Of course, the results atshort time for ρ s are better when one takes the initial ρ m = 0 . x s t ρ s Fig. 6 . Position x s ( t ) and density ρ s ( t ) of the solitonic edgeof the DSW induced by the initial profile (75) and (65). Theblack dots are obtained by numerical integration of Eq. (4).The colored solid lines are obtained from the analytic formulas(72). In each plot the green theoretical solid line is obtainedby taking ρ m = 1 and the red one by taking ρ m = 0 . and the results for x s almost do not depend on the choice ρ m = 1 or 0.9362.For our initial parameters, the analytic formula (74)gives the asymptotic propagation velocity of the rightsmall-amplitude edge of the DSW: dx r /dt = 1 .
11. Thenumerical value of the velocity is 1 .
12 for our choice of ini-tial condition, an agreement which should be consideredas very good since the position of the small-amplitudeedge is not easily determined, as is clear from Fig. 5, andwe evaluate it by means of an approximate extrapolationof the envelopes of the wave at this edge. (26)
We now consider the initial density profile (26) with u ( x,
0) = 0; that is, an initial state which is not of simplewave type, and for which the formulas of the previousSec. 4.1.1 are not directly applicable. In this case, thewave breaking of the right propagating part of the pulseoccurs in region
III , where r − is constant and equalto − r at ρ s ( t WB ) = ρ ( x s ( t WB ) , t WB ) = ρ . Hence Eq.(59) should be integrated with the boundary condition (cid:101) α ( ρ ) = 1. In a first stage of its development , the DSWmatches the dispersionless profile described by Eq. (43)where W (3) ( r + , r − ) is given by Eq. (40) and x ( r ) by (30).In the region of interest for us r − = − r , so all quan-tities in Eq. (43) depend on r + only. Differentiation It is explained below what happens in a second stage, seeEq. (91). with respect to r + yields, at the solitonic edge where dx s = V s dt :( V s − v + ) dtdr + − dv + dr + t = ∂ W (3) ∂r , (78)where all quantities are evaluated at r s = r + ( x s ( t ) , t ) and r − = − r . It follows from (13) that dr s = c ( ρ s ) dρ s /ρ s ,yielding dv + dr + = 1 + ρ s c ( ρ s ) dcdρ s = 3 + γρ s γρ s ) . (79)For explicitly evaluating the right hand side of Eq. (78) itis convenient to write R ( r , r ) = f ( r ) /f ( r ) in Eq. (40)[see expressions (38) and (39)] which leads to ∂W (3) ∂r + = − f (cid:48) ( r + − r − ) f ( r + − r − ) W (3) + x ( r + ) f (2 r + ) f ( r + − r − ) , (80)and ∂ W (3) ∂r = (cid:20) f (cid:48) ( r + − r − ) f ( r + − r − ) − f (cid:48)(cid:48) ( r + − r − ) f ( r + − r − ) (cid:21) W (3) − x ( r + ) f (cid:48) ( r + − r − ) f (2 r + ) − f ( r + − r − ) f (cid:48) (2 r + ) f ( r + − r − )+ dxdr + f (2 r + ) f ( r + − r − ) . (81)At the wave-breaking time t WB one has r s = r ; Eqs. (40),(80) and (81) then directly yield W (3) = 0, ∂W (3) /∂r + = x ( r ) = l and [using Eqs. (30) and (29)] ∂ W (3) ∂r = dxdr (cid:12)(cid:12)(cid:12)(cid:12) r = r = − l (1 + γρ ) √ ρ ρ m − ρ . (82)At the wave breaking V s = v + [ (cid:101) α ( ρ ) = 1 in (66)] andEq. (78) thus yields t WB = lρ m − ρ √ ρ (1 + γρ ) γρ . (83)This result is in agreement with the findings of Ref. [14]and reduces to the one obtained in Ref. [12] for the non-linear Schr¨odinger equation in the limit γ = 0. Forthe example considered in Fig. 1, expression (83) gives t WB = 12 . t WB one should solve Eq. (78) us-ing the generic form (81) of its right hand side. It is con-venient to multiply this equation by ρ s /c ( ρ s )( dr s /dρ s ) =1, which yields ρ s ( (cid:101) α ( ρ s ) − dtdρ s − γρ s γρ s ) t = ∂ W (3) ( r s , − r ) ∂r s . (84)The solution reads t ( ρ s ) = (cid:90) ρ s ρ ∂ W (3) ∂r G ( ρ s , ρ ) ρ ( (cid:101) α ( ρ ) − dρ, (85)3where the expression of G is given in (70) and r + and x should be computed as functions of ρ via Eqs. (52), (29)and (30). Once t ( ρ s ) is known, x s ( ρ s ) is determined fromEq. (43). One can check that expression (85) yields thecorrect result (83) for t WB : When ρ is close to ρ one getsfrom Eq. (70) and (59) (with (cid:101) α ( ρ ) = 1): G ( ρ s , ρ ) (cid:39) (cid:18) ρ − ρ ρ s − ρ (cid:19) / . (86)Inserting this expression back into (85) one obtains t ( ρ ) = − γρ )3 + γρ (cid:18) ∂ W (3) ∂r (cid:19) r + = r . (87)That is, as expected, t ( ρ ) = t WB , where t WB is given byEq. (83).The solution (85) is acceptable only up to a time whichwe denote as t | at which r s = r m , i.e., ρ s reaches amaximum which we denote as ρ M where [see Eq. (14)] ρ M = γ − tan [ √ γ ( r m + r ) / . (88)At times larger than t | the region III disappears andthe DSW is in contact with region II r of the dispersion-less profile in which Eq. (45) applies. Instead of (84)one should thus solve ρ s ( (cid:101) α ( ρ s ) − dtdρ s − γρ s γρ s ) t = ∂ W (2) ( r s , − r ) ∂r s , (89)with the initial condition t ( ρ M ) = t | . W (2) in the aboveequation is given by Eq. (41) which can be cast in theform W (2) ( r + , r − ) = 1 f ( r + − r − ) (cid:104) (cid:90) r m − r − x ( r ) f (2 r ) dr + (cid:90) r m r + x ( r ) f (2 r ) dr (cid:105) . (90)This yields ∂W (2) ∂r + = − f (cid:48) ( r + − r − ) f ( r + − r − ) W (2) − x ( r + ) f (2 r + ) f ( r + − r − ) , (91)and ∂ W (2) ∂r = (cid:20) f (cid:48) ( r + − r − ) f ( r + − r − ) − f (cid:48)(cid:48) ( r + − r − ) f ( r + − r − ) (cid:21) W (2) + 2 x ( r + ) f (cid:48) ( r + − r − ) f (2 r + ) − f ( r + − r − ) f (cid:48) (2 r + ) f ( r + − r − ) − dxdr + f (2 r + ) f ( r + − r − ) . (92) The equivalent time was denoted as t A | B in Ref. [12]. x s t ρ s ρ M Fig. 7 . Position x s ( t ) and density ρ s ( t ) of the solitonic edgeof the DSW induced by the initial profile (26) with ρ =0 . ρ m = 2 and l = 20. The continuous red curves are thetheoretical results. The big red dots locate the wave breakingtime t WB = 12 .
12, position x s ( t WB ) = 25 .
71 and intensity ρ ( x s ( t WB ) , t WB ) = ρ . The dotted curves are extracted fromnumerical simulations (see the text). The dashed gray line inthe upper panel corresponds to a soliton edge which wouldpropagate at the background sound velocity. The horizontalred line in the lower panel marks the position of the maximum ρ M as given by Eq. (88). The solution of (89) reads t ( ρ s ) = t | + (cid:90) ρ s ρ M ∂ W (2) ∂r G ( ρ s , ρ ) ρ ( (cid:101) α ( ρ ) − dρ. (93)Once t ( ρ s ) is determined from (93) x s ( ρ s ) is obtainedfrom (45). Formulas (85), (43), (93) and (45) give para-metric expressions of ρ s ( t ) and x s ( t ) for all t ≥ t WB .We note here that in the limit γ →
0, the solution ofEq. (59) reads (cid:101) α = − (cid:112) ρ /ρ which, from (66), yields V s = √ ρ s = ( r s − r − ), in agreement with the findings ofRef. [12]. Also in this limit, one gets the exact expression G ( ρ s , ρ ) = ( √ ρ − √ ρ ) / ( √ ρ s − √ ρ ) − / and formulas(85) and (93) reduce to the equivalent ones derived forthe NLS equation in Ref. [12].In the generic case γ (cid:54) = 0, the results for x s ( t ) and ρ s ( t )corresponding to the parameters of Fig. 1 are comparedin Fig. 7 with the values extracted from the numericalsolution of Eq. (4) (with γ = 1). The position x s ofthe solitonic edge is easily determined from the numer-ical simulation: it is located between the maximum ofthe smooth part of the intensity and the following zeroof ρ − ρ (respectively black and blues dots in the up-per panel of the figure). The intensity ρ s of the solitonicedge is more difficult to extract from the numerical sim-ulations in the present situation that in the case studiedin Fig. 6. At times close to the theoretical wave-breakingtime, one cannot exactly decide from the numerical in-tensity profile if the oscillations at the boundary of the4right propagating edge of the region of increased inten-sity are linear disturbances (due to dispersive effects) orcorrespond to the birth of a DSW. It is only after a cer-tain amount of time that the oscillations become clearlynonlinear. Even then, the precise location of the inten-sity of the solitonic edge is not easily determined: it iscomprised between the maximum intensity of the smoothpart of the spectrum and the following maximum inten-sity (respectively black and green dots in the lower panelof Fig. 7). The green dots initially yield a clear underestimate of ρ s , and, on the contrary, when the dispersiveshock if fully developed, they indicate a too large value.This is due to the fact that, at large time, the ampli-tude of the oscillations in the vicinity of the soliton edgeincrease within the DSW.We note also that expression (88) gives a simple andaccurate prediction for the maximum of ρ s ( ρ M = 1 inthe case considered in Fig. 7). At times large compareto the time t | at which this maximum is reached, onecan approximately evaluate the behavior of ρ s and x s asfollows: The integrands in Eqs. (85) and (93) are of theform ( ρ s − ρ ) − / F ( i ) ( ρ s , ρ ) where F is a non singularfunction [with i = 2 for Eq. (85) and i = 3 for Eq. (93)].Since at large time ρ s tends to ρ , formula (93) can beapproximated by t ( ρ s ) (cid:39) A ( ρ s − ρ ) / , (94)where A = (cid:90) ρ M ρ (cid:16) F (3) ( ρ , ρ ) − F (2) ( ρ , ρ ) (cid:17) dρ. (95)From there it follows that, at large t , ρ s tends to ρ as t − / and x s exceeds c ( ρ ) t by a factor of order t / .The fact that the soliton edge propagates at a velocityhigher that the speed of sound is illustrated in Fig. 7 bythe difference between x s ( t ) and the gray dashed line ofequation x = x s ( t WB ) + c ( ρ )( t − t WB ).A final remark is in order here. When ρ m gets largecompared with ρ , (cid:101) α ( ρ s )—solution of (59)—may becomenegative. For instance, when γ = 0 this occurs for ρ m ≥ ρ . As noticed in Ref. [32], this is linked to theoccurrence within the DSW of a vacuum point [57] atwhich the density cancels. In this case, one should useanother branch in the dispersion relations (48) and (55),but once this modification is performed, the approachremains perfectly valid and accurate.When γ >
0, a more serious problem occurs at largervalues of ρ m , when (cid:101) α ( ρ s ) reaches − . A rough estimatebased on the regime γ (cid:28) ρ m ≈ ρ . In this case the theoreti-cal approach leads to singularities [see, e.g., Eq. (59)]whereas the numerical simulations do not indicate a dras-tic change of behavior in the dynamics of the DSW. Thispoints to a probable failure of the Gurevich-Meshcherkin-El approach, however the study of this problem is beyondthe scope of the present work and we confine ourselves xρ ( x, ρ (a) ρ x ( ρ ) ρ (b) Fig. 8 . (a) Initial profile ρ ( x ) of a monotonous negative pulse.(b) Inverse function x ( ρ ). to a regime of initial parameters for which (cid:101) α ( ρ s ) remainslarger than − and the theory is applicable. In this section we consider an initial condition with aregion where the density is depleted with respect to thatof the background. For simplicity we only consider thecase where, as in Sec. 4.1.1, the initial Riemann invariant r − ( x,
0) is constant and equal to − r for all x .Let us first consider a monotonous initial pulse ρ ( x,
0) = (cid:40) ρ if x < ,ρ − (cid:101) ρ ( x ) if x ≥ . (96)This initial condition is sketched in Fig. 8. The small-amplitude edge propagates with the group velocity V r ( ρ r ) = dωdk = u ( ρ r ) + √ ρ r γρ r α ( ρ r ) − α ( ρ r ) , (97)where ρ r ( t ) is the density of the small amplitude edge.Following a reasoning analogous to the one already em-ployed in Sec. 4.1.1 for the solitonic edge, we parame-terize all relevant quantities at the small amplitude edgein term of ρ r . This leads to the following differentialequation √ ρ r γρ r α − α − α dtdρ r − γρ r √ ρ r (1 + γρ r ) t = dx ( ρ r ) dρ r , (98)which should be solved with the initial condition ρ r = ρ at t = 0. The solution reads t ( ρ r ) = (cid:90) ρ r ρ α ( ρ )(1 + γρ ) x (cid:48) ( ρ ) G ( ρ r , ρ ) √ ρ ( α ( ρ ) − α ( ρ )) dρ, (99)where x (cid:48) ( ρ ) = dx/dρ and G ( ρ r , ρ ) =exp (cid:18)(cid:90) ρ r ρ α ( ρ (cid:48) )(3 + γρ (cid:48) ) dρ (cid:48) ρ (cid:48) (1 + γρ (cid:48) )( α ( ρ (cid:48) ) − α ( ρ (cid:48) )) (cid:19) . (100)5 xρ ( x ) ρ ρ m x m (a) ρ x ( ρ ) ρ x m ρ m x x (b) Fig. 9 . (a) Initial profile ρ ( x ) of a non-monotonous negativepulse. (b) The inverse function x ( ρ ) is represented by twobranches x ( ρ ) and x ( ρ ). Here x m is the point where ρ takesits minimal value ρ m . Consequently, the position x r of the small amplitude edgeis given by x r ( ρ r ) = [ u ( ρ r ) + c ( ρ r )] t ( ρ r ) + x ( ρ r ) , (101)where u ( ρ ) is given by Eq. (52). Formulas (99) and(101) determine, in a parametric form, the coordinates x r ( t ) and ρ r ( t ) of the small-amplitude edge.In the case of a non-monotonous negative pulse, suchas the one represented in Fig. 9, the approach has tobe modified in a manner similar to that exposed in Sec.4.1.1 for the non-monotonous profile displayed in Fig. 4.We do not write down the corresponding formulas for notoverloading the paper with almost identical expressions.We only indicate that ρ m denotes now the minimal valueof the density in the initial distribution.For numerically testing our approach we consider anon-monotonous initial profile of the form ρ ( x ) = ρ − ( ρ m − ρ ) (cid:16) xx − (cid:17) if x ∈ [0 , x ] ,ρ elsewhere , (102)see the blue curve in Fig. 10. In this case one has x ( ρ ) =0 and x ( ρ ≤ ρ ) = x − x (cid:18) ρ − ρρ − ρ m (cid:19) / . (103)The initial velocity distribution u ( x,
0) is determinedfrom (102) by imposing that r − ( x,
0) = − r . The typicalwave pattern for the light evolved from such a pulse isshown in Fig. 10.The motion of the right, small amplitude edge is de-termined as previously explained. Eq. (63) should besolved with the initial condition α ( ρ m ) = 1. Then, in thenon-monotonous case we consider here, Eq. (98) shouldbe solved with an initial condition different from the oneused in the case of the initial profile (96). The appropri-ate initial value is here ρ r = ρ m at t = 0. As a result, thelower boundary of the integration region in (99) shouldbe changed from ρ to ρ m . The theoretical results com-pare very well with the one extracted from the numerical −
100 0 100 200 300 4000 . . . ρ ρ m x ρ Fig. 10 . Density profile ρ ( x, t = 300) (red solid curve), whichresults from the time evolution of the simple-wave initial con-dition (102) with ρ = 1, ρ m = 0 . x = 300 (blue solidcurve).
200 300 400 500 600 700 800200400 x r
200 300 400 500 600 700 8000 . . . . tρ r Fig. 11 . Coordinate x r ( t ) and density ρ r ( t ) of the small-amplitude edge of the DSW induced by the initial profile(102). The black dots are obtained by numerical integra-tion of Eq. (4). The green solid lines are obtained from theanalytic solution, i.e., from Eqs. (99) and (101). simulations, as illustrated in Fig. 11. Note that the de-termination of the small amplitude edge of the numericalDSW is a delicate task, which is simplified if the DSWcontains a large number of oscillations. This is the reasonwhy we had to chose here an initial profile of extensionlarger than the one considered in the previous section4.1.1 ( x = 300 instead of 50). Note however that thecase considered here presents an advantage compared tothe similar one considered in Sec. 4.1.1: the motion of6 . . . . . . t V s Analytical SolutionNumerical Calculations
Fig. 12 . Black dots: numerically determined velocity V s ( t )of the soliton edge. The initial conditions are specified inEq. (102) with ρ = 1, ρ m = 0 . x = 300 (blue).Continuous red line: fit of the numerical data by the for-mula: V s ( t ) = V num s + (cid:80) i =1 a i t − i . One obtains V num s =0 .
39, in good agreement with the theoretical prediction fromEq. (104): V s ( t → ∞ ) = 0 .
37 (blue dashed line). the small-amplitude edge is less dependent of the precisecharacteristics of ρ ( x ) at its extremum than the motionof the solitonic edge. As a result, we need not discusshere the precise numerical implementation of the ideal-ized discontinuous profile (102), contrarily with what hasbeen done in Sec. 4.1.1.At asymptotically large time the velocity of the trailingsoliton generated from a localized pulse is determinedfrom the knowledge of (cid:101) α ( ρ ), where (cid:101) α ( ρ ) is the solutionof equation (59) with the boundary condition (cid:101) α ( ρ m ) = 1.Within this approximation, the asymptotic velocity ofthe soliton edge reads V s ( t → ∞ ) = √ ρ γρ (cid:101) α ( ρ ) . (104)The comparison of the analytic prediction (104) with our numerical simulations is not difficult because the soli-ton edge is easily identified in the numerical density pro-files, and, indeed, a leading soliton is easily located atthis edge of the numerically determined DSW. As wesee in Fig. 12, its velocity tends at large t to a con-stant value. In this figure, the numerical result is fittedwith the empirical formula V s ( t ) = V num s + (cid:80) i =1 a i t − i ,where V num s and a i are fitting parameters. The trendis in excellent agreement with the prediction (104) sinceone obtains V num s = 0 .
39 whereas from (104) one expects V s ( t → ∞ ) = 0 .
5. CONCLUSION
In this work we have presented a theoretical study ofthe dynamics of spreading of a pulse of increased (ordecreased) intensity propagating over a uniform back-ground. The initial non-dispersive stage of evolution isdescribed by means of Riemann’s method, the result ofwhich compares very well with numerical simulations.After the wave breaking time, we have studied the be-havior of the resulting dispersive shock wave at its edges,by mean of the modification of El’s method presented inRef. [39]. Here also the results compare very well withthe ones of numerical simulations.Our work represents a comprehensive theoretical de-scription of the spreading, of the wave breaking and of thesubsequent formation of a dispersive shock in a realisticsetting for a system described by a nonintegrable nonlin-ear equation. In particular, our approach yields a simpleanalytic expression for the wave-breaking time, even insituations where the initial density and velocity profilesdo not reduce to a simple wave configuration. Also, inview of future experimental studies, we have devoted aspecial attention to the determination of the position andintensity of the solitonic edge of the DSW issued from thespreading of a region of increased light intensity.
ACKNOWLEDGMENTS
We thank T. Bienaim´e and M. Bellec for fruitful dis-cussions. J.-E.S thanks Moscow Institute of Physicsand Technology and Institute of Spectroscopy of RussianAcademy of Sciences for kind hospitality during his in-ternship there, when this research was started. S.K.I andA.M.K thank RFBR for financial support of this studyin framework of the project 20-01-00063. [1] Z. Dutton, M. Budde, C. Slowe, and L. V. Hau, Ob-servation of Quantum Shock Waves Created with Ultra-Compressed Slow Light Pulses in a Bose-Einstein Con-densate, Science , 663 (2001).[2] T. P. Simula, P. Engels, I. Coddington, V. Schweikhard,E. A. Cornell, and R. J. Ballagh, Observations on SoundPropagation in Rapidly Rotating Bose-Einstein Conden-sate, Phys. Rev. Lett. , 080404 (2005). [3] M. A. Hoefer, M. J. Ablowitz, I. Coddington, E. A. Cor-nell, P. Engels, and V. Schweikhard, Phys. Rev. A ,023623 (2006).[4] N. Ghofraniha, C. Conti, G. Ruocco, and S. Trillo,Shocks in nonlocal media, Phys. Rev. Lett. , 043903(2007).[5] N. Ghofraniha, S. Gentilini, V. Folli, E. DelRe, andC. Conti, Shock waves in disordered media, Phys. Rev. Lett. , 243902 (2012).[6] G. Xu, A. Mussot, A. Kudlinski, S. Trillo, F. Copie,and M. Conforti, Shock wave generation triggered by aweak background in optical fibers, Optics Letters , 11(2016).[7] G. Xu, M. Conforti, A. Kudlinski, A. Mussot, and S.Trillo, Dispersive dam-break flow of a photon fluid, Phys.Rev. Lett. 118, 254101 (2017).[8] G. Marcucci, M. Braidotti, S. Gentilini, and C. Conti,Time asymmetric quantum mechanics and shock waves:Exploring the irreversibility in nonlinear optics, Ann.Phys. (Berlin) , 1600349 (2017).[9] J. Nu˜no, C. Finot, G. Xu, G. Millot, M. Erkintalo, andJ. Fatome, Vectorial dispersive shock waves in opticalfibers, Commun. Phys. , 138 (2019).[10] M. G. Forest, C.-J. Rosenberg, and O. C. Wright III,On the exact solution for smooth pulses of the defocus-ing nonlinear Schr¨odinger modulation equations prior tobreaking, Nonlinearity , 2287 (2009)[11] S. K. Ivanov and A. M. Kamchatnov, Collision of rar-efaction waves in Bose-Einstein condensates Phys. Rev.A , 013609 (2019).[12] M. Isoard, A. M. Kamchatnov, and N. Pavloff, Wavebreaking and formation of dispersive shock waves in adefocusing nonlinear optical material, Phys. Rev. A ,053819 (2019).[13] M. Isoard, A. M. Kamchatnov, and N. Pavloff, Short-distance propagation of nonlinear optical pulses, in Compte-rendus de la 22 e rencontre du Non Lin´eaire , Eds.E. Falcon, M. Lefranc, F. P´etr´elis, C.-T. Pham, (Non-Lin´eaire Publications, Saint-´Etienne du Rouvray, 2019).[14] M. Isoard, A. M. Kamchatnov, and N. Pavloff, Disper-sionless evolution of inviscid nonlinear pulses, EPL ,64003 (2020).[15] G. B. Whitham, Non-linear dispersive waves, Proc. R.Soc. London, Ser. A , 238 (1965).[16] G. B. Whitham, Linear and Nonlinear Waves, (WileyInterscience, New York, 1974).[17] A. M. Kamchatnov,
Nonlinear Periodic Waves and TheirModulations—An Introductory Course , (World Scientific,Singapore, 2000).[18] G. A. El and M. A. Hoefer, Dispersive shock waves andmodulation theory, Physica D , 11 (2016).[19] M. G. Forest and J. E. Lee, Geometry and modulationtheory for periodic nonlinear Schr¨odinger equation, in
Oscillation Theory, Computation, and Methods of Com-pensated Compactness,
Eds. C. Dafermos et al. , IMAVolumes on Mathematics and its Applications , p. 35(Springer, N.Y., 1986).[20] M. V. Pavlov, Nonlinear Schr¨odinger equation and theBogolyubov-Whitham method of averaging, Teor. Mat.Fiz. , 351 (1987) [Theoret. Math. Phys. , 584(1987)].[21] A. V. Gurevich and L. P. Pitaevskii, Nonstationary struc-ture of a collisionless shock wave, Zh. Eksp. Teor. Fiz. ,590 (1973) [Sov. Phys. JETP , 291 (1974)].[22] A. V. Gurevich and A. L. Krylov, Dissipationless shockwaves in media with positive dispersion, Zh. Eksp. Teor.Fiz. 92, 1684 (1987) [Sov. Phys. JETP , 944 (1987)].[23] G. El, V. Geogjaev, A. Gurevich, and A. Krylov, Decayof an initial discontinuity in the defocusing NLS hydro-dynamics, Physica D , 186 (1995).[24] G. A. El, A. L. Krylov, General solution of the Cauchy problem for the defocusing NLS equation in the Whithamlimit, Phys. Lett. A
77 (1995).[25] R. Z. Sagdeev, Cooperative Phenomena and Shock Wavesin Collisionless Plasmas, in
Reviews of Plasma Physics,
Ed. M. A. Leontovich, Vol. 4, p. 23, (Consultants Bureau,New York, 1966).[26] G. Couton, H. Mailotte, and M. Chauvet, Self-formationof multiple spatial photovoltaic solitons, J. Opt. B: Quan-tum Semiclassical Opt. , S223 (2004).[27] W. Wan, S. Jia, and J. W. Fleischer, Dispersivesuperfluid-like shock waves in nonlinear optics, Nat.Phys. , 46 (2007).[28] A. V. Gurevich and A. P. Meshcherkin, Expanding self-similar discontinuities and shock waves in dispersive hy-drodynamics, Zh. Eksp. Teor. Fiz. Fluid Mechanics , (Perg-amon, Oxford, 1987).[30] G. A. El, Resolution of a shock in hyperbolic systemsmodified by weak dispersion, Chaos , 037103 (2005).[31] G. A. El, R. H. J. Grimshaw, N. F. Smyth, Unsteady un-dular bores in fully nonlinear shallow-water theory, Phys.Fluids , 027104 (2006).[32] G. A. El, A. Gammal, E. G. Khamis, R. A. Kraenkel, andA. M. Kamchatnov, Theory of optical dispersive shockwaves in photorefractive media, Phys. Rev. A , 053813(2007).[33] G. A. El, R. H. J. Grimshaw, N. F. Smyth, Transcriti-cal shallow-water flow past topography: finite-amplitudetheory, J. Fluid Mech. , 187 (2009).[34] J. G. Esler, J. D. Pearce, Dispersive dam-break and lock-exchange flows in a two-layer fluid, J. Fluid Mech. ,555 (2011).[35] M. A. Hoefer, Shock waves in dispersive Eulerian fluids,J. Nonlinear Sci. , 525 (2014).[36] T. Congy, A. M. Kamchatnov, and N. Pavloff, Disper-sive hydrodynamics of nonlinear polarization waves intwo component Bose-Einstein condensates, SciPost Phys. , 006 (2016).[37] M. A. Hoefer, G. A. El, A. M. Kamchatnov, Obliquespatial dispersive shock waves in nonlinear Schr¨odingerflows, SIAM J. Appl. Math. , 1352 (2017).[38] X. An, T. R. Marchant and N. F. Smyth, Dispersiveshock waves governed by the Whitham equation andtheir stability, Proc. Roy. Soc. London A , 20180278(2018).[39] A. M. Kamchatnov, Dispersive shock wave theory fornonintegrable equations, Phys. Rev. E , 012203 (2019).[40] S. K. Ivanov and A. M. Kamchatnov, Evolution of wavepulses in fully nonlinear shallow-water theory, Physics ofFluids , 057102 (2019).[41] T. Bienaim´e, private communication.[42] L. D. Landau and E. M. Lifshitz, Electrodynamics ofContinupous Media , (Pergamon, Oxford, 1984).[43] J.L. Coutaz and M. Kull, Saturation of the nonlinearindex of refraction in semiconductor-doped glass, J. Opt.Soc. Am. B , 95 (1991).[44] Y. S. Kivshar and G. P. Agrawal, Optical Solitons, FromFibers to Photonic Crystals (Academic Press, San Diego,2003).[45] R. W. Boyd,
Nonlinear Optics (Academic Press, SanDiego, 1992).[46] Y. Wang and M. Saffman, Experimental study of non-linear focusing in a magneto-optical trap using a Z-scan technique, Phys. Rev. A , 013801 (2004).[47] N. ˇSanti´c, A. Fusaro, S. Salem, J. Garnier, A. Picozzi,and R. Kaiser, Nonequilibrium Precondensation of Clas-sical Waves in Two Dimensions Propagating throughAtomic Vapors, Phys. Rev. Lett. , 055301 (2018).[48] Q. Fontaine, T. Bienaim´e, S. Pigeon, E. Giacobino, A.Bramati, and Q. Glorieux, Observation of the BogoliubovDispersion in a Fluid of Light, Phys. Rev. Lett. ,183604 (2018).[49] S. A. Akhmanov, A. P. Sukhorukov, R. V. Khokhlov,Self-focusing and diffraction of light in a Nonlinearmedium, Usp. Fiz. Nauk. , 19 (1967) [Sov. Phys. Usp. , 609 (1968)].[50] A. Sommerfeld, Partial Differential Equations in Physics (Academic Press, N. Y., 1949).[51] R. Courant and D. Hilbert,
Methods of Mathematical Physics , Vol. II (Interscience Publishers, N. Y., 1962).[52] G. S. S. Ludford, On an extension of Riemann’s methodof integration, with applications to one-dimensional gasdynamics, Proc. Cambridge Philos. Soc. , 499 (1952).[53] B. A. Dubrovin and S. P. Novikov, Hydrodynamics ofSoliton Lattices, Sov. Sci. Rev. Sect. C, Math. Phys. Rev. , 1 (1993).[54] A. M. Kamchatnov, New approach to periodic solutionsof integrable equations and nonlinear theory of modula-tional instability, Phys. Rep. , 199 (1997).[55] G. G. Stokes, The Outskirts of the Solitary Wave, in Mathematical and Physical Papers, vol. V, p. 163 (UP,Cambridge, 1905).[56] H. Lamb,
Hydrodynamics, (UP, Campridge, 1994).[57] G. A. El, V. V. Geogjaev, A. V. Gurevich, and A. L.Krylov, Decay of an initial discontinuity in the defocusingNLS hydrodynamics, Physica D87