PHYSICAL REVIEW E 109, 015102 (2024)

Undular bore theory for the modified Korteweg–de Vries–Burgers equation

1 1,2,3
L. F. Calazans de Brito and A. M. Kamchatnov
1
Higher School of Economics, 20 Myasnitskaya Ulitsa, Moscow 101000, Russia
2
Institute of Spectroscopy, Russian Academy of Sciences, Troitsk, Moscow 108840, Russia
3
Skolkovo Institute of Science and Technology, Skolkovo, Moscow 143026, Russia

(Received 22 August 2023; accepted 19 December 2023; published 16 January 2024)

We consider nonlinear wave structures described by the modified Korteweg–de Vries equation, taking into
account a small Burgers viscosity for the case of steplike initial conditions. The Whitham modulation equa-
tions are derived, which include the small viscosity as a perturbation. It is shown that for a long enough time
of evolution, this small perturbation leads to the stabilization of cnoidal bores, and their main characteristics
are obtained. The applicability conditions of this approach are discussed. Analytical theory is compared with
numerical solutions and good agreement is found.

DOI: 10.1103/PhysRevE.109.015102

I. INTRODUCTION The Whitham modulation equations for the mKdV case were
derived in Ref. [18], but their application to the theory of
The modified Korteweg-de Vries (mKdV) equation
dispersive shock waves turned out to be quite a difficult task
ut − 6αu2 ux + uxxx = 0 (1) even in the case of an initial discontinuity problem. The reason
for this difficulty is that the mKdV equation is not genuinely
appeared first in the study of the famous KdV equation nonlinear. (This notion was introduced by Lax in Ref. [19] for
ut + 6uux + uxxx = 0 (2) hyperbolic systems of first-order partial differential equations
and it plays an important role in the classification of wave
related to Eq. (1) by the Miura transformation [1]. The ex- structures evolving from initial discontinuities in dispersive
istence of such a transformation allowed the pioneers of the nonlinear systems; see, e.g., Ref. [20].) This means that in
inverse scattering transform method to discover this method the dispersionless approximation, the nonlinear velocity 6αu2
[2–4] for the KdV equation, and it was extended later to many has an extremal (minimal for α > 0) value at u = 0, whereas
other equations, including the mKdV equation [5,6] (see also, in the case of the genuinely nonlinear KdV equation, the
e.g., the books in [7–9] and references therein). The mKdV nonlinear velocity 6u is everywhere a monotonic function of
equation is almost as widely used in physical applications as the wave amplitude u. As a result, in the KdV case an initial
the KdV equation. Actually, the Gardner equation discontinuity can only evolve into two different structures
ut + 6βuux − 6αu2 ux + uxxx = 0, (3) (rarefaction waves or cnoidal undular bores), whereas in the
mKdV case an initial discontinuity evolves into eight different
combining the nonlinear terms of the KdV and mKdV equa- wave structures depending on the parameters of the initial
tions, can be transformed into Eq. (1) by a simple change of jump of u. Some particular results in this direction were ob-
variables. In addition, in physical applications, it often hap- tained in Ref. [21] and the full solution was given in Ref. [22]
pens that the coefficient β is very small and can be neglected, in the context of the Gardner equation (3).
so Eq. (3) reduces directly to the equation. The Gardner In Gurevich-Pitaevskii theory, dispersive shock waves
equation and its simplified mKdV version find applications are wave structures that expand with time, so in initial
to the theory of nonlinear waves in stratified fluids, for ex- discontinuity-type problems, the change of modulation pa-
ample, for the description of large-amplitude internal waves rameters per unit length decreases with time and can become,
[10–12]. at large enough time, smaller than some other physical pa-
One of the most important and universal phenomena in rameters that were neglected in the derivation of Eq. (1) or
nonlinear physics is the formation and evolution of dispersive (2). For such large values of time, the neglected effects must
shock waves (see, e.g., review articles in [13,14] and refer- be taken into account in the modulation theory. For example,
ences therein). They are called undular bores in water wave small dissipation stops the infinite expansion of undular bores
physics and they were observed in both surface and inter- and their length is stabilized at some value inversely propor-
nal waves. Their theory was initially developed by Gurevich tional to the viscosity coefficient in accordance with the early
and Pitaevskii [15], who represented such structures as mod- ideas of Refs. [23,24] about the structure of undular bores in
ulated nonlinear periodic waves governed by the Whitham water-wave physics and plasma. The corresponding modified
modulation equations [16,17]. They gave two typical exam- Whitham equations for the KdV theory with weak Burgers
ples of solutions that describe dispersive shock waves: the dissipation were derived in Refs. [25,26] and were applied in
evolution of an initial discontinuity and the formation of a these papers to the description of stationary dispersive shocks
shock after generic wave breaking for the KdV equation case. whose characteristic length is defined by the small viscosity

2470-0045/2024/109(1)/015102(11) 015102-1 ©2024 American Physical Society

DE BRITO AND KAMCHATNOV PHYSICAL REVIEW E 109, 015102 (2024)

