Modulational Instability of Nonlinear Wave Packets within (2+4) Korteweg–de Vries Equation

Abstract

The higher-order nonlinear Schrödinger equation with combined nonlinearities is derived by an asymptotic reduction from the (2+4) Korteweg–de Vries model for weakly nonlinear wave packets for the context of interfacial waves in a three-layer symmetric media. Focusing properties and modulation instability effects are discussed for the considered physical context. Instability growth rate, maximum of the increment and the boundaries of the instability interval are derived in terms of three-layer density stratification, their structure on the parameter planes of relative layer depth, carrier wavenumber and envelope amplitude, are considered in detail.

1. Introduction

This paper considers the extended modified Korteweg–de Vries equation, called also (2+4) KdV equation (cubic-quintic mKdV or mKdV 3–5), containing a combination of nonlinear terms of the third and fifth degree in the same order of smallness:

$$\frac{\partial \zeta }{\partial t}+c\frac{\partial \zeta }{\partial x}+{\mathsf{\alpha }}_2{\zeta }^2\frac{\partial \zeta }{\partial x}+{\mathsf{\alpha }}_4{\zeta }^4\frac{\partial \zeta }{\partial x}+\mathsf{\beta }\frac{{\partial }^3\zeta }{\partial x^3}=0$$

(1)

where x is a horizontal coordinate, t is time, ζ(x, t) describes the wavefield. It belongs to the class of KdV-type equations, subclass of generalized modified KdV equations, representing dynamics of weakly nonlinear weakly dispersive waves. Equations of this class have been widely studied in recent decades [1,2,3,4,5]. This equation clarifies the character of wave dynamics near the point of zero cubic nonlinearity due to fifth-order nonlinearity and is an analogue of the Gardner equation, which is valid in various branches of physics [6,7,8,9,10], especially for describing situations with a change in the sign of the coefficient of quadratic nonlinearity [11,12]. Equation (1) was, in particular, derived for waves at interfaces in a three-layer fluid flow that is symmetrical relative to the half-depth [13,14]. It also naturally arises in modeling the propagation of longitudinal waves in thin rods, plates and shells when the dependence of the stress intensity on the strain intensity has the form of a fifth-degree polynomial [15,16]. The classification of its exact solitary-wave and periodic solutions is given in [17]. Long linear wave speed, c, coefficients of nonlinearity, α₂, α₄ and the dispersion coefficient, β, are determined by a specific physical situation and depend on what physical field is described by the variable ζ (in a hydrodynamic context, for example, it can serve as displacement of density isolines, velocity potential, one of the components of velocity or vorticity, pressure, etc.) and on the properties of the medium in which the waves propagate. We will assume that the physical quantity ζ represents a perturbation of one of the interfaces between layers in a fluid in the context of internal waves in a layered medium; the assumptions when performing the asymptotic expansion will partly depend on this. From the point of view of wave dynamics, the signs of the coefficients of various terms in Equation (1) are also important.

When studying the equations of KdV hierarchy, the main focus is usually paid on solitonic solutions [18,19,20,21,22,23,24,25]: their properties, stability and the nature of interactions. Questions of the wave group dynamics within the framework of KdV-type equations are raised much less frequently. In the case of weak nonlinearity, the stability of weakly modulated quasi-monochromatic wave packets described by some higher-order KdV equations can be studied using the equivalent modified nonlinear Schrödinger equation for the envelope of the wave field [26,27,28], and determining its type (focusing/defocusing). For the classical KdV equation with quadratic nonlinearity, the corresponding Schrödinger equation is defocusing and the wave packets are stable [29]. In the case of the mKdV equation (with cubic nonlinearity), the corresponding Schrödinger equation is focusing, and wave packets are unstable to modulation, which leads to the formation of large amplitude waves [30,31,32]. For the generalized KdV equation with nonlinear term in a general form ∂f(ζ)/∂x, some considerations are given in [33], for the case of one nonlinear term of arbitrary integer degree, a similar analysis was carried out in [34]. For the Gardner equation, or extended KdV equation with combined quadratic and cubic nonlinearity, which has much richer dynamics, it has been demonstrated [30], that both focusing and defocusing cases are realized depending on the combination of the signs of coefficients in the original equation. In such general cases, it is necessary to obtain the condition for modulation instability in terms of the coefficients of nonlinearities and dispersion of the original KdV-type equation (their signs and combinations of signs are especially important), as well as in terms of the characteristics of the initial wave field.

The aim of this research is to study the modulation instability of wave packets within the framework of the (2+4) KdV equation using the equivalent modified nonlinear Schrödinger equation, which is derived explicitly using an asymptotic expansion in Section 2, and its form is specified with the help of rescaling in Section 3. Furthermore, in Section 4, for the case under consideration, the mechanism of the occurrence of modulation instability is analytically described in detail and explicit expressions are given for the growth rate of instability, the maximum growth rate and the boundaries of the instability interval. In Section 5, these parameters are analyzed in the hydrodynamic context of waves at interfaces in a symmetric three-layer fluid in the Boussinesq approximation. In the Conclusion (Section 6), we list the results and discuss their possible applicability.

2. Derivation of the Higher-Order Nonlinear Schrödinger Equation

A reduction of Equation (1) to a nonlinear Schrödinger equation can be obtained for weakly nonlinear wave packets. We will follow the classical procedure for deriving the Schrödinger equation, see, for example [35]. The assumption of weak modulation of waves (narrowness of the spectrum) allows us to look for a solution to Equation (1) in the form of a superposition of harmonics

$$\zeta (x,t)=\mathsf{\mu }{\displaystyle \sum _{m=-\infty }^{+\infty }{\zeta }_m(x,t)E^m}$$

(2)

where ζ_m(x, t) are the amplitudes of the harmonics

$$E^m=\exp (im(kx-\mathsf{\omega }t))$$

(3)

in this case, the fundamental (carrier) harmonic (m = ±1) of the wavefield has a frequency ω and wavenumber k. In the expansion (2), μ << 1 is a small parameter responsible for the approximation of weak nonlinearity (small, but finite amplitude). For the wavefield ζ(x, t) to be real, it is necessary to satisfy the relation:

