Wave Breaking in Dispersive Fluid Dynamics of the Bose-Einstein Condensate
A. M. Kamchatnov




 
Institute of Spectroscopy,
  Russian Academy of Sciences, Troitsk, Moscow, 108840, Russia
Moscow Institute of Physics and Technology, Institutsky
  lane 9, Dolgoprudny, Moscow region, 141701, Russia



The problem of wave breaking during its propagation in the Bose-Einstein condensate
to a stationary medium is considered for the case when the initial profile at the
breaking instant can be approximated by a power function of the form (-x)^1/n.
The evolution of the wave is described by the Gross-Pitaevskii equation
so that a dispersive shock wave is formed as a result of breaking; this wave can be
represented using the Gurevich-Pitaevskii approach as a modulated periodic solution
to the Gross-Pitaevskii equation, and the evolution of the modulation parameters is
described by the Whitham equations obtained by averaging the conservation laws over
fast oscillations in the wave. The solution to the Whitham modulation equations is
obtained in closed form for n = 2,3, and the velocities of the dispersion shock
wave edges for asymptotically long evolution times are determined for arbitrary
integer values n > 1. The problem considered here can be applied for
describing the generation of dispersion shock waves observed in experiments with
the Bose-Einstein condensate.








Wave Breaking in Dispersive Fluid Dynamics
of the Bose-Einstein Condensate[JETP, 127, 903-911 (2018).]
    A. M. Kamchatnov
    
======================================================================================================






Contribution for the JETP special issue in honor of L. P. Pitaevskii's 85th birthday



§ INTRODUCTION


It is well known that with the disregard of viscosity
and dispersion effects, nonlinear waves experience
“breaking,” i.e., after a certain critical instant,
the formal solution to corresponding evolution equations
becomes multi-valued (see, for example, <cit.>).
In classical gas dynamics, this problem is eliminated by taking
into account weak dissipation effects so that instead of
the multi-valuedness domain, a shock wave (i.e., a narrow region
of transition from the flow with some values of parameters
characterizing the flow to a flow
with other values of parameters) appears in the solution.
The width of this transition region is proportional
to coefficients characterizing dissipative processes; in
real conditions, this width is usually of the order of the
mean free path of molecules in the gas. For this reason,
it can be assumed in the macroscopic theory that
this region is a discontinuity in the parameters of the
flow of the medium, and when the medium passes
through the discontinuity, the mass, momentum, and
energy conservation laws must hold. The theory of
shock waves formulated on this basis has been profoundly
developed and has found numerous applications (see, for
example, <cit.>).

In modern physics, however, flows of the medium
in which dissipation processes can be disregarded in
the first approximation are often considered, and then
nonlinear wave breaking is eliminated by taking into
account dispersion effects which lead to formation of
dispersive shock waves (DSWs) (i.e., the evolving
regions of the nonlinear flow of the medium) instead of
the multi-valued domains. Such effects were studied
for the first time in the theory of undular bores in a
shallow water flow (see, e.g., <cit.>), and the general
nature of this phenomenon was realized by Sagdeev <cit.>,
who indicated that wave breaking in dispersive
wave systems leads to the formation of an extended
wave structures connecting different states of the flow
like a transition in a shock wave connects different
states of the medium flow with predominant dissipation.
In typical cases, a dispersive shock wave occupies
a spatial region expanding with time so that this wave
is a sequence of solitons at one of its edge and degenerates
into a small-amplitude harmonic wave propagating with the
corresponding group velocity at the
other edge. The main theoretical approach to the
description of DSWs was proposed in the classical
work by Gurevich and Pitaevskii <cit.> based on the
Whitham theory of modulation of nonlinear waves <cit.>.
In this approach, a DSW is represented in the form of
a modulated periodic solution to the corresponding
nonlinear wave equation, and the slow evolution of the
modulation parameters obeys the Whitham equations
obtained by averaging of the conservation laws over
fast oscillations of physical variables in the wave.
Gurevich and Pitaevskii considered two typical problems
of wave breaking, the evolution of the wave is
described by the Korteweg-de Vries (KdV) equation.
First, a complete analytic solution of the discontinuity
decay was obtained in the case when the initial distribution
at the breaking instant has a sharp jump. Second, they found
the main characteristics of the DSW
in vicinity of the breaking point when the initial distribution
is described by a cubic parabola. Later, Potemin <cit.>
obtained a full analytic solution to this problem
(see also <cit.>). The Gurevich and
Pitaevskii approach to the DSW theory was developed
further and was extended to other equations (see, for
example, review <cit.>).