coefficient γ in the KdV-Burgers equation regions an evident self-similar solution of Eq. (6),
⎧
⎪
⎨u− , x < xl
ut + 6uux + uxxx = γ uxx . (4)
u(x, t ) = ± x
−6αt
, xl < x < xr (7)
⎪
⎩
The extension of this theory to the mKdV-Burgers (mKdVB) u+ , x > xr .
equation
Obviously, such a solution exists only if both boundary values
u± are lying in either of the monotonicity intervals 0 < u+ <
ut − 6αu ux + uxxx = γ uxx
2
(5)
u− or 0 > u+ > u− . In both cases, these rarefaction waves
propagate to the left with the left edge speed sl = −6αu− 2
was discussed qualitatively in Ref. [27]; however, the modi-
smaller than the right edge speed sr = −6αu+ 2
. If the bound-
fied Whitham equations were not obtained for this case and
ary values u± are lying in the different monotonicity intervals,
the quantitative theory was not developed. The main aim of
then more complicated combined structures are generated,
this paper is to derive the Whitham modulation equations for
which we study in the following sections.
the mKdVB case (5) and apply them to the theory of undular
bores. To this end, we will use the direct Whitham method
[16,18] developed further for the perturbed KdV equation in B. Periodic solutions
Ref. [28]. Its advantage is that it does not require the devel- If the boundary values u± do not satisfy the above con-
opment of quite involved methods of the inverse scattering dition of belonging to the same monotonicity intervals, then
transform (see Ref. [29]). We obtain analytical formulas for the wave breaks and an undular bore forms. In the Gurevich-
the main characteristics of shock waves and confirm them by Pitaevskii approach [15], they are represented by modulated
numerical solutions of Eq. (5). periodic solutions of Eq. (5), so first we have to describe the
We confine ourselves to the case of the defocusing mKdV nonmodulated solutions for zero dissipation.
equation with α > 0, when dispersive shock wave are modula- We look for traveling-wave solutions u = u(ξ ), ξ = x −
tionally stable. The focusing mKdV equation with α < 0 has V t, of Eq. (5) with γ = 0 and after two integrations we get
modulationally unstable periodic solutions (see, e.g., Ref. [30]
and references therein) and their slow modulations are also de- uξ2 = αu4 + Vu2 + 2Bu − 2A, (8)
scribed by the Whitham modulation equations, but this theory where A and B are constants of integration. We assume that
is beyond the scope of the present paper. the polynomial on the right-hand side has four real roots
νi , i = 1, 2, 3, 4, which are ordered according to inequalities
ν1  ν2  ν3  ν4 , so Eq. (8) can be rewritten in the form
II. ELEMENTARY WAVE STRUCTURES
IN mKdVB EQUATION THEORY uξ2 = α(u − ν1 )(u − ν2 )(u − ν3 )(u − ν4 ). (9)
Wave structures evolved from an initial discontinuity are The constants in these two equations are related by the
typically combined from several types of elementary wave expressions
structures and we consider them briefly. For definiteness, we
confine ourselves to the case of a positive coefficient α > 0, V = α(ν1 ν2 + ν1 ν3 + ν1 ν4 + ν2 ν3 + ν2 ν4 + ν3 ν4 ),
although a similar theory can be developed for the case of α
negative α. Naturally, the viscosity coefficient γ is positive. B = − (ν1 ν2 ν3 + ν1 ν2 ν4 + ν1 ν3 ν4 + ν2 ν3 ν4 ),
2
α
A = − ν1 ν2 ν3 ν4 (10)
A. Rarefaction waves
2
and the roots νi are not independent of each other but
First we consider situations when a wave connects two
connected by the formula
trivial solutions u = u− on the left and u = u+ on the right
from the initial discontinuity and assume that during the evo- ν1 + ν2 + ν3 + ν4 = 0. (11)
lution the wave remains a smooth function of x. Then we
can neglect dispersive and dissipative effects proportional to Periodic real solutions can only exist when u oscillates be-
higher-order derivatives of x and describe such a wave in the tween two consecutive roots where the potential curve is
simplest approximation taking into account only nonlinear positive, that is, ν2  u  ν3 , as shown in Fig. 1. Integration
effects proportional to the first-order space derivative of Eq. (9) with the initial condition u = ν3 at ξ = ξ0 gives
 ν3
du
ut − 6αu2 ux = 0. (6) ξ − ξ0 = √ (12)
u α(u − ν1 )(u − ν2 )(u − ν3 )(u − ν4 )
The boundary conditions suggest that there are two character- and standard calculation yields the expression
istic functions, one for the sound wave propagating along the ν3 (ν4 − ν2 ) − ν4 (ν3 − ν2 )sn2 (θ ; m)
plateau u = u− , which has the characteristic xl = −6αu− 2
t, u= , (13)
(ν4 − ν2 ) − (ν3 − ν2 )sn2 (θ ; m)
and the other for the sound wave propagating along the
plateau u+ , so this edge moves according to the equation xr = where sn(θ , m) is the Jacobi elliptic sinus function, with
−6αu+ 2
t. Consequently, the solution consists of three parts: 
u = u− for x < xl , u = u+ for x > xr , and between these two θ = 21 α(ν3 − ν1 )(ν4 − ν2 )ξ (14)

015102-2
UNDULAR BORE THEORY FOR THE MODIFIED … PHYSICAL REVIEW E 109, 015102 (2024)

C. Kink
Here we find the kink solution of Eq. (5) with γ = 0. As
usual, we look for a traveling-wave solution u = u(ξ ), ξ =
x − V t, and assume that u → u− as ξ → −∞. Then trivial
integration taking into account our boundary condition gives
uξ ξ = γ uξ + V (u − u− ) + 2α(u3 − u−
3
). (22)
We also have u → u+ as ξ → +∞, as it should be for a
kink solution. Then we immediately get the expression for the
velocity
V = −2α(u−
2
+ u− u+ + u+
2
); (23)
FIG. 1. Periodic solutions correspond to oscillations in the
interval ν2  u  ν3 where f (u)  0. substitution of this expression into Eq. (22) gives
uξ ξ = γ uξ + 2α(u − u− )(u − u+ )(u + u− + u+ ). (24)
and
Now, following Ref. [31], we assume that this equation has an
(ν4 − ν1 )(ν3 − ν2 ) integral in the form
m= . (15)
(ν4 − ν2 )(ν3 − ν1 )
uξ = a(u − u− )(u − u+ ),
The expression (14) allows us to define the wave number and
the frequency of the periodic wave in terms of parameters νi , that is,
 duξ du
k = α(ν3 − ν1 )(ν4 − ν2 ), ω = kV, (16) uξ ξ = = a2 (2u − u− − u+ )(u − u− )(u − u+ ).
du dξ
where V is given by Eq. (10). Substitution of these expressions into Eq. (24) yields
The cnoidal wave solution (13) reduces to important par-
ticular solutions in special limits. When ν1 → ν2 , and thus a2 (2u − u− − u+ ) = γ a + 2α(u + u− + u+ ).
m → 1 and sn(θ ; m) → tanh θ , we arrive at the bright soliton A comparison of the coefficients before u gives a2 = α or
ν3 − ν1 √
u(ξ ) = ν1 + , (17) a = ± α. (25)
cosh θ − νν43 −ν
2
−ν1
1
sinh2 θ
Then the remaining terms give
propagating along a constant background u = ν1 . When ν3 → γ
u− + u+ = ∓ √ . (26)
ν4 , we obtain the dark-soliton solution 3 α
ν4 − ν2 Finally, elementary integration of the equation
u(ξ ) = ν4 − , (18)
cosh2 θ − νν44 −ν 2
sinh2 θ √
−ν1 uξ = ± α(u − u− )(u − u+ ) (27)
propagating along a constant background u = ν4 . When ν3 → yields
ν2 , we get m → 0, so the elliptical sinus becomes the trigono-
metric one, sn(θ ; 0) = sin θ , and we obtain a harmonic wave u(ξ ) = 21 {u− + u+ ± (u− − u+ ) tanh[α(u− − u+ )ξ )]}. (28)
solution oscillating with very small amplitude around u = ν2 , As one can see, the upper sign corresponds to the growing
√
u(ξ ) = ν2 + 1
(ν − ν2 ) cos 2θ . (19) kink with u+ > u− and u+ + u− = γ /3 α and the lower sign
2 3 corresponds√to the decreasing kink with u+ < u− and u+ +
Finally, if we have simultaneously ν1 → ν2 and ν3 → ν4 , it is u− = −γ /3 α. These are the separatrix solutions joining the
convenient to change the initial condition in such a way that stationary solutions of the second-order equation (24), so there
the integral (12) takes the form are no other kink solutions of this equation.
 ν3