$${\zeta }_m(x,t)=\overline{{\zeta }_{-m}}(x,t)$$

where the overbar indicates complex conjugation.

Next, we introduce a slow coordinate x₁, representing weak dispersion (the slow changes of the envelope compared to the carrier wavelength, or the ratio between the envelope wavenumber and the carrier wavenumber), which is characterized by the parameter δ << 1:

$$\frac{\partial }{\partial x}=\frac{\partial }{\partial x_0}+\mathsf{\delta }\frac{\partial }{\partial x_1}$$

(4)

and slow time coordinates—t_j:

$$\frac{\partial }{\partial t}=\frac{\partial }{\partial t_0}+\mathsf{\epsilon }\frac{\partial }{\partial t_1}+{\mathsf{\epsilon }}^2\frac{\partial }{\partial t_2}+\dots$$

(5)

We will look for a solution in the form of expansion of the amplitude ζ_m(x, t) of each harmonic in series (2) as a series by small parameter ε, characterizing the difference in time scales:

$${\zeta }_m(x,t)={\displaystyle \sum _{n=0}^{\infty }{\mathsf{\epsilon }}^n}{\zeta }_{mn}(x_0,x_1,t_0,t_1,\dots )$$

(6)

Small parameters of nonlinearity μ and dispersion δ may relate differently to each other at different ways of scaling, while a small parameter ε, appearing as a marker of terms of different time scales in the asymptotic series (6), must be chosen depending on the relationship between μ and δ. Here, to begin with, we will assume the classical relation of small parameters μ~δ~ε, and start with δ = μ and ε = μ, further using only the parameter μ. This will not prevent one from rescaling the resulting Schrödinger equation if necessary.

Now expression (2), taking into account (6), can be substituted into Equation (1), taking into account (4), (5) and collecting the terms for each harmonic component E^m (3) in turn at fixed value of m and for each order of smallness (at μⁿ). We number the steps of the procedure with a double index {m, n}.

In order {1, 0} we obtain the condition that determines the dispersion relation for the carrier wave:

$$\mathsf{\omega }=ck-\mathsf{\beta }{k}^3$$

(7)

The next order {1, 1} gives the linear equation

$$\frac{\partial {\zeta }_{10}}{\partial t_1}+V\frac{\partial {\zeta }_{10}}{\partial x_1}=0$$

(8)

describing the propagation of a linear wave in a dispersive medium with the group velocity

$$V=c_{gr}=\frac{d\mathsf{\omega }}{dk}=c-3\mathsf{\beta }{k}^2$$

(9)

determined by the dispersion relation (7).

The order {1, 2} leads to the nonlinear evolution equation:

$$i\frac{\partial {\zeta }_{10}}{\partial t_2}-{\mathsf{\alpha }}_{2}k{\left|{\zeta }_{10}\right|}^2{\zeta }_{10}-3\mathsf{\beta }k\frac{{\partial }^2{\zeta }_{10}}{\partial x_1^2}=0$$

(10)

Orders {0, m}, m = 2…5, yield, respectively:

$$\frac{\partial {\zeta }_{0m}}{\partial t_1}+c\frac{\partial {\zeta }_{0m}}{\partial x_1}=0$$

which corresponds to the free solution for the zero harmonic. These terms can be set equal to zero.

The second order evolution equation for the first harmonic is obtained in the order {1, 3} of our expansion:

$$\frac{\partial {\zeta }_{10}}{\partial t_3}+2{\mathsf{\alpha }}_2{\left|{\zeta }_{10}\right|}^2\frac{\partial {\zeta }_{10}}{\partial x_1}+{\mathsf{\alpha }}_2{{\zeta }_{10}}^2\frac{\partial \overline{{\zeta }_{10}}}{\partial x_1}+\mathsf{\beta }\frac{{\partial }^3{\zeta }_{10}}{\partial x_1^3}=0$$

(11)

The expression for the third harmonic amplitude is obtained in the order {3, 2}, and it has the form:

$${\zeta }_{32}=\frac{{\mathsf{\alpha }}_2}{24\mathsf{\beta }{k}^2}{{\zeta }_{10}}^3$$

(12)

The order {1, 4} yields the equation describing the third-order correction for the amplitude of the first harmonic:

$$\frac{\partial {\zeta }_{10}}{\partial t_4}+i\frac{{{\mathsf{\alpha }}_2}^2+48{\mathsf{\alpha }}_4\mathsf{\beta }{k}^2}{24\mathsf{\beta }k}{\left|{\zeta }_{10}\right|}^4{\zeta }_{10}=0$$

(13)

In the order {3, 3} we have a correction to the amplitude of the third harmonic:

$${\zeta }_{33}=i\frac{{\mathsf{\alpha }}_2}{12\mathsf{\beta }{k}^3}{{\zeta }_{10}}^2\frac{\partial {\zeta }_{10}}{\partial x_1}$$

(14)

The remaining corrections in series (6) turn out to be zero during the expansion process:

ζ_m0 = 0 ∀ m ≠ ± 1    ζ_mn = 0 ∀ |m| > n + 1    ζ_1n = 0 ∀ n ≥ 1

The order {1, 5} leads to a fourth order evolution equation for the first harmonic amplitude:

$$\frac{\partial {\zeta }_{10}}{\partial t_5}+\frac{1}{24\mathsf{\beta }{k}^2}{\left|{\zeta }_{10}\right|}^2\left(\left(2{{\mathsf{\alpha }}_2}^2+96{\mathsf{\alpha }}_4\mathsf{\beta }{k}^2\right){\zeta }_{10}\frac{\partial \overline{{\zeta }_{10}}}{\partial x_1}+\left({{\mathsf{\alpha }}_2}^2+144{\mathsf{\alpha }}_4\mathsf{\beta }{k}^2\right)\overline{{\zeta }_{10}}\frac{\partial {\zeta }_{10}}{\partial x_1}\right)=0$$

(15)

Combining Equations (8), (10), (11), (13) and (15), and taking into account (4) and (5), we obtain the following equation:

$$\begin{array}{l}i\left(\frac{\partial {\zeta }_{10}}{\partial t}+V\frac{\partial {\zeta }_{10}}{\partial x}\right)+\mathsf{\mu }\left({\tilde{\mathsf{\alpha }}}_1{\left|{\zeta }_{10}\right|}^2{\zeta }_{10}+{\tilde{\mathsf{\beta }}}_1\frac{{\partial }^2{\zeta }_{10}}{\partial x^2}\right)+\\ +i{\mathsf{\mu }}^2\left({\tilde{\mathsf{\alpha }}}_{21}{\left|{\zeta }_{10}\right|}^2\frac{\partial {\zeta }_{10}}{\partial x}+{\tilde{\mathsf{\alpha }}}_{22}{{\zeta }_{10}}^2\frac{\partial \overline{{\zeta }_{10}}}{\partial x}+{\tilde{\mathsf{\beta }}}_2\frac{{\partial }^3{\zeta }_{10}}{\partial x^3}\right)+\\ +{\mathsf{\mu }}^3{\tilde{\mathsf{\alpha }}}_3{\left|{\zeta }_{10}\right|}^4{\zeta }_{10}+i{\mathsf{\mu }}^4{\left|{\zeta }_{10}\right|}^2\left({\tilde{\mathsf{\alpha }}}_{41}{\zeta }_{10}\frac{\partial \overline{{\zeta }_{10}}}{\partial x}+{\tilde{\mathsf{\alpha }}}_{42}\overline{{\zeta }_{10}}\frac{\partial {\zeta }_{10}}{\partial x}\right)=0\end{array}$$

(16)

where the coefficients are expressed as follows:

$${\tilde{\mathsf{\alpha }}}_1=-{\mathsf{\alpha }}_{2}k, {\tilde{\mathsf{\beta }}}_1=-3\mathsf{\beta }k {\tilde{\mathsf{\alpha }}}_{21}=2{\mathsf{\alpha }}_2 {\tilde{\mathsf{\alpha }}}_{22}={\mathsf{\alpha }}_2, {\tilde{\mathsf{\beta }}}_2=\mathsf{\beta }$$

(17)

$${\tilde{\mathsf{\alpha }}}_3=-\frac{{{\mathsf{\alpha }}_2}^2+48{\mathsf{\alpha }}_4\mathsf{\beta }{k}^2}{24\mathsf{\beta }k} {\tilde{\mathsf{\alpha }}}_{41}=-\frac{2{\tilde{\mathsf{\alpha }}}_3}{k}=\frac{{{\mathsf{\alpha }}_2}^2+48{\mathsf{\alpha }}_4\mathsf{\beta }{k}^2}{12\mathsf{\beta }{k}^2} {\tilde{\mathsf{\alpha }}}_{42}=\frac{{{\mathsf{\alpha }}_2}^2+144{\mathsf{\alpha }}_4\mathsf{\beta }{k}^2}{24\mathsf{\beta }{k}^2}$$

(18)

Equation (16) is a generalization of the classical NLS equation taking into account terms up to the fourth order of smallness. It describes weak modulations of waves within the framework of the dynamics reproduced by Equation (1)—in cases in the equations of the KdV hierarchy, due to the symmetry of the medium, the terms of even nonlinearity and nonlinear dispersion completely degenerate, and the coefficient of cubic nonlinearity in Equation (1) is close to zero and can change sign (transition zone).

The physical field ζ in the same order of accuracy is determined by the formula:

$$\zeta =\mathsf{\mu }\mathrm{Re}\left[{\zeta }_{10}\exp (i(\mathsf{\omega }t-kx))\right]+{\mathsf{\mu }}^3\mathrm{Re}\left[({\zeta }_{32}+\mathsf{\epsilon }{\zeta }_{33})\exp (3i(\mathsf{\omega }t-kx))\right]+O({\mathsf{\mu }}^5)$$

(19)

Modifications of Equation (16) up to the third order of accuracy in amplitude were derived, for example, in [30] for the general case of long waves, and in [36]—in the context of surface gravity waves in fluid. In both cases, the derivation was produced with other than our initial equations.

3. Rescaling

It is easy to see that the assumption of equivalence of the contributions of nonlinearity and dispersion (μ~δ) assumed in Section 2 to the formation of motion patterns is not valid within the context described by Equation (1). In such situations, it is first necessary to establish which from the higher-order terms provides the largest contribution to the motion. As the presence of solitonic solutions in the system in question normally is associated with a specific balance of nonlinear and dispersive terms, it is natural to keep this balance also in the higher-order modifications to the NLS equation. The appearance of Equation (16) suggests that such a balance in symmetric stratified media is possible if μ²~δ. This property actually means that the higher-order nonlinear terms provide a relatively large contribution to the governing equation compared with the dispersive terms. Introducing variables,

$${\zeta }_{10}=\mathsf{\mu }\mathsf{\eta }, \tilde{x}=\mathsf{\delta }x, \tilde{t}=\mathsf{\delta }t$$

Equation (16) can be presented as follows (for simplicity we omit tildes):

$$\begin{array}{l}i\mathsf{\mu }\mathsf{\delta }\left(\frac{\partial \mathsf{\eta }}{\partial t}+V\frac{\partial \mathsf{\eta }}{\partial x}\right)+\left({\mathsf{\mu }}^3{\tilde{\mathsf{\alpha }}}_1{\left|\mathsf{\eta }\right|}^2\mathsf{\eta }+\mathsf{\mu }{\mathsf{\delta }}^2{\tilde{\mathsf{\beta }}}_1\frac{{\partial }^2\mathsf{\eta }}{\partial x^2}\right)+\\ +i\left({\mathsf{\mu }}^3\mathsf{\delta }{\tilde{\mathsf{\alpha }}}_{21}{\left|\mathsf{\eta }\right|}^2\frac{\partial \mathsf{\eta }}{\partial x}+{\mathsf{\mu }}^3\mathsf{\delta }{\tilde{\mathsf{\alpha }}}_{22}{\mathsf{\eta }}^2\frac{\partial \overline{\mathsf{\eta }}}{\partial x}+\mathsf{\mu }{\mathsf{\delta }}^3{\tilde{\mathsf{\beta }}}_2\frac{{\partial }^3\mathsf{\eta }}{\partial x^3}\right)+\\ +{\mathsf{\mu }}^5{\tilde{\mathsf{\alpha }}}_3{\left|\mathsf{\eta }\right|}^4\mathsf{\eta }+i{\mathsf{\mu }}^4\mathsf{\delta }{\left|\mathsf{\eta }\right|}^2\left({\tilde{\mathsf{\alpha }}}_{41}\mathsf{\eta }\frac{\partial \overline{\mathsf{\eta }}}{\partial x}+{\tilde{\mathsf{\alpha }}}_{42}\overline{\mathsf{\eta }}\frac{\partial \mathsf{\eta }}{\partial x}\right)=0\end{array}$$