One of important applications of the DSW theory
is the dynamics of the Bose-Einstein condensate,
which is described by the Gross-Pitaevskii equation <cit.>;
for simplicity, we write here this equation in
the standard dimensionless form for a 1D flow of
the condensate:

    iψ_t+1/2ψ_xx-|ψ|^2ψ=0,

where ψ is the “wave function” of the condensate
flowing along the x coordinate; we assume that the
interaction between atoms is repulsive, which ensures
the stability of its homogeneous state. The theory of
Eq. (<ref>) was considered in a huge number of publications.
In particular, its solution in the form of a dark
soliton was obtained in <cit.>, and periodic solutions
were obtained, for example, in <cit.>. The integrability
of Eq. (<ref>) by the inverse scattering transform method
was established in <cit.> and this approach was used in
<cit.> for deriving the modulation equations.
Finally, the problem of the initial discontinuity decay
was analyzed in <cit.>, and typical wave breaking was
studied in <cit.>. The apparatus developed in these
works was applied to the description of the DSW
dynamics in the Bose-Einstein condensate.

The DSWs in the condensate were observed experimentally for the first time
in <cit.>, where a shock
wave was formed under the action of laser radiation
repelling the condensate. The interpretation of these
observations as DSWs was reported in <cit.>,
and this experiment was described in <cit.>
using the theory formulated in <cit.>,
where it was assumed that a DSW is
formed as a result of emergence of a discontinuity.
Although such a discontinuity can be formed when the
flow of the condensate is induced by a piston moving
at a constant velocity <cit.>,
such a case is nevertheless
quite specific, and the wave propagating to the bulk of
a stationary medium in typical situations breaks from
a profile with a certain root singularity rather than a
sharp discontinuity. For instance, a singularity in the
form of a square root appears in the case of the uniformly
accelerated motion of the piston <cit.>, and it is
clear even from this example that the actual flow of the
condensate may have a quite arbitrary singularity at
the instant of wave breaking. Here, we consider the
DSW formation during wave breaking with an initial
singularity of the type (-x)^1/n. A detailed theory will be
developed for n=2 and n=3, and important DSW
characteristics (such as the laws of motion of its edges)
will be obtained for an arbitrary integer n>1.  The theory
developed here forms the basis for describing quite
general forms of condensate flow with wave breaking.



§ GUREVICH-PITAEVSKII METHOD