du
ξ= √ (20) III. WHITHAM MODULATION EQUATIONS
u α(u − ν2 )(u − ν4 ) FOR mKdVB THEORY
and elementary integration yields According to Whitham [16,17], the modulation theory
√ can be based on averaging of the conservation laws for
ν2 + ν4 exp[∓ α(ν2 − ν4 )(ξ − ξ0 )]
u= √ . (21) the equation under consideration over fast oscillations in
1 + exp[∓ α(ν2 − ν4 )(ξ − ξ0 )]
the slightly modulated cnoidal wave. The perturbed theory
It is important that due to Eq. (11) the parameters are related of the Whitham modulation method for the mKdVB equation
by the formula ν2 + ν4 = 0 and therefore the left and right can be performed in the same way as it was done for the KdVB
limiting values of u have opposite signs and their absolute equation [28].
values are equal to each other. It is remarkable that an exact Due to condition (11), in this theory there are three inde-
solution of this type exists for the full Eq. (5), taking into pendent parameters that can be chosen arbitrarily from the
account dissipation [31]; we consider this modification of this set νi , i = 1, 2, 3, 4. Therefore, we have to average three
so-called kink solution in the next section. conservation laws. However, it is convenient to replace one

015102-3
DE BRITO AND KAMCHATNOV PHYSICAL REVIEW E 109, 015102 (2024)

of them by the universal law of conservation of the number we can express all averaged quantities in terms of the above
of waves [16,17]. Indeed, a slightly modulated wave can be expressions and arrive at
considered locally as a uniform one with the wave number
and the frequency defined by the expressions (−kWB )t + (−kV WB + B)x = R,
(−kWV )t + (−kV WV + A)x = uR,
k = θx , ω = −θt . (29)
(WA )t − V (WA )x = WAVx . (38)
Consequently, they satisfy the conservation law
These equations can be rewritten in a more convenient form
kt + ωx = 0, (30) D
with the use of the differential operator Dt = ∂t∂ + V ∂x
∂
,
where k plays the role of the density of waves and ω is their 
DWB ∂B
flux. They are still expressed in terms of the local values of = WA − R ,
Dt ∂x
the modulation parameters νi by Eqs. (16). Averaging can 
be performed over a wavelength due to the smallness of DWV ∂A
= WA − uR ,
modulations, Dt ∂x
 ∂V
1 L 1 φ(x, t ) DWA
= WA .
φ = φ dx = √ du, (31) Dt ∂x
(39)
L 0 L f (u)
As we mentioned in the Introduction, the mKdV equa-
where L = k −1 is the wavelength and f (u) = ux2 = α tion is not genuinely nonlinear. Therefore, as in the case of
(u − νi ). Thus, the averaged Eq. (30) can be written as the Gardner equation [22], the relationship between physical
kx + ωt = 0, (32) parameters νi and the most convenient modulation parameters
used in the Whitham equations is not single valued. In a
and it is easy to find two other conservation laws for the similar way, the Whitham equations (39) can be transformed
perturbed mKdV equation to the diagonal form for the Riemann invariants r1 , r2 , and r3
which are related to ν1 , ν2 , ν3 , and ν4 by two different sets of
ut − 6αu2 ux + uxxx = R, (33)
formulas. To obtain one such set of formulas, we choose ν1 ,
so in the averaged form they read ν2 , and ν3 as dependent variables, so ν4 is given by Eq. (11)
and dν4 = −(dν1 + dν2 + dν3 ). Then differentials dV , dA,
ut + −2αu3 + uxx x = R, and dB of the modulation parameters used in Eqs. (39) are
u2 t + − 3αu4 + 2uuxx − ux2 x = 2uR. (34) equal to

In the case of Burgers friction we have R = γ uxx [see Eq. (5)], dV = α[(ν4 − ν1 )dν1 + (ν4 − ν2 )dν2 + (ν4 − ν3 )dν3 ],
but to stress the generality of our derivation we keep it α
dB = − [(ν4 − ν1 )(ν2 + ν3 )dν1 + (ν4 − ν2 )(ν1 + ν3 )dν2
unspecified here. 2
Following Refs. [16–18,28], we express all averaged + (ν4 − ν3 )(ν1 + ν2 )dν3 ],
functions in terms of
α
 dA = − [ν2 ν3 (ν4 − ν1 )dν1 + ν1 ν3 (ν4 − ν2 )dν2
W (A, B, V ) = − uξ du = − f (u)du 2
+ ν1 ν2 (ν4 − ν3 )dν3 ]. (40)

=− αu4 + Vu2 + 2Bu − 2A du (35) Introducing the variables wi = ν4 − νi , we write Eq. (39) in
the form
so that

3
Dνi
du WA,νi = αWA (w1 ν1,x + w2 ν2,x + w3 ν3,x ),
WA = √ = dx = L = k −1 , Dt
f (u) i=1

udu 
3
Dνi α
WB = − √ , WB,νi = − WA [w1 (ν2 + ν3 )ν1,x + w2 (ν1 + ν3 )ν2,x
f (u) Dt 2
i=1
1 u2 du
WV = − √ . (36) + w3 (ν1 + ν2 )ν3,x ] − WA R,
2 f (u)