(20)

Making use of assumption μ²~δ~ε and taking into account that Equation (1) in the context under consideration describes a situation when the coefficient of cubic nonlinearity α₂ in it changes its sign (that is, it is less than or of the order of ε), which means, by virtue of (17), the coefficients are also small ${\tilde{\mathsf{\alpha }}}_1$, ${\tilde{\mathsf{\alpha }}}_{21}$, ${\tilde{\mathsf{\alpha }}}_{22}$, leads to the following equation:

$$\begin{gathered} i\left(\frac{\partial \mathsf{\eta }}{\partial t}+V\frac{\partial \mathsf{\eta }}{\partial x}\right)+{\mathsf{\mu }}^2{\tilde{\mathsf{\alpha }}}_1{\left|\mathsf{\eta }\right|}^2\mathsf{\eta }+{\mathsf{\mu }}^2{\tilde{\mathsf{\beta }}}_1\frac{{\partial }^2\mathsf{\eta }}{\partial x^2}+{\mathsf{\mu }}^2{\tilde{\mathsf{\alpha }}}_3{\left|\mathsf{\eta }\right|}^4\mathsf{\eta } \\ +i{\mathsf{\mu }}^3{\left|\mathsf{\eta }\right|}^2\left({\tilde{\mathsf{\alpha }}}_{41}\mathsf{\eta }\frac{\partial \overline{\mathsf{\eta }}}{\partial x}+{\tilde{\mathsf{\alpha }}}_{42}\overline{\mathsf{\eta }}\frac{\partial \mathsf{\eta }}{\partial x}\right)+i{\mathsf{\mu }}^4\left({\tilde{\mathsf{\alpha }}}_{21}{\left|\mathsf{\eta }\right|}^2\frac{\partial \mathsf{\eta }}{\partial x}+{\tilde{\mathsf{\alpha }}}_{22}{\mathsf{\eta }}^2\frac{\partial \overline{\mathsf{\eta }}}{\partial x}\right)+i{\mathsf{\mu }}^4{\tilde{\mathsf{\beta }}}_2\frac{{\partial }^3\mathsf{\eta }}{\partial x^3}=0 \end{gathered}$$

(21)

Then, taking into account only the first corrections to the wave equation, we obtain the modified equation:

$$i\left(\frac{\partial \mathsf{\eta }}{\partial t}+V\frac{\partial \mathsf{\eta }}{\partial x}\right)+\mathsf{\epsilon }{\tilde{\mathsf{\alpha }}}_1{\left|\mathsf{\eta }\right|}^2\mathsf{\eta }+\mathsf{\epsilon }{\tilde{\mathsf{\beta }}}_1\frac{{\partial }^2\mathsf{\eta }}{\partial x^2}+\mathsf{\epsilon }{\tilde{\mathsf{\alpha }}}_3{\left|\mathsf{\eta }\right|}^4\mathsf{\eta }=0$$

(22)

Equation (22) reproduces the nonlinear dynamics of the envelope of the wave field described by Equation (1), in particular, internal waves in a medium with a vertically symmetrical change in density. Displacements of density isolines or interfaces, if the liquid is layered, are described by Expression (19).

For the convenience of further transformations, we replace the variables x and t by new ones (moving reference frame), where the amplitude of the envelope changes slowly due to the effects of nonlinearity and dispersion:

$$x^{\prime }=x-Vt, t^{\prime }=\mathsf{\epsilon }t$$

In the new variables, Equation (22) is reduced to the extended nonlinear Schrödinger equation (the primes on the independent variables are omitted):

$$i\frac{\partial \mathsf{\eta }}{\partial t}+{\tilde{\mathsf{\alpha }}}_1{\left|\mathsf{\eta }\right|}^2 \mathsf{\eta }+{\tilde{\mathsf{\alpha }}}_3{\left|\mathsf{\eta }\right|}^4 \mathsf{\eta }+{\tilde{\mathsf{\beta }}}_1\frac{{\partial }^2\mathsf{\eta }}{\partial x^2}=0$$

(23)

This equation is also known as the modified NLS equation with the cubic-quintic (CQ) nonlinearity, see, e.g., [37,38] in the context of waves in optical media. In general, Equation (23) is valid for a wide range of physics [39,40,41,42]. Both nonlinear terms in it have the same order of magnitude, and, according to (17), (18) the signs of the nonlinear coefficients ${\tilde{\mathsf{\alpha }}}_1$, ${\tilde{\mathsf{\alpha }}}_3$ in CQ NLS (23) are determined by the values and signs of the nonlinearity coefficients α₂, α₄ and dispersion parameter β (always positive in the context of gravity waves in fluids) in the original Equation (1), that is, any combination of signs of nonlinearity coefficients (${\tilde{\mathsf{\alpha }}}_1$, ${\tilde{\mathsf{\alpha }}}_3$) can occur in CQ NLS (23), and dispersion coefficient ${\tilde{\mathsf{\beta }}}_1$ in (23) always has the opposite sign to the sign of β in Equation (1). Therefore, this equation provides an example of competing self-focusing and self-defocusing nonlinearities [43,44]. The physical context of internal waves in symmetrically stratified three-layer fluid which is of our special interest, is discussed in detail further in Section 5.

4. Modulation Instability

By replacing the unknown function in Equation (23)

$$\mathsf{\eta }=a\exp (i\phi )$$