Let us first write the basic relations of the Gurevich-Pitaevskii 
method in the DSW theory as applied
to the dynamics of the Bose-Einstein condensate,
which obeys the Gross-Pitaevskii equation (<ref>).
It is convenient to represent the periodic solutions to this
equation in terms of more transparent physical variables by
performing the substitution

    ψ(x,t)=√(ρ(x,t))exp(i∫^x u(x',t)dx'),

so that after the separation of the real and imaginary
parts, we obtain the system

    ρ_t+(ρ u)_x=0,
       u_t+uu_x+ρ_x+[ρ_x^2/8ρ^2
       -ρ_xx/4ρ]_x= 0.

Here, ρ(x,t)=|ψ(x,t)|^2  is the condensate density and u(x,t)
is the condensate flow velocity. The periodic
solution can be written in the form

    ρ =   1/4(_4-_3-_2+_1)^2+ (_4-_3)×
       ×(_2-_1)  sn^2(√((_4-_2)(_3-_1)) θ,m)  ,


    u=V - j/ρ  ,

where

    j   =1/8 (-λ_1 - λ_2 + λ_3 +
    λ_4)×
       ×(-λ_1 + λ_2 - λ_3 + λ_4)
    (λ_1 - λ_2 - λ_3 + λ_4),
    θ   =x-Vt,     V=1/2∑_i=1^4_i,
    
    m   =(_2-_1)(_4-_3)/(_4-_2)(_3-_1).

This solution depends on four parameters _1≤_2≤_3≤_4,
in terms of which the main characteristics of
the wave can be expressed. In particular, the wavelength is given by

    L= 2K(m)/√((_4-_2)(_3-_1)),

where K(m)  is the complete elliptic integral of the first
kind. In the limit _3→_2,  when m→1 and L→∞,
the periodic wave is transformed into the soliton solution

    ρ=ρ_0 - a_s/cosh^2(√(a_s)(x-V_st)) ,

where background density ρ_0, along which the dark
soliton propagates, its amplitude a_s  and velocity V_s are
given by

    ρ_0   =1/4(_4 - _1)^2,     a_s=(_4 - _2)(_2 - _1) ,
    
        V_s   =1/2(_1+2_2+_4) .

In the opposite limit _3→_4, when m→0, the wave
amplitude tends to zero, and it is transformed into a
linear harmonic wave propagating over the background with constant density ρ_0.

DSW parameters _i become slow functions of x and t, which change little
in one wavelength L. Therefore, we can average the conservation laws for
Eq. (<ref>) over fast oscillations in the wave and obtain as
a result the Whitham equations for modulation parameters _i.
These equations can be written in the form

    ∂_i/∂
    t+v_i()∂_i/∂ x=0,      i=1,2,3,4,

where the characteristic velocities are given by

    v_i()   =(1-L/∂_iL∂_i)V ,    i=1,2,3,4   ,
    ∂_i   ≡∂/∂_i .

The substitution of expressions (<ref>) and (<ref>) into these
formulas gives the following expressions for velocities:

    v_1   =1/2∑_i
    -(_4-_1)(_2-_1)K/(_4-_1)K-(_4-_2)E,
    
    v_2   =1/2∑_i
    +(_3-_2)(_2-_1)K/(_3-_2)K-(_3-_1)E,
    
    v_3   =1/2∑_i
    -(_4-_3)(_3-_2)K/(_3-_2)K-(_4-_2)E,
    
    v_4   =1/2∑_i
    +(_4-_3)(_4-_1)K/(_4-_1)K-(_3-_1)E,

where  E=E(m) is the elliptic integral of the second
kind. Variables _i are known as Riemann invariants of
the system of the Whitham modulation equations. We
will also need the limiting expressions for these velocities
at the DSW edges. At the soliton edge, where _3=_2 and m=1, we have

    v_1= 3/2_1 + 1/2_4  ,   
             v_4= 3/2_4 + 1/2_1, 
       v_2=v_3=1/2(λ_1+ 2λ_2 + λ_4),

while at the small-amplitude edge for _3=_4 and m=0 we have

    v_1= 3/2_1 + 1/2_2  ,   
    v_2= 3/2_2 + 1/2_1 , 
        v_3=v_4=2λ_4 +
    (λ_2-λ_1)^2/2(λ_1+_2-2_4)  ,

(we will not need the expressions for the analogous
limit _1=_2).

In the generalized hodograph method <cit.> the
solutions to Eqs. (<ref>)
are sought in the form

    x-v_i()t=w_i(),   i=1,2,3,4,

where v_i() are velocities (<ref>) and w_i() are the sought
functions. If these functions have been determined, x=x() and t=t()
turn out to be the functions of parameters _i. Since these functions
must be inverted and the modulation parameters must
become functions _i=_i(x,t),
the functions w_i cannot be independent of one another. Differentiating
Eq. (<ref>) with respect to _j, j≠ i, and eliminating t from
all pair combinations of the resultant relations, we
arrive at the system of the Tsarev equations

    1/w_i-w_j w_i/_j=1/v_i-v_j v_i/_j,
           i≠ j.

In view of their symmetry in v_i and w_j, it is natural to
seek their solution in the form analogous to (<ref>) (see <cit.>):

    w_i()=(1-L/∂_iL∂_i)W ,    i=1,2,3,4 .

Then Eqs. (<ref>) are transformed into the system of
Euler-Poisson equations (i≠ j)

    ^2W/_i_j-1/2(_i-_j)( W/_i-
         W/_j)=0.

For our purposes, it is sufficient to know the set of
solutions obtained from the generating function

    W=^2/√(∏_i=1^4(-_i)),

that satisfies Eqs. (<ref>) for any . The expansion of
function (<ref>)
in inverse powers of   gives

    W=∑_k=0^∞W^(k)/^k=1+W^(1)1/
        +W^(2)1/^2+…,

where W^(k)=W^(k)(_1,_2,_3,_4) are the required
particular solutions to system (<ref>). As a result, using
expression (<ref>)  we obtain the set of functions w_i^(k),
which give the solutions to Whitham equations (<ref>):

    w_i^(k)=W^(k)+(2v_i-s_1)_iW^(k),

so that w_i^(0)=1, w_i^(1)=v_i.  The Euler-Poisson equation
is linear in W, like expressions (<ref>) that are linear
in W^(k); therefore, any of their linear combinations also
gives the solution

    x-v_i()t=∑_k=0^n A_kw_i^(k)(),   i=1,2,3,4,

where the number of terms n and constant coefficients A_k are
chosen in accordance with the conditions of the problem.

Let us now show that the above expression of the
DSW theory in the Gurevich and Pitaevskii approach
make it possible to solve the problem of the DSW formation during
wave breaking in the Bose-Einstein condensate.



§ DISPERSIONLESS LIMIT


Until the instant of breaking, the distributions of
density ρ(x,t) and flow velocity
u(x,t) are smooth functions of spatial coordinate x. Moreover, even after
breaking, the DSW occupies a finite spatial region and
its edges at the matching points with the smooth distributions
must be determined as a part of the solution
of the wave breaking problem. In the case of quite
smooth functions ρ and u  the terms with a large number
of derivatives in system (<ref>) can be omitted, which
means the disregard of the dispersion effects; therefore,
the evolution of smooth distributions can be
described by the dispersionless limit equations

    ρ_t+(ρ u)_x=0,   u_t+uu_x+ρ_x=0.

These equations coincide with the “shallow water”
equations (see <cit.>), equivalent to the gas dynamic
equations with adiabatic exponent γ=2. Therefore, their
solutions can be obtained using well-known classical methods.

Equations (<ref>) can be transformed to diagonal
form by introducing the Riemann invariants

    _±=u/2±√(ρ),

so that

    _+/ t+v_+(_+,_-)_+/ x=0,
       _-/ t+v_-(_+,_-)_-/ x=0,

where

    v_+=3/2_++1/2_-,   v_-=1/2_++3/2_-,

here, ρ and u can be expressed in terms of _± by the
formulas

    ρ=1/4(_+-_-)^2,   u=_++_-.


Generally, both Riemann invariants _± are functions of x and t.
We are interested, however, in the
problem in which a wave propagates to the bulk of the
condensate at rest with constant density ρ_0. It is known <cit.>,
that only a flow in the form of a simple wave, in
which one of the Riemann invariants is constant, can
border on such a state of the gas. Assuming for definiteness that the wave
propagates to the right, we can
conclude that Riemann invariant _- must be constant
and, hence, must have the same value as in the stationary medium bordering the wave:

    _-=-√(ρ_0).

Then the second equation in (<ref>) is satisfied automat-
ically, while the first equation is transformed into the
Hopf equation

    _+/ t+(3/2_+-1/2√(ρ_0))
      _+/ x=0

with the well-known general solution

    x-(3/2_+-1/2√(ρ_0))t=w(_+).

If function w(_+) is known, this solution is defined by
the dependence _+=_+(x,t),
and this function must
be joined at the boundary with the stationary condensate with value
_+=√(ρ_0) in this spatial region.

We are interested in the situation when the smooth
solution for _+ at the breaking instant tends to its
boundary value √(ρ_0) as a root function of x. Choosing
the coordinate system and its origin so that breaking of
Riemann invariant _+ occurs at instant t=0 at the origin x=0,
we obtain the dependence

    x-(3/2_+-1/2√(ρ_0))t=-(_+-√(ρ_0))^n,

where, to simplify calculations, the units of measurements of length and
time are chosen so that the coefficient on the right-hand side be equal to unity.
Therefore, the dependence of _+ on x for t<0
has no singularities, while, at t=0 the root singularity appears,

    _+=√(ρ_0)+(-x)^1/n,

and the dependence of _+ on x for t>0
becomes multi-valued (see Fig. <ref>).




§ DISPERSIVE SHOCK WAVE


At instant t the DSW occupies the spatial region

    x_-(t)≤ x≤ x_+(t),

matching to the smooth solution (<ref>) at its boundary
point. Comparison of velocity v_+ in expression (<ref>)
with limiting expression (<ref>) shows that the DSW is
transformed into expression (<ref>) for _4=_+
at boundary x=x_-(t); coefficients A_k in this case must be
chosen so that the right-hand side of relation (<ref>)
with i=4 be equal to the right-hand side of relation (<ref>).
Further, the solution to the Whitham equations is transformed
into a harmonic wave at the small-amplitude
edge x=x_+(t), if _2=_+=√(ρ_0) along the entire DSW,
and we have _3=_4 at point x_+(t). Since _-=-√(ρ_0)
at both edges of the SDW, the condition of matching of _1 to _-
at the DSW edges can be satisfied by setting _1=-√(ρ_0)
along the DSW. Therefore, Whitham equations (<ref>) with i=1,2
are satisfied by constant solutions

    _1=-√(ρ_0),    _2=√(ρ_0),

and only two Riemann invariants _3 and _4, which
satisfy the boundary conditions

    _4=_+,  _3=_2,   w_4=-(_4-_2)^n,
       (m=1),   x=x_-(t)

and

    _3=_4,     (m=0),   x=x_+(t)

vary along the DSW. These conditions define the solution completely.
As a result, the dependence of the
Riemann invariants on coordinate x for a fixed value of t
has the form shown in Fig. <ref>.
It should be noted that
waves with two variable Riemann invariants were
called quasi-simple and were studied for the first time
in <cit.> in the theory of the KdV equation. Taking
relations (<ref>) into account, we can find the first three
coefficients in Eq. (<ref>) in the form

    W_0=1,   W_1=1/2(_3+_4),
        W_3=1/8(4_2^2+3_3^2+2_3_4+3_4^2),
        W_3=1/16(_3+_4)(4_2^2+5_3^2-2_3_4+5_4^2),

the knowledge of these coefficients is sufficient for
analyzing typical cases with n=2 and n=3.




 §.§ The case with n=2


Since formulas (<ref>) are polynomial, quadratic
function -(_4-_2)^2 can be written in the form of a
linear combination of the first three expressions (<ref>)
with coefficients

    A_0=-4/15_2^2,   A_1=16/15_2,   A_2=-8/15.

Then formulas (<ref>) with i=3,4 and with these values
of the coefficients define implicitly the dependencies
of _3 and _4 on x and t, which solves in principle the
problem in this particular case. At the soliton edge,
these formulas are transformed into

    x-1/2(_4+_2)t=-1/5(_4-_2)^2,
        x-1/2(3_4-_2)t=-(_4-_2)^2,

which gives the relation between t and _4 at this
boundary:

    t=4/5(_4-_2),    _4=_2+5/4t.

Substituting the resultant value of _4 into any relation
from (<ref>) we obtain the law of motion of the soliton edge:

    x_-(t)=_2t+5/16t^2.

At the small-amplitude edge for _3=_4 both formulas(<ref>) with i=3,4
are transformed into the same relation

    x_+-(2_4-_2^2/_4)t=
      -8/15(_4-_2)^2(3+2_2/_4).

This edge moves with the group velocity corresponding to wavenumber
k=2π/L with L=2√(_4^2-_2^2), which gives for the Bogoliubov dispersion law
ω=k√(_2^2+k^2/4) the expression dω/dk=2_4^2-_2^2/_4.
Therefore, the differentiation of expression (<ref>) with
respect to t with account of dx_+/dt=dω/dk
determines the dependence of t on the value of _4 at this
boundary. Introducing parameter y=_4/_2 instead of
_4, we can write this dependence in the form

    t=16_2/15·(y-1)(3y^2+y+1)/2y^2+1,   y≥1,

and its substitution into relation (<ref>) gives

    x_+=8_2^2/15·(y^2-1)(6y^2-1)/2y^2+1.

Formulas (<ref>) and (<ref>) define parametrically the law
of motion x_+=x_+(t) of the small-amplitude edge. For
t∼_0y≫_2^2=ρ_0 this law of motion asymptotically
takes the simple form

    x_+(t)≈5/8t^2,   t≫ρ_0.

It should be noted that analogous expressions
obtained by solving the problem of motion of the condensate under
the action of a uniformly accelerated
piston <cit.> can be transformed to the expression
obtained above after the transfer of the breaking point
to the origin and subtracting the breaking time from t.



 §.§ Case with n=3


In this case, the calculations are performed analogously.
The right-hand sides of formulas (<ref>) with
i=3,4 now contain function W_3 and condition (<ref>) gives
the values of the coefficients

    A_0=-8/35_2^3,   A_1=-8/7_2^2,
        A_2=48/35_2,   A_3=-16/35.

At the soliton edge, solution (<ref>) is transformed into

    x-1/2(_4+_2)t=-1/7(_4-_2)^3,
        x-1/2(3_4-_2)t=-(_4-_2)^3,

which gives

    t=6/7(_4-_2)^2,    _4=_2+√(7t/6).

The substitution of this relation into (<ref>) gives the law
of motion of the soliton edge:

    x_-=_2t+1/3√(7/6) t^3/2.

At the small-amplitude edge, formulas  (<ref>) with i=3.4 are transformed into

    x_+-(2_4-_2^2/_4)t=
      -16/35(_4-_2)^3(4+3_2/_4).

The matching condition for the law of motion with the
group velocity of a linear wave gives

    t=48_2^2/35·(y-1)^2(4y^2+2y+1)/2y^2+1,   y≥1,

and the substitution into expression (<ref>) gives

    x_+=16_2^3/35·(y-1)^2(16y^3+14y^2-4y-5)/2y^2+1.

For asymptotically long times t∼_0^2y^2≫_2^2=ρ_0^2 we
obtain

    x_+(t)≈1/3√(35/6) t^3/2,   t≫ρ_0^2.

The laws of motion of the DSW edges as functions of
time for n=2,3  is shown in Fig. <ref>.
For short times, the
motion with the velocity of sound √(ρ_0),
prevails, while,
for long times, a transition to asymptotic laws occurs.


General formulas (<ref>) together with the specific
expressions for functions w_i^(k), obtained above with the help
of formulas (<ref>), (<ref>) and coefficients (<ref>) (n=2)
and (<ref>) (n=3), make it possible to calculate
_3 and _4 as functions of x and t, so that their substitution into expressions
(<ref>) and (<ref>)  gives the density and
velocity distribution profiles in a DSW (Fig. <ref>). The
envelopes of the condensate density in the DSW are
calculated by the formulas

    ρ_max=1/4(_4-_3+2√(ρ_0))^2,  ρ_min=1/4(_4-_3-2√(ρ_0))^2.



Although the formulas are complicated with
increasing exponent n,
important DSW characteristics
such as the laws of motion of the edges can be determined without
detailed analysis of the complete solution (at least,
in the asymptotic limit t≫ρ_0^n).
In the next sections, we will consider this problem.



§ LAW OF MOTION OF THE SOLITON EDGE


Combining relations (<ref>) with the limiting expression of formula (<ref>),
we can write the boundary condition at the soliton edge in the form of a differential
equation for the function W=W(-_2,_2,_2,_4),
which depends only on _4:

    W+2(_4-_2)dW/d_4=-(_4-_2)^n,

the solution to this equation is

    W(-_2,_2,_2,_4)=-1/2n+1(_4-_2)^n,

where the integration constant is chosen so that time t
in subsequent formulas tends to zero for _4→_2. Since w_3=W at this boundary,
formulas (<ref>) give

    x-1/2(_4+_2)t=-1/2n+1(_4-_2)^n,
        x-1/2(3_4-_2)t=-(_4-_2)^n,

which is in conformity with relations (<ref>) (n=2) and (<ref>) (n=3). This gives

    t=2n/2n+1(_4-_2)^n-1,

and

    x_-(t)=_2t+n-1/2n(2n+1/2n)^1/n-1t^n/n-1.

These expressions generalize the formulas obtained
above to arbitrary integer values of n>1.

It should be noted that, in fact, we can find the law
of motion of the soliton edge for an arbitrary monotonic dependence
of the initial distribution of Riemann invariant _+ of the form
w=w(_+-√(ρ_0)), by resorting to the considerations used in
<cit.> for deriving
the law of motion of the small-amplitude edge for
wave breaking in the theory of the KdV equation.
Indeed, Whitham equations (<ref>) with i=3,4
in the classical hodograph method are transformed into linear differential
equations for functions x=x(_3,_4) and t=t(_3,_4),
one of which for  _3→_2 becomes

    x/_4-1/2(_2+_4) t/_4=0.

On the other hand, the solution at this boundary must
match to the smooth solution, which gives

    x-1/2(3_4-_2)t=w(_4-_2).

The differentiation of this relation with respect to _4
leads to one more differential equation

    x/_4-3/2t-1/2(3_4-_2) t/_4=
       w/_4.

Eliminating x/_4 from Eqs. (<ref>) and (<ref>), we obtain
the differential equation for t, the integration of which
gives (z=_4-_2)

    t=t(z)=-z^-3/2∫_0^z z^1/2dw/dzdz,

and, hence,

    x_-=x_-(z)=(32z+_2)t(z)-w(z).

These formulas specify the parametric dependence x_-=x_-(t).



§ LAW OF MOTION OF THE SMALL-AMPLITUDE EDGE


Formulas (<ref>) (n=2) (<ref>) and (n=3) have a simple structure
leading to the assumption that for m→0 (_3→_4) and integer n functions w_3 and w_4
must pass to the right-hand side of the relation

    x_+-(2_4-_2^2/_4)t=
      A_n(_4-_2)^n(n+1+n_2/_4).

Let us prove this formula in the asymptotic limit t→∞ (_4→∞),
when the terms with _2/_4 can be
neglected. We note that in the limit _2→0
generating function (<ref>) can be reduced to the generating function
of the Legendre polynomials (see, for example, <cit.>)

    W   ≈1/√(1-2(_3+_4/2√(_3_4))
      √(_3_4)/+(√(_3_4)/)^2)=
       =∑_nP_n(_3+_4/2√(_3_4))(_3_4)^n/2/^n,

i.e.,

    W_n≈ (_3_4)^n/2P_n(_3+_4/2√(_3_4)).

Using the recurrent formula for the derivative of the
Legendre polynomial (see <cit.>) we can easily prove that

    W_n/_4=n/_4-_3(W_n-_3W_n-1)

in this approximation. To evaluate function w_4^(n) in the
limit m→0 we will prove that the following relation
holds for  _3→_4:

    W_n-_3W_n-1≈(_4-_3)·1/2_4^n-1.

For this purpose, we note that the argument of the
Legendre polynomial in expression (<ref>), which is the
ratio of the arithmetic mean to the geometric mean,
attains its maximal value for _3=_4
and, hence, is
quadratic in the small difference _4-_3=, so that P_n(1)=1, i.e.,
W_n(_4,_4)=_4^n. Consequently, in the
first order in  we obtain

    W_n(_4-,_4)-(_4-)W_n-1(_4-,_4)≈
       ≈(_4-)^n/2_4^n/2-_4(_4-
      )^n-1/2_4^n-1/2+_4^n-1≈
       ≈1/2_4^n-1.

In addition, we note that for _2→0 only the highest
term with  k=n is left in the sum in expression (<ref>)
since coefficients A_k for k<n  contain powers of _2 as
factors. Therefore, with account of relation.
2(V-v_4)≈-2_4 we obtain

    w_4≈ A_n(n+1)_4^n.

For determining w_4 completely for m=0 it remains
to find A_n, which can easily be done using the
expansion of the higher term in the Legendre polynomial into a series
for _3→_2=0 (see <cit.>):

    W_n   ≈(_2_4)^n/2·(2n)!/2^n(n!)^2(_2+_4/2√(_2_4))^n≈
       ≈(2n)!/4^n(n!)^2 _4^n.

Then the condition of matching to the smooth solution gives

    w_4≈ A_n·(2n+1)!/4^n(n!)^2 _4^n =-_4^n,

i.e.

    A_n=-4^n(n!)^2/(2n+1)!.

Differentiating the expression

    x-2_4t=w_4=A_n(n+1)_4^n

with respect to _4 provided that x/_4=0  for a fixed t,
which is equivalent to the matching condition with
the group velocity at the small-amplitude edge, we
obtain

    t=-1/2n(n+1)A_n_4^n-1.

Substituting _4 obtained from this expression into (<ref>)
we obtain the law of motion of the small-amplitude edge in the asymptotic regime:

    x_+=2(n-1)/n((2n+1)!/2^2n-1n(n+1)(n!)^2)^1/n-1t^n/n-1.

This formula naturally reproduces the above asymptotic laws
(<ref>) for n=2 and (<ref>) for n=3.

Expression (<ref>) confirms the validity of formula (<ref>)  in the limit
_2→0. Assuming that this formula
also holds for a finite _2,
we obtain the law of motion
of the small-amplitude boundary in parametric form:

    t   =nA_n_2^n-1(y-1)^n-1/2y^2+1[(n+1)y^2+(n-1)y+1],
    
      x   =A_n_2^n(y-1)^n-1/2y^2+1[2(n^2-1)y^3+
       +2(n^2-n+1)y^2-(n-1)^2y-n^2+n+1].


It should be noted that in contrast to the theory of
the KdV equations, the self-similar regime of motion
of the boundaries is realized only asymptotically for
long times
t≫ρ_0^n-1, when the velocity of motion is
much higher than the velocity of sound in the background distribution.
However, the limiting transition to ρ_0→0 is impossible in the expressions describing
the wave profile since the magnitude of m in elliptic
functions vanishes for _2=-_1→0.



§ CONCLUSION


Thus, the approach developed by Gurevich and
Pitaevskii makes it possible to analyze in detail the
process of DSW formation during wave breaking in the
Bose-Einstein condensate, the dynamics of which
obeys the Gross-Pitaevskii equation. The developed
theory is applicable to the initial stage of the process,
in which the smooth part of the profile can be treated
as a monotonic function of the coordinate. It should
be noted, however, that the theory of quasi-simple
waves must also describe the asymptotic stage of the
evolution of a finite-duration pulse since, analogously
to the simple wave theory, the initial pulse splits with
time into two pulses, in each of which two of four Riemann 
invariants again remain constant. Therefore, the
Gurevich-Pitaevskii approach supplemented with the
generalized hodograph method and modern method
for deriving the Whitham modulation equations
remains a powerful tool for investigating dispersion
shock waves, which are of considerable interest for
modern nonlinear physics.


  I thank M. Isoard, S. K. Ivanov and N. Pavloff for useful
  discussions.





LL-6 L. D. Landau, E. M. Lifshitz, Fluid Mechanics, Pergamon,
 Oxford (1987).

CF-50 R. Courant,  K. O. Friedrichs, Supersonic Flow and Shock Waves,
Interscience Publishers,  New York (1948).

bl-54 T. B. Benjamin, M. J. Lighthill, Proc. Roy. Soc. London, A 224, 448 (1954).

sagdeev R. Z. Sagdeev, Cooperative phenomena and shock waves in collisionless
plasmas, Rev. Plasma Phys. 4,  23 (1966).

gp-73 A. V. Gurevich and L. P. Pitaevskii,
  Zh. Eksp. Teor. Fiz. 65, 590 (1973) [Sov. Phys. JETP 38,
  291 (1974)].

whitham G. B. Whitham, Proc. Roy. Soc. London, A 283, 238 (1965).

potemin G. V. Potemin, Russian Math. Surveys, 43, 252 (1988).

kamch A. M. Kamchatnov, Nonlinear Periodic Waves and Their Modulations.
An Introductory Course, World Scientific, Singapore (2000).

eh-16 G. A. El, M. A. Hoefer, Physica D 333, 11 (2016).

gross E. P. Gross.

 Nuovo Cimento, 20, 454 (1961).

pit L. P. Pitaevskii,  
Zh. Eksp. Teor. Fiz. 40, 646 (1961) [Sov. Phys. JETP 3, 451.

tsuzuki-1971 T. Tsuzuki,

J. Low Temp. Phys., 4, 441 (1971).

gk-87  A.V. Gurevich and A.L. Krylov, 

Zh. Eksp. Teor. Fiz. 92, 1684 (1987) [Sov. Phys. JETP,  65,  944 (1987)].

zs-73 V. E. Zakharov and A. B. Shabat, Zh. Eksp. Teor. Fiz., 64, 1627 (1973)
[Sov. Phys. JETP, 37, 823 (1973).

fl-86  M. G. Forest and J. E. Lee,


in Oscillation Theory,
Computation, and Methods of Compensated Compactness, Eds. C. Dafermos et al,
IMA Volumes on Mathematics and its Applications 2, (Springer, N.Y., 1986).

pavlov-87 M. V. Pavlov, 
Teor. Mat. Fiz. 71, 351 (1987) [Theor. Math. Phys., 71, 584 (1987)].

eggk-95  G. A. El, V. V. Geogjaev, A. V. Gurevich, A. L. Krylov,

Physica D, 87, 186 (1995).

kku-02 A. M. Kamchatnov, R. A. Kraenkel, B. A. Umarov,

Phys. Rev. E,  66,  036609 (2002).

cornell E. A. Cornell, Talk at NATO Advanced Workshop “Nonlinear
Waves: Classical and Quantum Aspects,” Lisbon, 2003.

hoefer M. A. Hoefer, M. J. Ablowitz, I. Coddington,
E. A. Cornell, P. Engels, V. Schweikhard,

 Phys. Rev. A   74, 023623 (2006).

kgk-04  A. M. Kamchatnov,  A. Gammal, R. A. Kraenkel,


Phys. Rev. A  69, 063605 (2004).

hae-08 M. A. Hoefer, M. J. Ablowitz, P. Engels,
Phys. Rev. Lett. 100, 084504 (2008).

kk-10 A. M. Kamchatnov, S. V. Korneev,
Zh. Eksp. Teor. Fiz., 137, 191 (2010) [JETP, 110, 170 (2010)].

tsarev S. P. Tsarev, Math. USSR Izv. 37, 397 (1991).

gke-92  A. V. Gurevich, A. L. Krylov and G. A. El,
  Zh. Eksp. Teor. Fiz. 101, 1797 (1992) [Sov. Phys. JETP
  74, 957 (1992)].

wright O. Wright, Commun. Pure Appl. Math., 46, 421 (1993).

tian F. R. Tian, Commun. Pure Appl. Math., 46, 1093 (1993).

gkm-89 A. V. Gurevich, A. L. Krylov, and N. G. Mazur,
  Zh. Eksp. Teor. Fiz. 95, 1674 (1989) [Sov. Phys. JETP 68, 966 (1989)].

ww E. T. Whittaker and G. N. Watson, A Course of Modern Analysis,
CUP, Cambridge (1927).