3
Dνi α
Consequently, we get WV,νi = − WA (ν2 ν3 w1 ν1,x + ν1 ν3 w2 ν2,x
i=1
Dt 2
udu
u = k √ = −kWB ,
f (u) + ν1 ν2 w3 ν3,x ) − WA uR. (41)
k u2 du To diagonalize the last system, we multiply the first, second,
1 2
2
u = √ = −kWV . (37) and third lines by the constant parameters p, q, and r, respec-
2 f (u)
tively, sum the resulting equations, and choose p, q, and r in
df
In view of the relation uxx = 21 du , we have uxx  = 0. After such way that the coefficient of ν1,x on the right-hand side
simple transformations with the use of the mKdV equation, vanishes and the coefficients of ν2,x and ν3,x are equal to each

015102-4
UNDULAR BORE THEORY FOR THE MODIFIED … PHYSICAL REVIEW E 109, 015102 (2024)

other. These conditions determine p, q, and r up to a numerical The Riemann invariants ri are positive and we assume that
factor, and we take the following values: they are ordered according to the inequalities 0 < r1 r2 r3 .
Then the parameters νi are ordered as
p = −(ν2 + ν3 )(ν1 ν4 + ν2 ν3 ),
ν1  ν2  ν3 < 0 < ν4 . (50)
q = −2(ν1 ν4 − ν2 ν3 ),
r = −4(ν2 + ν3 ). (42) The phase velocity V and elliptic modulus m reduce to
r3 − r2
After elementary transformations the resulting right-hand side V = −2α(r1 + r2 + r3 ), m= (51)
r3 − r1
of the sum takes the form
and the wavelength is given by the formula
∂ (ν2 + ν3 )
WA α(ν2 − ν1 )(ν3 − ν1 )(ν4 − ν2 )(ν4 − ν3 ) 2
∂x L=√ K (m), (52)
 α(r3 − r1 )
+ 2(ν1 ν4 − ν2 ν3 )R + 4(ν2 + ν3 )uR . (43) with K (m) the complete elliptic integral of the first kind. The
integral (46) can also be expressed in terms of the Riemann
Calculation of the coefficient before Dν1 /Dt gives invariants,
K1 = pWA,ν1 + qWB,ν1 + rWV,ν1 √ √ √ √ ∂L
I1 = 2( r2 − r1 )( r3 − r1 ) , (53)
ν4 − ν1 (p − qu − ru2 /2)du ∂r1
=−  and similar expressions can be obtained for its counterparts for
2 α(u − ν1 )3 (u − ν2 )(u − ν3 )(u − ν4 )3
 equations derived from Eq. (47) by cyclic permutations of ν1 ,
d (u − ν2 )(u − ν3 ) ν2 , and ν3 . As a result, we arrive at the form of the Whitham
= −(ν4 − ν1 ) = 0. (44) equations for the perturbed mKdV theory,
du α(u − ν1 )(u − ν4 )
√
∂ri ∂ri L r1 r2 r3 R − ri uR
Similar calculation of the coefficient before Dν2 /Dt gives + vi = , (54)
∂t ∂x ∂L/∂ri j=i (ri − r j )
K2 = pWA,ν2 + qWB,ν2 + rWV,ν2
where
= (ν4 − ν2 )(ν4 − ν3 )I1 , (45) 
L ∂ 2αL
vi = 1 − V =V + (55)
where ∂L/∂ri ∂ri ∂L/∂ri

u − ν1 are the standard Whitham velocities for the unperturbed
I1 = . (46)
α(u − ν2 )(u − ν3 )(u − ν4 )3 mKdV equation [18,32].
The definitions (48) and (49) of the Riemann invariants
As one can see, this expression is symmetrical with respect imply that in this case a modulated wave oscillates in the
to the interchange of ν2 and ν3 , so K3 = pWA,ν3 + qWB,ν3 + region ν2  u  ν3 < 0 of its amplitude [see Eq. (50)]. To
rWV,ν3 = K2 . Consequently, we obtain one of the modulation get modulation equations for bores with positive values of
equations in the form the amplitude, it is convenient to take ν2 , ν3 , and ν4 as the
 modulation parameters, so that ν1 = −(ν2 + ν3 + ν4 ), and to
∂ (ν2 + ν3 ) ∂ (ν2 + ν3 )
(ν4 − ν2 )(ν4 − ν3 )I1 +V define the Riemann invariants by the formulas
∂t ∂x
∂ (ν2 + ν3 ) r1 = 41 (ν2 + ν3 )2 , r2 = 41 (ν2 + ν4 )2 , r3 = 41 (ν3 + ν4 )2
= WA α(ν2 − ν1 )(ν3 − ν1 )(ν4 − ν2 )(ν4 − ν3 ) (56)
∂x
 and
+ 2(ν1 ν4 − ν2 ν3 )R + 4(ν2 + ν3 )uR , (47) √ √ √ √ √ √
ν1 = − r1 − r2 − r3 , ν2 = r1 + r2 − r3 ,
√ √ √ √ √ √
and the other two equations can be obtained by cyclic permu- ν3 = r1 − r2 + r3 , ν4 = − r1 + r2 + r3 . (57)
tations of ν1 , ν2 , and ν3 .
For 0 < r1  r2  r3 the parameters νi are ordered according
The terms that do not depend on R have a diagonal form
to
with respect to derivatives, so the three values of any function
of ν1 + ν2 , ν1 + ν3 , and ν2 + ν3 can serve as the Riemann ν1 < 0 < ν2  ν3  ν4 (58)
invariants of the resulting Whitham modulation equations. It
is convenient to define them as and the variable u takes positive values in the interval
0 < ν2  u  ν3 . (59)
r1 = 14 (ν2 + ν3 )2 , r2 = 14 (ν1 + ν3 )2 , r3 = 41 (ν1 + ν2 )2
(48) The Whitham equations (54) for this definition of the Riemann
invariants remain the same. Consequently, one solution of
and the Whitham modulation equations describes two different
√ √ √ √ √ √ modulated wave structures, which is a characteristic feature
ν1 = r1 − r2 − r3 , ν2 = − r1 + r2 − r3 ,
√ √ √ √ √ √ of not-genuinely nonlinear wave equations (other examples of
ν3 = − r1 − r2 + r3 , ν4 = r1 + r2 + r3 . (49) such a behavior can be found in Refs. [22,33,34]).

015102-5
DE BRITO AND KAMCHATNOV PHYSICAL REVIEW E 109, 015102 (2024)