(24)

one can move on to a set of equations for real functions—amplitude a(x, t) and phase φ(x, t) of the envelope of quasiharmonic wave:

$$\frac{\partial \phi }{\partial t}+{\tilde{\mathsf{\beta }}}_1{\left(\frac{\partial \phi }{\partial x}\right)}^2-\frac{{\tilde{\mathsf{\beta }}}_1}{a}\frac{{\partial }^{2}a}{\partial x^2}={\tilde{\mathsf{\alpha }}}_1{a}^2+{\tilde{\mathsf{\alpha }}}_3{a}^4$$

(25)

$$\frac{\partial a}{\partial t}+{\tilde{\mathsf{\beta }}}_1\left(2\frac{\partial a}{\partial x}\frac{\partial \phi }{\partial x}+a\frac{{\partial }^2\phi }{\partial x^2}\right)=0$$

(26)

Let us introduce the new dependent variable instead of the phase, φ:

$$U(x,t)=2{\tilde{\mathsf{\beta }}}_1\frac{\partial \phi }{\partial x}$$

then system (25), (26) after (1) multiplication (25) by 2${\tilde{\mathsf{\beta }}}_1$ and (2) differentiation with respect to x and multiplication (26) by 2a—reduces to:

$$\frac{\partial U}{\partial t}+U\frac{\partial U}{\partial x}=2{\tilde{\mathsf{\alpha }}}_1{\tilde{\mathsf{\beta }}}_1\frac{\partial a^2}{\partial x}+2{\tilde{\mathsf{\alpha }}}_3{\tilde{\mathsf{\beta }}}_1\frac{\partial a^4}{\partial x}+2{{\tilde{\mathsf{\beta }}}_1}^2\frac{\partial }{\partial x}\left(\frac{1}{a}\frac{{\partial }^{2}a}{\partial x^2}\right)$$

(27)

$$\frac{\partial a^2}{\partial t}+\frac{\partial }{\partial x}\left(a^{2}U\right)=0$$

(28)

From Equations (27) and (28), it follows that the behavior of nonlinear envelope waves is determined by the parameters

$$q_1=2{\tilde{\mathsf{\alpha }}}_1{\tilde{\mathsf{\beta }}}_1\mathrm{and} q_2=2{\tilde{\mathsf{\alpha }}}_3{\tilde{\mathsf{\beta }}}_1$$

whose signs are determined by the signs of the nonlinearity and dispersion coefficients in Equation (23), while the dispersion parameter ${\tilde{\mathsf{\beta }}}_1$ is always negative.

The terms $q_1\frac{\partial a^2}{\partial x}$ and $q_2\frac{\partial a^4}{\partial x}$ in the right hand side of Equation (27) are associated with the nonlinearity of the medium, and the term $2{{\tilde{\mathsf{\beta }}}_1}^2\frac{\partial }{\partial x}\left(\frac{1}{a}\frac{{\partial }^{2}a}{\partial x^2}\right)$ takes into account the “nonlinear” dispersion of the group velocity of the wave packet. These terms can have different orders of smallness. For example, in the dynamic approximation of space-time geometric optics, it is assumed that the modulation is so slow that the group velocity dispersion described by the last term in the right hand side of (27) can be neglected. Then the equations for amplitude and frequency modulation take the form:

$$\frac{\partial U}{\partial t}+U\frac{\partial U}{\partial x}=q_1\frac{\partial a^2}{\partial x}+q_2\frac{\partial a^4}{\partial x}$$

(29)

$$\frac{\partial a^2}{\partial t}+\frac{\partial }{\partial x}\left(a^{2}U\right)=0$$

(30)

Equations (29) and (30) at q₂ = 0 have solutions in the form of Riemann waves, having a local connection between the variables a² and U. Also for this particular case, the effects of group velocity splitting and an increase in the steepness of the leading and trailing fronts of the envelope are known [35]. In general case q₂ ≠ 0, these effects still need to be studied.

Let us consider a monochromatic wave in a physical field with an envelope of the form (24), where a = a₀ = const and U = U₀ = 0, and study its stability with respect to small modulating disturbances. To complete this, we linearize the system of Equations (29) and (30) in the vicinity of the stationary value:

a = a₀ + a₁(x, t),  a₁/a₀ << 1
U = U₁(x, t),  U₁ << 1

Linearized system has the form:

$$\frac{\partial U_1}{\partial t}+2q{a}_0\frac{\partial a_1}{\partial x}=0$$

(31)

$$\frac{\partial a_1}{\partial t}+\frac{1}{2}{a}_0\frac{\partial U_1}{\partial x}=0$$

(32)

where

$$q=-(q_1+2q_2{a}_0^2)=-2{\tilde{\mathsf{\beta }}}_1({\tilde{\mathsf{\alpha }}}_1+2{\tilde{\mathsf{\alpha }}}_3{a}_0^2)$$

(33)

The set of Equations (31) and (32) has two families of characteristics:

$$\frac{dx}{dt}=\pm \sqrt{q{a}_o^2}K$$

which are real at q > 0 (i.e., the system (31), (32) is of hyperbolic type), and they are imaginary at q < 0 (i.e., the system (31), (32) is elliptic).

The solution for (a, U) can be sought in the form:

$$(a_1,U_1)\sim \exp (i(Kx-\mathsf{\Omega }t))$$

with real K. Then one can find from (31), (32), that the frequency Ω and wavenumber K of the envelope satisfy the dispersion relation:

$${\mathsf{\Omega }}_{1,2}=\pm \sqrt{q{a}_0^2}K$$

(34)

This implies that at q > 0, frequencies Ω_1,2 are real, and small initial disturbance does not grow, which means that the wave in this case is stable. The situation is opposite at q < 0. At this case, the frequencies Ω_1,2 are imaginary, and small initial disturbances of the envelope are growing in time like $\exp \left(\sqrt{\left|q\right|}{a}_{0}t\right)$, i.e., the monochromatic wave is unstable, and self-modulation occurs. Thus, the condition for modulation instability in our case has the form:

$$q < 0 \mathrm{or} {\tilde{\mathsf{\beta }}}_1({\tilde{\mathsf{\alpha }}}_1+2{\tilde{\mathsf{\alpha }}}_3{a}_0^2)>0$$

(35)

Returning with the help of (17), (18) to the coefficients of the original Equation (1), let us rewrite this inequality in the following form:

$$12{\mathsf{\alpha }}_2\mathsf{\beta }{k}^2+a_0^2\left({\mathsf{\alpha }}_2^2+48{\mathsf{\alpha }}_4\mathsf{\beta }{k}^2\right)>0$$

(36)

or

$$12\mathsf{\beta }{k}^2\left({\mathsf{\alpha }}_2+4{\mathsf{\alpha }}_4{a}_0^2\right)+{\mathsf{\alpha }}_2^2{a}_0^2>0$$

(37)

or expressing the amplitude a₀:

$$a_0^2>a_{0*}^2=-\frac{{\tilde{\mathsf{\alpha }}}_1}{2{\tilde{\mathsf{\alpha }}}_3}=-\frac{12{\mathsf{\alpha }}_2\mathsf{\beta }{k}^2}{{\mathsf{\alpha }}_2^2+48{\mathsf{\alpha }}_4\mathsf{\beta }{k}^2}$$

(38)

when ${\tilde{\mathsf{\alpha }}}_3$ < 0, or when the denominator is positive (${\mathsf{\alpha }}_2^2+48{\mathsf{\alpha }}_4\mathsf{\beta }{k}^2$ > 0), and

$$a_0^2<a_{0*}^2=-\frac{{\tilde{\mathsf{\alpha }}}_1}{2{\tilde{\mathsf{\alpha }}}_3}=-\frac{12{\mathsf{\alpha }}_2\mathsf{\beta }{k}^2}{{\mathsf{\alpha }}_2^2+48{\mathsf{\alpha }}_4\mathsf{\beta }{k}^2}$$

(39)

when ${\tilde{\mathsf{\alpha }}}_3$ > 0, or when the denominator is negative (${\mathsf{\alpha }}_2^2+48{\mathsf{\alpha }}_4\mathsf{\beta }{k}^2$ < 0). The conditions for the implementation of these inequalities in the context of internal wave dynamics in a three-layer symmetric fluid is discussed in detail in Section 5.

Let us now take into account the effect of dispersive spreading of the wave packet (the last term in the right hand side of (27)). Consider the behavior of small disturbances a₁(x, t) and U₁(x, t) of the monochromatic wave, which are described by a linearized system of equations following from (27), (28):

$$\frac{\partial U_1}{\partial t}+2q{a}_0\frac{\partial a_1}{\partial x}=\frac{2{{\tilde{\mathsf{\beta }}}_1}^2}{a_0}\frac{{\partial }^3{a}_1}{\partial x^3}$$

(40)

$$\frac{\partial a_1}{\partial t}+\frac{1}{2}{a}_0\frac{\partial U_1}{\partial x}=0$$

(41)

where q is determined by Expression (33). Assuming $(a_1,U_1)\sim \exp (i(Kx-\mathsf{\Omega }t))$, we obtain the dispersion relation linking the frequency and wavenumber of the modulation:

$${\mathsf{\Omega }}_{1,2}=\pm a_{0}K\sqrt{q+\frac{{{\tilde{\mathsf{\beta }}}_1}^2}{a_0^2}{K}^2}$$

(42)

It differs from the dispersion relation (34) for the case of nonlinear space-time geometric optics approximation by the presence of an additional term under the radical. Expression (42) looks the same as for the classical NLS equation, but an analogue of the Lighthill parameter q (33) contains additionally the term of squared amplitude a₀ multiplied by the coefficient of higher nonlinearity ${\tilde{\mathsf{\alpha }}}_3$. Thus, taking into consideration the terms of the higher nonlinearity can change the type of Equation (23) for the envelope from focusing to defocusing and vice versa. The type of equation from the point of view of modulation instability and the behavior of the disturbance in this case also depends on the sign of the parameter q. At q > 0, the monochromatic wave is stable, and at q < 0, there exists an interval in the envelope wavenumbers

$$K<K_{*}=\sqrt{\left|q\right|a_0^2/{{\tilde{\mathsf{\beta }}}_1}^2}$$

(43)

where the frequency (42) is imaginary, the monochromatic wave is unstable and the disturbances grow up with the increment:

$$\mathrm{Im}\mathsf{\Omega }=a_{0}K\sqrt{\left|q\right|-{\left(\frac{{\tilde{\mathsf{\beta }}}_1}{a_0}K\right)}^2}$$

(44)

Consequently, taking into account the dispersion of the wave packet envelope does not eliminate the effect of modulation instability, but only limits the range of modulation wavenumbers K, where it takes place. The functional dependence (44) is non-monotonic, it has an extremum at $K_m=K_{*}/\sqrt{2}$. This extremum is maximum, and the maximal value of the increment is given by the following expression:

$${(\mathrm{Im}\mathsf{\Omega })}_{\max }=a_0\frac{\sqrt{\left|q\right|}}{2}{K}_{*}=\frac{\left|q\right|}{2\left|{\tilde{\mathsf{\beta }}}_1\right|}{a}_0^2$$

Substituting q from (33), we obtain:

$${(\mathrm{Im}\mathsf{\Omega })}_{\max }=\left|{\tilde{\mathsf{\alpha }}}_1+2{\tilde{\mathsf{\alpha }}}_3{a}_0^2\right|a_0^2$$

(45)

and one can see that the maximal increment does not depend on the dispersion parameter ${\tilde{\mathsf{\beta }}}_1$. Modulation with the wavelength

$${\lambda }_m=\frac{2\pi }{K_m}=\frac{2\pi }{a_0}\sqrt{\frac{\left|{\tilde{\mathsf{\beta }}}_1\right|}{\left|{\tilde{\mathsf{\alpha }}}_1+2{\tilde{\mathsf{\alpha }}}_3{a}_0^2\right|}}$$

(46)

grows most quickly; it is inversely related to the amplitude. All these quantities can be given in terms of the coefficients of the original Equation (1) using Relations (17) and (18).