IV. STATIONARY BORES IN mKdVB THEORY a variable the modulus m. The Riemann invariants can be
expressed as functions of m in the following way. The first
As mentioned in the Introduction, after a long enough time
and third of Eqs. (65) give r1 and r2 as functions of r3 :
of evolution, small dissipation stops expansion of undular 
bores and they acquire stationary profiles. The correspond- r1 = 21 [σ1 − r3 − (σ1 − r3 )2 − 4σ3 /r3 ],
ing theory for the KdV-Burgers equation was developed in 
Refs. [25,26,28]. Here we obtain similar solutions for the r2 = 21 [σ1 − r3 + (σ1 − r3 )2 − 4σ3 /r3 ]. (67)
case of the mKdVB theory, following mainly the method of Then, with the use of Eq. (51) for m, we find the formula
Ref. [28]. 
A stationary bore propagates with constant velocity V 3r3 − σ1 − (σ1 − r3 )2 − 4σ3 /r3
m=  , (68)
without a change of the profile determined by the modulation 3r3 − σ1 + (σ1 − r3 )2 − 4σ3 /r3
variables ri = ri (ξ ), ξ = x − V t. Such a stationary profile is
supported by the difference of the values of the wave variable which defines in an implicit form the function r3 = r3 (m)
u at two infinities, so that substitution of this function into Eqs. (67) gives the
 functions r1 = r1 (m) and r2 = r2 (m). Differentiation of m
u as x → −∞ by ξ and substitution of Eq. (64) with R = 0 yield the
u(x, 0) → − (60)
u+ as x → +∞. equation for m:
If there were no dispersion effects, we would get a jumplike dm
viscous shock with velocity determined by the Rankine- = − (m). (69)
dξ
Hugoniot conditions (see, e.g., Ref. [17]). Dispersion effects
transform a jumplike transition between two levels of the u Consequently, we obtain the solution in the implicit form
variable into an oscillatory bore, but the Rankine-Hugoniot  1
dm
conditions are still applicable [27]. Following Whitham’s ξ − ξ0 = , (70)
m (m)
theory of weak shocks [17], we introduce the flux function
Q = −2αu3 so that the dispersionless limit of the mKdV with
equation takes the form the conservation law r1 (r2 − r3 )2 + r2 (r1 − r3 )2 + r3 (r1 − r2 )2
(m) = uR,
ut + Qx = 0 (61) 2α(r1 − r2 )(r1 − r3 )3 (r2 − r3 )
(71)
and then a shock wave propagates with velocity
where uR can also be expressed in terms of the Riemann
Q(u− ) − Q(u+ ) invariants, that is, as a function of m (ξ0 is the position of
V = = −2α(u−
2
+ u− u+ + u+
2
). (62)
u− − u+ the soliton edge of the bore with m = 1 at the initial moment
[It is worth noting that it coincides with the velocity of kinks of time). This completes, in principle, the solution of the
(23) calculated taking into account viscosity, which confirms Whitham equations for a stationary bore. When the func-
the generality of the above argumentation.] This velocity must tion m = m(ξ ) is found, it means that the dependence of the
coincide with the constant velocity V of the bore given by Riemann invariants r1 , r2 , and r3 on ξ is also known. Substitu-
Eq. (51), tion of these functions into the two sets (49) and (57) gives us
two different dependences of the parameters νi , i = 1, 2, 3, 4,
V = −2α(r1 + r2 + r3 ). (63) on ξ . This means that their substitution into the solution (13)
Thus, in stationary solutions, the sum of three Riemann invari- yields two different modulated bores. We distinguish the cor-
ants is constant and Eq. (54) reduces to rect solution by the boundary conditions. Thus, now we are in
√ a position to classify all possible wave structures supported by
dri r1 r2 r3 R − ri uR boundary conditions at infinities in the mKdV theory taking
= , i = 1, 2, 3. (64)
dξ 2α i= j (r j − ri ) into account the small Burgers viscosity.
It is convenient to introduce symmetric functions of the
V. CLASSIFICATION OF WAVE STRUCTURES
Riemann invariants,
FOR JUMPLIKE BOUNDARY CONDITIONS
σ1 = r1 + r2 + r3 , σ2 = r1 r2 + r1 r3 + r2 r3 , σ3 = r1 r2 r3 .
In the region of applicability of the Gurevich-Pitaevskii
(65)
theory based on the Whitham method of slow modulations of
It is not hard to derive equations for them, periodic solutions of the mKdV equation, the general diagram
√ of possible wave structures coincides qualitatively with the
dσ1 dσ2 1 dσ3 σ3
= 0, = uR, = R. (66) diagram obtained in Ref. [22] for the related Gardner equa-
dξ dξ 2α dξ 2α tion without viscosity (see also Ref. [27]). Taking viscosity
Consequently, σ1 is an integral of motion, as it should be. into account leads to two modifications: (i) Undular bores
The theory greatly simplifies if R = 0. In particular, it takes become stationary and (ii) the kinks’ parameters are slightly
place for the Burgers viscosity uxx  = (1/L)(ux )|L0 = 0 due to changed, as shown in Sec. II C. The resulting diagram is
periodicity of u in the main approximation. Then σ3 = const is shown in Fig. 2 and here we derive analytical formulas for the
also an integral of motion and we get an ordinary differential main characteristics of the wave structures and compare them
equation for a sole dependent variable σ2 or any other variable with numerical solutions of the mKdVB equation. All numer-
changing along the bore. It is convenient to choose as such ical simulations in this section are performed with the use

015102-6
UNDULAR BORE THEORY FOR THE MODIFIED … PHYSICAL REVIEW E 109, 015102 (2024)

u− At the soliton edge we have m = 1 and r2 = r1 , that is, ν1 =

√
ν2 = − r3 = u+ (or r3+ = u+ 2
), so
γ
u− = −u+ − √
3 α
u− = u+ σ1 = 2r2+ + u+
2
, σ3 = (r2+ )2 u+
2
. (73)
The values of these two constants of motion must be the
same at both edges of the bore, so simple calculations give
the limiting expressions for the Riemann invariants at the
small-amplitude edge,
3 4
2 5 r1− = u−
2
, r2− = r3− = 21 u+ (u+ + u− ), (74)
and at the soliton edge,
1 6
u+ r1+ = r2+ = 21 u− (u− + u+ ), r3+ = u+
2
. (75)
8 7 Naturally, their substitution into Eq. (63) reproduces the ex-
pression (62) for the velocity of the bore. In addition, we
obtain the necessary expressions for the constants of motion
σ1 = u−
2
+ u− u+ + u+
2
, σ3 = 41 u− u+ (u− + u+ )2 .
2 2
(76)
γ To average the Burgers friction term with uR = γ uuxx , it
u− = −u+ + √
3 α
is convenient to make the replacement u → 2v − s1 , where
√ √ √
s1 = r1 + r2 + r3 . The variable v oscillates in the inter-
FIG. 2. Wave structures supported by the boundary conditions √ √
val r2  v  r3 , so we obtain the expression
u− as x → −∞ and u+ as x → +∞.
 √r3 
16
uuxx  = − Q(v)dv, (77)
L √r2
of the split-step Fourier technique with fourth-order accuracy
√ √ √
similar to the approach in Refs. [35,36]. To avoid interference where Q(v) = α(v − r1 )(v − r2 )(v − r3 )(v − s1 ). The
caused by the periodicity of the boundary conditions in this integral here can be expressed in terms of the Jacobi elliptic
method, we use a very large x domain. integrals, but it is convenient enough for practical calculations
In regions 1 and 5 in Fig. 2, we get just undular bores to keep it in this nonintegrated form.
of different polarities. Let us consider first region 1 where To find the criterion of applicability of our theory, we note
u+ < u− < 0 so that u oscillates in the negative interval ν2  that it is correct as long as the length l of the whole bore is
u  ν3 < 0. Correspondingly, we have to use the formulas much greater than a typical local wavelength L inside it. To
(48) and (49) relating νi and r j . In the small-amplitude limit estimate these two parameters, we turn to the small-amplitude
x → −∞ we have ν2 = ν3 = u− and m → 0, that is, r2 → r3 . limit ξ → −∞ where the Riemann invariants are given by the
Consequently, we get at the left edge of the bore r1− = u−2
and formulas (74). Then Eq. (69) reduces to
− −
r2 = r3 , that is,
dm
= 4γ m, m ∝ exp(4γ ξ ), (78)
dξ
σ1 = 2
u− + 2r2− , σ3 = 2
u− (r2− )2 . (72)
so the bore’s length can be estimated as
l ∼ (4γ )−1 . (79)
Substitution of Eqs. (74) into Eq. (52) gives, according to the
standard definition L = 2π /k, the wavelength
√
π 2π
L=√ =√ . (80)
r3 − r1 (u− − u+ )|u+ + 2u− |
Then the condition L l can be written in the form
32π 2 γ 2
u− − u+ . (81)
|u+ + 2u− |
√
On the axis u− = 0 we get −u+ uc = 4 2π γ , and for
|u− | uc we obtain
uc2 γ2
FIG. 3. Riemann invariants for the bores in regions 1 and 5 u− − u+ ∼ . (82)
3|u+ | |u+ |
and the boundary conditions u− = −0.1 and u+ = −0.5 in region
1 and u− = 0.1 and u+ = 0.5 in region 5. The parameters of the Thus, the applicability region is separated from the line u+ =
equations are equal to α = 0.2 and γ = 0.01. u− by a narrow strip formed by the hyperbola boundary (81).

015102-7
DE BRITO AND KAMCHATNOV PHYSICAL REVIEW E 109, 015102 (2024)

FIG. 4. Numerical evolution of undular bores (blue solid lines) and the stationary asymptotic solution (70) (red dashed lines) at (a) t = 50,
(b) t = 100, (c) t = 2000, and (d) t = 3000. The parameters of the equations are equal to α = 0.2 and γ = 0.01 and the boundary conditions
are u− = −0.1 and u+ = −0.5.

In a similar way, in region 5, where u oscillates in the resulting undular bore wave structures for regions 1 and 5 are
positive interval 0 < ν2  u  ν3 , we have to use the for- shown in Figs. 5(a) and 5(b), respectively.
mulas (56) and (57) relating the Riemann invariants to the As was shown in Ref. [22], we cannot join the boundaries
physical parameters of the wave. We obtain the same for- u− > 0 and u+ < 0 by a single undular bore solution because
mulas (74) and (75) for the limiting values of the Riemann the mKdV equation is not genuinely nonlinear. In this case,
invariants, but to average the viscosity term, we make the the wave structure must contain a kink solution as shown
replacement u = −2v + s1 and obtain again the same formula in Fig. 2 for region 2 and for the symmetrical region 6. In
(77). region 2, we have a decreasing kink joining the right boundary
If we take symmetrical boundary conditions in regions 1 u+ < 0 with the intermediate plateau
and 5 that differ only by signs, then in both cases we get the
same function m = m(ξ ) [see Eq. (70)] and the same plots γ
u∗ = −u+ − √ > u− . (83)
of the Riemann invariants r1 (ξ ), r2 (ξ ), and r3 (ξ ) shown in 3 α
Fig. 3. Their substitution into Eqs. (49) or (57) gives the
dependences νi = νi (ξ ), i = 1, 2, 3, 4, for the modulation pa- This plateau is connected with the left boundary u− < u∗ by
rameters of the bores in regions 1 and 5, correspondingly. the negative undular bore, whose profile can be found in the
These functions νi = νi (ξ ) substituted into Eq. (15) yield the same way as above by means of the replacement u+ → u∗ . In
profiles of bores in these two regions shown in Fig. 4 by red particular, velocities of the kink and the bore are equal to
dashed lines. They are compared with numerical solutions of
the mKdVB equation for different times of evolution and quite
good agreement is found for large values of time, especially Vkink = −2α(u∗2 + u∗ u+ + u+
2
),
for the positions and amplitudes of the leading solitons. In Vbore = −2α(u−
2
+ u− u∗ + u∗2 ). (84)
agreement with the qualitative estimates of Ref. [26], the
leading soliton reaches its stationary state at the characteristic For separation of these two constituents of the whole wave
time of order of magnitude of approximately γ −1 . At the same structure in space, the difference
time, quite a slow convergence to the stationary profile is
observed at the small-amplitude tail of the dispersive shock. 
The numerically found velocity of the shock agrees very well γ
Vkink − Vbore = 2α(u− − u+ ) u− − √
with the analytical formula (62) for the asymptotic state. The 3 α

015102-8
UNDULAR BORE THEORY FOR THE MODIFIED … PHYSICAL REVIEW E 109, 015102 (2024)