5. Context of Interfacial Waves in a Symmetric Three-Layer Fluid at Near-Critical Situation

Many studies have been devoted to the modulation instability of surface water waves to explain the rogue-wave formation [45,46,47,48,49,50], but in the present study we discuss waves in the interior of the fluid, where the nonlinear instabilities were also predicted [51,52,53,54]. Equation (1) was derived, in particular, for long lowest-mode interfacial waves in a three-layer fluid with an almost symmetric configuration [13,14] in which the lower-order nonlinear terms vanish and higher-order contributions govern the behavior of wave phenomena in the system. A simple symmetric situation corresponds to the equal thicknesses h of the uppermost and the lowermost layers provided that the density differences Δρ between the layers are also equal (Figure 1). Density stratification of fluid becomes critical from the point of view of long nonlinear internal waves at a certain relative layer thickness [13,55]:

$$h_{cr}/H=9/26$$

(47)

when the cubic nonlinearity coefficient α₂ in Equation (1) vanishes. The coefficients in Equation (1) are as follows in this context [13]:

$$\begin{gathered} {\mathsf{\alpha }}_2{H}^2/c=-\frac{3(26l-9)}{8l^3} {\mathsf{\alpha }}_4{H}^4/c=-\frac{9(1324l^3-1508l^2+513l-45)}{128l^7} \\ \mathsf{\beta }/(c{H}^2)=-\frac{4l^2-3l}{12} \end{gathered}$$

(48)

where l = h/H. Expanding the above Expressions (48) for the coefficients of Equation (1) into a Taylor series near the point h_cr (δ = (h − h_cr)/H) yields

$${\mathsf{\alpha }}_2{H}^2/c=-\frac{57122}{243}\mathsf{\delta }+O\left({\mathsf{\delta }}^2\right)={\stackrel{⌢}{\mathsf{\alpha }}}_2\mathsf{\delta }+O\left({\mathsf{\delta }}^2\right)$$

(49)

$$\mathsf{\beta }/(c{H}^2)=\left(\frac{63}{1352}+\frac{1}{52}\mathsf{\delta }\right)+O\left({\mathsf{\delta }}^2\right)=\stackrel{⌢}{\mathsf{\beta }}+O\left(\mathsf{\delta }\right)$$

(50)

$${\mathsf{\alpha }}_4{H}^4/c=-\frac{47411260}{59049}+\frac{12084101978}{531441}\mathsf{\delta }+O\left({\mathsf{\delta }}^2\right)={\stackrel{⌢}{\mathsf{\alpha }}}_4+O\left({\mathsf{\delta }}^2\right)$$

(51)

As follows from Expressions (48)–(51), the coefficient α₂ at the cubic nonlinearity is sign-variable in the vicinity of h_cr/H, the coefficient α₄ at the quintic nonlinearity is negative in this region, and the dispersion coefficient β is always positive for internal waves.

Let us write down the criterion for modulation instability of the wave field and its other characteristics for this physical context. Using Relations (33), (17) and (18), we obtain

$$q=-2{\tilde{\mathsf{\beta }}}_1({\tilde{\mathsf{\alpha }}}_1+2{\tilde{\mathsf{\alpha }}}_3{a}_0^2)=-6\mathsf{\beta }k\left({\mathsf{\alpha }}_{2}k+\frac{{{\mathsf{\alpha }}_2}^2+48{\mathsf{\alpha }}_4\mathsf{\beta }{k}^2}{12\mathsf{\beta }k}{a}_0^2\right)$$

(52)

and see, that the sign of growth rate q is determined by the sign of the multiplier in the bracket (coincides with it), and it, in turn, depends on the signs of the nonlinearity coefficients in CQ NLS equation, Equation (23), which are determined through the coefficients of the original Equation (1) through Relations (17) and (18). Let us analyze the signs ${\tilde{\mathsf{\alpha }}}_1$ and ${\tilde{\mathsf{\alpha }}}_3$ depending on the parameters of the environment and the envelope. Cubic nonlinearity coefficient ${\tilde{\mathsf{\alpha }}}_1$ changes its sign at l = l_cr = 9/26. The coefficient of the higher nonlinearity ${\tilde{\mathsf{\alpha }}}_3$ goes to zero and changes sign at

$$\stackrel{\wedge }{k}={\stackrel{\wedge }{k}}_{\pm }=\pm \frac{26l-9}{\sqrt{2\left(5296l^4-10004l^3+6576l^2-1719l+135\right)}}$$

(53)

where $\widehat{k}=kH$ is a nondimensional quantity, characterizing the ratio of fluid total depth to wavelength, i.e., dispersion (the waves are long when $\widehat{k}<1$). Combinations of signs of nonlinearity parameters ${\tilde{\mathsf{\alpha }}}_1$, ${\tilde{\mathsf{\alpha }}}_3$ in Equation (23) in the context of interfacial waves in a three-layer fluid are shown in Figure 2 on the plane $(l,\widehat{k})$. The first sign in a pair corresponds to the sign of ${\tilde{\mathsf{\alpha }}}_1$, the second—to the sign of ${\tilde{\mathsf{\alpha }}}_3$. Combination of signs «− −» provides q < 0 for any values of wave amplitude a₀. Combination «− +» means that the case (39) is realized, and at this situation, q < 0 if the condition $a_0^2<a_{0*}^2$ for amplitude is satisfied. In case when both signs are pluses, q > 0 for any values of wave amplitude a₀, and modulation instability cannot occur in this zone. Finally, the combination of the signs “+ −” gives q < 0 when $a_0^2>a_{0*}^2$. Normalized values of $a_{0*}/H$ on the plane $(l,\widehat{k})$ are given in Figure 2b.

Substituting Expressions (48) for the coefficients into expression for q, (52) we obtain:

$$\stackrel{\wedge }{q}=q\frac{H^2}{c^2}=\left(Q_1(l){\widehat{k}}^2+Q_2(l)\right){\widehat{a}}_0^2+Q_3(l){\widehat{k}}^2$$

(54)

where ${\widehat{a}}_0=a_0/H$ is a nondimensional wave amplitude (normalized by the total fluid depth H)—this quantity characterizes nonlinearity of the wave field,