FIG. 5. Wave structures calculated numerically (blue solid lines) and analytically (red dashed lines) for α = 0.2, γ = 0.01, and the
evolution time t = 3000. The boundary conditions correspond to different regions in the diagram of Fig. 2: (a) region 1, u− = −0.1 and
u+ = −0.5; (b) region 5, u− = 0.1 and u+ = 0.5; (c) region 2, u− = 0.1 and u+ = −0.8; (d) region 6, u− = −0.1 and u+ = 0.8; (e) region
3, u− = 1.0 and u+ = −0.6; (f) region 7, u− = −1.0 and u+ = 0.6; (g) region 4, u− = 1.0 and u+ = 0.3; and (h) region 8, u− = −1.0 and
u+ = −0.3.

015102-9
DE BRITO AND KAMCHATNOV PHYSICAL REVIEW E 109, 015102 (2024)

must be positive. Hence, to realize such a structure the left Eq. (54) becomes applicable. Cylindrical or spherical mKdV
boundary must satisfy the additional condition equations were derived, for example, in Ref. [37] and they can
γ be written in the form
u− > √ . (85)
3 α d
ut − 6αu2 u + uxxx = − u, (90)
If this condition is not fulfilled, then a combined rarefaction 2(t + t0 )
wave matched with a kink is formed (see the discussion of where d = 1 or 2 for cylindrical or spherical geometry, re-
such situations in Ref. [27]). spectively. For a large enough time of evolution t0 1 the
In region 6 with u− < 0 and u+ > 0 we get a structure with perturbative right-hand side term is small, so the dispersive
a growing kink, so the intermediate plateau has the amplitude shock wave solutions to this equation can be approximated by
γ periodic solutions of the standard mKdV equation with slowly
u∗ = −u+ + √ < u− , (86)
3 α changing parameters, whose evolution is governed by Eq. (54)
and such a structure is realized for with
γ d udu
u− < − √ . (87) R = − k √ ,
3 α 2(t + t0 ) f (u)
We compare analytical and numerical solutions for regions d u2 du
2 and 6 in Figs. 5(c) and 5(d), respectively. Again, quite uR = − k √ . (91)
2(t + t0 ) f (u)
satisfactory agreement is observed.
It is clear that when u− reaches the level u− = u∗ , the The integrals here can be expressed in terms of standard
cnoidal bore disappears and the wave structure reduces to Jacobi elliptic integrals of first, second, and third kinds, so we
a sole kink. After a further increase of u− we get into arrive quite easily at the Whitham equations derived earlier
region 3 where the left boundary u− is joined with the by different methods in Ref. [38] for cylindrical cases and in
plateau u∗ by a rarefaction wave (7). Its left edge propagates Ref. [39] for spherical cases, respectively. Thus, the Whitham
−
with velocity Vrw = −6αu− 2
and its right edge propagates with equations (54) can find various applications besides consider-
+
velocity Vrw = −6αu∗ , which must be smaller than the kink’s
2 ation of the effects of the small viscosity.
velocity. This gives the condition
2γ γ VII. CONCLUSION
u+ < − √ or 0 > u+ > − √ (88)
3 α 6 α The above theory confirms the general statement that
for the realization of such a structure in region 3. A similar weak dissipative effects stabilize the expanding evolution of
structure in the symmetrical region 7 is realized for dispersive shock waves, so after a long enough time, they
converge to stationary structures characterized by some finite
2γ γ
u+ > √ or 0 < u+ < √ . (89) length, which is inversely proportional to the viscosity coef-
3 α 6 α ficient. The appearance of the new parameter leads to some
As one can see in Figs. 5(e) and 5(f), the analytical theory limitations for the applicability of the Whitham method used
agrees very well with the numerical solutions for these two in the Gurevich-Pitaevskii approach to description of bores.
regions. In particular, the condition that the size of the whole shock
Finally, in regions 4 and 8 the boundary values u± have is much greater than the typical wavelength inside the shock
the same signs, so they are connected by standard rarefac- demands that the jump between the boundary conditions is
tion waves with negligible influence of the Burgers friction large enough. Since the mKdV equation is not genuinely non-
[see Figs. 5(g) and 5(h)]. This completes the classification linear, we get combined wave structures consisting of a kink
of possible wave structures supported by different boundary and a cnoidal bore or a rarefaction wave. Small viscosity leads
conditions in the theory of the mKdVB equation. to modification of the kink solution found in Ref. [31] and
the condition that the two structural elements of a combined
VI. WHITHAM EQUATIONS FOR CYLINDRICAL
structure propagate separately from each other also leads to
AND SPHERICAL mKdV EQUATIONS
some limitations for boundary conditions. Although in the
case of small viscosity these restrictions are not essential, one
We obtained the Whitham modulation equations in a quite should keep in mind their existence in the practical application
general form (54) where the expression for the perturba- of the theory.
tion term R in Eq. (33) was not specified. Therefore, this
universal form of the Whitham equations can be applied to
ACKNOWLEDGMENTS
other problems of the dynamics of mKdV dispersive shock
waves. In particular, when we consider cylindrical or spherical This research was funded by the research Project No.
dispersive shock waves whose width is much smaller than FFUU-2021-0003 of the Institute of Spectroscopy of the Rus-
the radius of the whole wave structure, the curvature of the sian Academy of Sciences (Secs. II and III) and by the RSF
shock can be treated as a small parameter of the theory and Grant No. 19-72-30028 (Secs. IV and V).

015102-10
UNDULAR BORE THEORY FOR THE MODIFIED … PHYSICAL REVIEW E 109, 015102 (2024)