$$\begin{gathered} Q_1=-9\left(5296l^4-10004l^3+6576l^2-1719l+135\right)/\left(64l^6\right) \\ {Q}_2=-9\left(676l^2-468l+81\right)/\left(128l^6\right)\hspace{1em}{Q}_3=-3\left(104l^2-114l+27\right)/\left(16l^2\right) \end{gathered}$$

Then, expanding (54) into a Taylor series near the critical point l = 9/26 (47), we have:

$$\stackrel{\wedge }{q}=\frac{1963780}{2187}{\left(\widehat{k}{\widehat{a}}_0\right)}^2+\frac{169}{13122}\left(5103-1945684{{\stackrel{\wedge }{a}}_0}^2\right){\widehat{k}}^2\mathsf{\delta }+O({\mathsf{\delta }}^2)$$

(55)

where δ = (h − h_cr)/H. From expansion (55) it is clear, that $\widehat{q}>0$ in the main order; thus, a wave is stable at the critical point and at its nearest vicinity. The behavior of the nondimensional modulation instability parameter $\stackrel{\wedge }{q}(l,\widehat{k},{\widehat{a}}_0)$, determined by the normalized version of (54), is illustrated in Figure 3, depending on the parameter of the undisturbed medium l (which determines the dimensionless thickness of the upper and lower layers in a symmetric three-layer liquid) and the parameter $\widehat{k}$, characterizing dispersion in the range of long waves $\widehat{k}<1$, for different values of nondimensional wave amplitude: ${\widehat{a}}_0=0.05$ (Figure 3a) and ${\widehat{a}}_0=0.1$ (Figure 3b). This is fully consistent with the physical limitations of the medium (the wave amplitude should be less than the layer thicknesses) and the weakly nonlinear approximation (the wave amplitude should be significantly less than the layer thicknesses), i.e., is in the range of applicability of equation (1). The region in l is chosen to be limited taking into account the fact that equation (1) is valid only in the vicinity of the critical point l = h_cr/H = 9/26 (for convenience, this critical value of the layer thickness is marked with a magenta dash-dotted line in Figure 3). The boundary of sign change for parameter $\widehat{q}$ (or the boundary of instability region) is shown in Figure 3 by a white dotted line, instability zone q < 0 lies under this line, it consists of two sub-regions: $\stackrel{\wedge }{k}<{\widehat{k}}_{-}^{*}(l)$ when l < 9/26 and $\stackrel{\wedge }{k}<{\widehat{k}}_{+}^{*}(l)$ when l > 9/26. Curve ${\widehat{k}}_{-}^{*}(l)$ does not fall below than ${\widehat{k}}_{-}(l)$ (53) when a₀ increases, and curve ${\widehat{k}}_{+}^{*}(l)$, in opposite, rises no higher than ${\widehat{k}}_{+}(l)$ (53) at the same time. That is, for h < h_cr for more intense waves the threshold of modulation instability shifts to the region of longer waves, and the amplitude here is limited from above by the value a_0*, and for h > h_cr, on the contrary, the threshold of modulation instability shifts to the region of shorter waves, and the critical amplitude value a_0* must be exceeded for instability to occur.

Thus, the phenomenon of modulation instability can be realized under certain combinations of environmental conditions, wavelengths and amplitudes of waves at the interfaces in a symmetrical three-layer liquid.

The upper limit of the interval of dimensionless wave numbers of the envelope $K_{*}H$ (43) as a function of l (the dimensionless thickness of the upper and lower layers in a symmetrical three-layer fluid) and $\widehat{k}$ (dimensionless wave number of the carrier wave) is shown in Figure 4 for ${\widehat{a}}_0=0.05$ (Figure 4a) and ${\widehat{a}}_0=0.1$ (Figure 4b) in the zone of the modulation instability q < 0. Near the boundary of the instability zone q = 0, values of $K_{*}$ are small, and they increase with distance from the boundary. At smaller values of the dimensionless amplitude ${\widehat{a}}_0$, values $K_{*}H$ are also smaller ceteris paribus. Values of $K_{*}H$ also grow when carrier wave number $\widehat{k}$ decreases. It is also necessary to take into account the implied ratio of carrier and envelope wavelengths k/K >> 1, with which constraint (43) must be compatible.

Maximal increment ${(\mathrm{Im}\mathsf{\Omega })}_{\max }H/c$ (45) as a function of l (the dimensionless thickness of the upper and lower layers in a symmetrical three-layer fluid) and $\widehat{k}$ (dimensionless wave number of the carrier wave) is shown in Figure 5 for ${\widehat{a}}_0=0.05$ (Figure 5a) and ${\widehat{a}}_0=0.1$ (Figure 5b) in the zone of the modulation instability q < 0. The increment grows with increasing amplitude, as well as with increasing carrier wavelength.

6. Conclusions

The conditions for the occurrence of modulation instability of weakly nonlinear quasi-monochromatic wave fields are described within the framework of the extended modified KdV equation (or (2+4)-KdV equation), which previously arose in the context of the hydrodynamics of a stratified fluid. For this purpose, an extended nonlinear Schrödinger equation of higher order associated with (2+4)-KdV was derived, which, taking into account the specifics of the considered physical context of internal waves in a three-layer liquid with symmetric stratification, where the sign of cubic nonlinearity changes, is reduced to cubic-quintic NLS equation at such a near-critical situation. Thus, it is possible to construct an approximate solution of the (2+4)-KdV equation (wave field envelope) by solving the equivalent CQ NLS equation. Furthermore, within the framework of the obtained CQ NLS equation, the instability growth rate, maximum of the increment and the boundaries of the instability interval, as well as the corresponding threshold values of the envelope amplitude are obtained analytically. It is shown that the phenomenon of modulation instability can be realized under certain combinations of environmental conditions, wavelengths and amplitudes of waves at the interfaces in a symmetrical three-layer liquid; these relationships are given analytically and examined in detail on the problem parameter planes. These results can be useful in planning numerical and laboratory experiments, as well as in explaining observed field data (within the limits of applicability of the considered approximate models).