[1] R. M. Miura, Korteweg–de Vries equation and generalizations. [22] A. M. Kamchatnov, Y. H. Kuo, T. C. Lin, T. L. Horng, S. C.
I. A remarkable explicit nonlinear transformation, J. Math. Gou, R. Clift, G. A. El, and R. H. Grimshaw, Undular bore
Phys. 9, 1202 (1968). theory for the Gardner equation, Phys. Rev. E 86, 036605
[2] C. S. Gardner, J. M. Greene, M. D. Kruskal, and R. M. Miura, (2012).
Method for solving the Korteweg–de Vries equation, Phys. Rev. [23] T. B. Benjamin and M. J. Lighthill, On cnoidal waves and bores,
Lett. 19, 1095 (1967). Proc. R. Soc. London Ser. A 224, 448 (1954).
[3] R. M. Miura, C. S. Gardner, and M. D. Kruskal, Korteweg–de [24] R. Z. Sagdeev, in Reviews of Plasma Physics, edited by M. A.
Vries equation and generalizations. II. Existence of conser- Leontovich (Consultants Bureau, New York, 1966), Vol. 4,
vation laws and constants of motion, J. Math. Phys. 9, 1204 p. 23.
(1968). [25] A. V. Gurevich and L. P. Pitaevskii, Averaged description of
[4] P. D. Lax, Integrals of nonlinear equations of evolution and waves in the Korteweg–de Vries–Burgers equation, Zh. Eksp.
solitary waves, Commun. Pure Appl. Math. 21, 467 (1968). Teor. Fiz. 93, 871 (1987) [Sov. Phys. JETP 66, 490 (1987)].
[5] M. Wadati, The exact solution of the modified Korteweg–de [26] V. V. Avilov, I. M. Krichever, and S. P. Novikov, Evolution of
Vries equation, J. Phys. Soc. Jpn. 32, 1681 (1972). Whitham’s zone in Kortewe–de Vries theory, Dokl. Akad. Nauk
[6] M. Wadati, The modified Korteweg–de Vries equation, J. Phys. SSSR 295, 345 (1987) [Sov. Phys. Dokl. 32, 564 (1987)].
Soc. Jpn. 34, 1289 (1973). [27] G. A. El, M. A. Hoefer, and M. Shearer, Dispersive and
[7] V. E. Zakharov, S. V. Manakov, S. P. Novikov, and L. P. diffusive-dispersive shock waves for nonconvex conservation
Pitaevskii, The Theory of Solitons: The Inverse Scattering laws, SIAM Rev. 59, 3 (2017).
Method (Nauka, Moscow, 1980). [28] A. M. Kamchatnov, Whitham theory for perturbed Korteweg–
[8] M. J. Ablowitz and H. Segur, Solitons and the Inverse Scattering de Vries equation, Physica D 333, 99 (2016).
Transform (SIAM, Philadelphia, 1981). [29] A. M. Kamchatnov, On Whitham theory for perturbed inte-
[9] A. C. Newell, Solitons in Mathematics and Physics (SIAM, grable equations, Physica D 188, 247 (2004).
Philadelphia, 1985). [30] A. Chowdury, A. Ankiewicz, and N. Akhmediev, Periodic and
[10] R. Grimshaw, E. Pelinovsky, and T. Talipova, The modified rational solutions of modified Korteweg–de Vries equation,
Korteweg–de Vries equation in the theory of large-amplitude Eur. Phys. J. D 70, 104 (2016).
internal waves, Nonlinear Process. Geophys. 4, 237 (1997). [31] D. Jacobs, B. McKinney, and M. Shearer, Traveling wave so-
[11] K. R. Helfrich and W. K. Melville, Long nonlinear internal lutions of the modified Korteweg–de Vries–Burgers equation,
waves, Annu. Rev. Fluid Mech. 38, 395 (2006). J. Differ. Equations 116, 448 (1995).
[12] J. R. Apel, L. A. Ostrovsky, Y. A. Stepanyants, and J. F. Lynch, [32] A. M. Kamchatnov, A. Spire, and V. V. Konotop, On dis-
Internal solitons in the ocean and their effect on underwater sipationless shock waves in a discrete nonlinear Schrödinger
sound, J. Acoust. Soc. Am. 121, 695 (2007). equation, J. Phys. A: Math. Gen. 37, 5547 (2004).
[13] G. A. El and M. A. Hoefer, Dispersive shock waves and modu- [33] S. K. Ivanov and A. M. Kamchatnov, Riemann problem for the
lation theory, Physica D 333, 11 (2016). photon fluid: Self-steepening effects, Phys. Rev. A 96, 053844
[14] A. M. Kamchatnov, Gurevich-Pitaevskii problem and its de- (2017).
velopment, Usp. Fiz. Nauk. 191, 52 (2021) [Phys. Usp. 64, 48 [34] S. K. Ivanov, A. M. Kamchatnov, T. Congy, and N. Pavloff,
(2021)]. Solution of the Riemann problem for polarization waves in
[15] A. V. Gurevich and L. P. Pitaevskii, Nonstationary structure of a two-component Bose-Einstein condensate, Phys. Rev. E 96,
a collisionless shock wave, Zh. Eksp. Teor. Fiz. 65, 590 (1973) 062202 (2017).
[Sov. Phys. JETP 38, 291 (1974)]. [35] G. M. Muslu and H. A. Erbay, A split-step Fourier method
[16] G. B. Whitham, Non-linear dispersive waves, Proc. R. Soc. for the complex modified Korteweg–de Vries equation,
London Ser. A 283, 238 (1965). Comput. Math. Appl. 45, 503 (2003).
[17] G. B. Whitham, Linear and Nonlinear Waves (Wiley Inter- [36] D. Ralf, G. Roland, and O. Hilde, A reliable split-step
science, New York, 1974). Fourier method for the propagation equation of ultra-fast pulses
[18] C. F. Driscoll, and T. M. O’Neil, Modulational instability of in single-mode optical fiber, J. Lightw. Technol. 31, 2008
cnoidal wave solutions of the modified Korteweg–de Vries (2013).
equation, J. Math. Phys. 17, 1196 (1976). [37] A. Mannan and A. A. Mamun, Nonplanar dust-acoustic Gard-
[19] P. D. Lax, Hyperbolic systems of conservation laws II, ner solitons in a four-component dusty plasma, Phys. Rev. E 84,
Commun. Pure Appl. Math. 10, 537 (1957). 026408 (2011).
[20] Y. Kodama, V. U. Pierce, and F.-R. Tian, On the Whitham [38] G. Aslanova, S. Ahmetolan, and A. Demirci, Nonlinear mod-
equations for the defocusing complex modified KdV equation, ulation of periodic waves in the cylindrical Gardner equation,
SIAM J. Math. Anal. 40, 1750 (2009). Phys. Rev. E 102, 052215 (2020).
[21] T. R. Marchant, Undular bores and the initial-boundary [39] N. Ozdemir, A. Demirci, and S. Ahmetolan, Undular bores
value problem for the modified Korteweg–de Vries equation, in the (3+1) dimensional mKP equation, Phys. Lett. A 483,
Wave Motion 45, 540 (2008). 129051 (2023).

015102-11