Formation of the dynamic energy cascades in quartic and quintic generalized KdV equations

: In this study we investigate the formation of dynamical energy cascades in higher order KdV-type equations. In the beginning we recall what is known about the dynamic cascades for the classical KdV (quadratic) and mKdV (cubic) equations. Then, we investigate further the mKdV case by considering a richer set of initial perturbations in order to check the validity and persistence of various facts previously established for the narrow-banded perturbations. Afterwards we focus on higher order nonlinearities (quartic and quintic) which are found to be quite different in many respects from the mKdV equation. Throughout this study we consider both the direct and double energy cascades. It was found that the dynamic cascade is always formed, but its formation is not necessarily accompanied by the nonlinear stage of the modulational instability. Direct cascade structure remains invariant regardless the size of the spectral domain. In contrast, the double cascade shape can depend on the size of the spectral domain, even if the total number of cascading modes remains invariant. Results obtained in this study can be potentially applied to plasmas, free surface and internal wave hydrodynamics.

(here ℓ is characteristic size of a vortex). A similar approach applied to weakly nonlinear wave systems resulted in the development of the so-called Weak (or wave) Turbulence Theory (WTT) [3] where the spectrum is not universal anymore and it depends on the dispersion function. However, the spectrum shape still remains a powel law as in the strong turbulence theory. In terms of the nonlinearity parameter ε the strong turbulence theory corresponds to ε ST ∼ O(1) while the WTT is situated in the range ε W ∼ O(10 −2 ).
The model of the dynamic energy cascade for the case of intermediate (moderate) nonlinearities ε M ∼ O(10 −1 ) was proposed quite recently in 2012 [4] and studied analytically in [5,6]. Since the dynamic energy cascade model is based on the Modulational Instability (MI) mechanism, the first developments were made the the case of the Nonlinear SCHRÖDINGER (NLS) equation and its modifications. The spectrum shape is deduced from the Incremental Chain Equations (ICE) which strongly relies on the knowledge of the instability interval. Fortunately this information is readily available for NLS-type equations in the literature.
On the other hand, the MI can be encountered in other types of model equations. It is important to study the dynamic energy cascade formation in these models as well. The present article is entirely devoted to some generalized KORTEWEG-DE VRIES (KdV)-type equations. The celebrated original KdV equation was derived for the first time for surface waves in shallow waters [7,8] in order to explain the phenomenon of solitary waves: However, it is known that the classical KdV equation 1 does not possess the MI. Later some extended versions of the KdV equation were derived for both surface [9] and internal [10] waves. By making appropriate changes of variables, it can be reduced to the following modified (m)KdV equation: It is known that equation (2) is integrable and it possesses the MI for β > 0. For negative values of parameter β < 0 its behaviour is similar to the KdV in the respect of the MI. Formally all these model equations can be encompassed into a single gKdV framework: gKdV(u) . = u t + u xxx + (αu + βu 2 + δu 3 + γu 4 )u x = 0, which generates the KdV, mKdV, quartic and quintic gKdV equations for particular choices of coefficients α, β, δ and γ. Quartic and quintic gKdV equations will be studied extensively numerically below.
In general, there are only few known mathematical facts about higher order KdV-type models. The global well-posedness of the quartic gKdV equation in some critical spaces was shown in [11], while the stability of periodic travelling waves in the gKdV equation was studied in [12,13].
As we already mentioned above, the shape of the energy spectrum was computed in the NLS-type equations based on the explicit knowledge of the instability interval with finite measure. In the case of the KdV-type equations and nonlinearity of order nonlinearities O(10 −1 ) such finite intervals do not exist, since any mode with the wave vector k 0 > 0 is modulationally unstable (if the equation has the MI, of course). This fact prevents the direct transposition of NLS predictions to the family of KdV-type equations. Nevertheless, we decided to study the dynamic energy cascades in several KdV-type models numerically.
In general, the MI is an effect which describes a special type of the instability of a single spectral component subject to narrow band excitation in an experiment (regardless whether it is numerical or laboratory). In our previous investigations [14,15] we followed this ideology and all the numerical experiments toward the observation of the direct and inverse cascades had only two spectral harmonics in the initial condition (i.e. the base wave k 0 along with its modulation K 0 ≪ k 0 ).

of 21
On the other hand, in real-world laboratory experiments it is extremely difficult to excite precisely only two spectral harmonics. It is always a finite range which appears. Consequently, it is interesting to study this situation in numerical experiments as well in order to determine the limits of the survivability of the dynamical energy cascade. Additional motivation comes from the fact that the so-called BENJAMIN-FEIR Index (BFI) [16], which characterizes the likelihood of freak wave formation in an open ocean, is also based on the effect of the MI in the presence of a fairly broad spectrum. There is a direct connection between the BFI and the increment chain equations which describe the cascade formation [17]. It should be noted that the BFI was introduced, strictly speaking, for surface waves in deep water which are described by the NLS-type equations. However, this notion might be transposed to the mKdV case as well, since the analytical bridges between NLS and (m)KdV are well studied [18] (and the references therein).
Even if historically the KdV equation (and related models) was proposed first as a model equation for surface waves, later it found applications in many other fields such as plasma physics, acoustic and internal waves in the ocean. That is why the study of this equation along with its various modifications is important. The big advantage of the classical KdV (1) and mKdV (2) equations consists in their exceptional mathematical propertiesi.e. the integrability, which allows a deep mathematical study. However, even in this case it does not allow to describe exactly the properties of the dynamical energy cascade. If we had a closed form analytical solution u(x, t) than we could develop it in a FOURIER series Then, the FOURIER spectrum is given by the graph of the application k → |c k (t)| 2 . Unfortunately, even for integrable equations we do not always have closed form solutions (with notable excpetion of solitons and cnoidal waves). The modern notion of integrability is related to the construction of the LAX's pair or of conserved quantities (in the number equal to the number of the degrees of freedom, i.e. infinity for PDEs). There is no general method to find a LAX's pair. Thus, the integrability does not always yield the solvability of an equation in practice. Several examples can be found in [19]. Consequently, the priveleged method of studying the non-integrable and even integrable, but not solvable equations is the numerical simulation. The formation of the dynamical energy cascade in the integrable mKdV equation (2) was studied numerically in our previous papers [14,15]. However, it was done for narrow initial perturbations in the FOURIER space, since originally it was discovered for the NLS equation which is valid only in this narrow band approximation [4]. On the other hand, the KdV equation does not have limitations on the spectral band, which allows us to study this equation more extensively with broader perturbations.
The present article is organized as follows. In Section 2 the mKdV equation is studied subject to broader initial perturbations. Quintic and quartic higher order modifications of the KdV equation are considered in Section 3. Finally, the main conclusions of the present study are outlined in Section 4.

mKdV equation with broad spectrum initial perturbation
In order to investigate the broader perturbations we considered the following set of numerical experiments which generalize those presented in [14,15]. Here again, in order to solve numerically the mKdV equation on a periodic domain we use the classical FOURIER-type pseudo-spectral method [20,21]. All the derivatives are computed in the FOURIER space, while all the nonlinear products are performed in the physical space (with the linear CPU-time, i.e. O(N), where N is the number of collocation points necessarily equal to the number of FOURIER modes). Thanks to the FFT algorithm [22][23][24] the passage between these two representations is done in the super-linear time O(N log(N)), which determines the overall algorithm complexity (per time step) O(mN log(N)), m ∼ 1 being the number of FFTs needed to evaluate the right-hand side. For the dealiasing we used the classical 2/3-rule [21] which was combined (when necessary) with the FOURIER smoothing method proposed in [25] for more delicate treatment of higher frequencies. The discretization in time was done with the embedded adaptive 5 th order CASH-KARP RUNGE-KUTTA scheme [26] with the adaptive PI step size control [27]. We note that for the numerical simulation of the gKdV equation one needs to set more stringent tolerances, because of the presence of highly nonlinear products. The accuracy of dynamical computations was checked by following, for example, the L 2 norm of the solution u(x, t) which is one of the invariants of the (g,m)KdV dynamics. In all the experiments this norm was conserved up to the 8 th digit with smaller oscillations around the mean value.

Direct cascade
The numerical set-up used in this study is a direct generalization of the test case used also earlier in [28]. We consider the following initial condition posed on a periodic domain [−ℓ, ℓ] = −π/k 0 , π/k 0 : where a is the base wave amplitude, δ is the perturbation magnitude and the wavenumbers k 0 , K 0 and j are chosen such that their ratio k 0 /(jK 0 ) ∈ Z. The last constraint comes from the periodicity conditions on the initial data. Moreover, the spectral gap between the base wave k 0 and its modulation K 0 contains two orders of magnitude (see Table 1 for the values of parameters). When we perturb the base wave according to (4) with a few modes we have to keep in mind that the pertubation has to be different from the base wave at least with one order of magnitude. Taking into account this constraint along with k 0 /(jK 0 ) ∈ Z, we obtain that for our numerical set-up j is roughly limited by 10. The case j ≡ 1 corresponds to the computations already discussed in [14,15,28]. The values of parameters are given in the Table 1. On Figure 1 we show the corresponding pertubed initial conditions. The case Figure 1(a) corresponds to j = 1. The FOURIER spectrum on the bottom panel shows three points which correspond to the base wave k 0 (the most energetic mode) and two side bands k 0 ± K 0 . On Figure 1(b) we depict the perturbation by k 0 ± K 0 , k 0 ± 2K 0 , k 0 ± 3K 0 and k 0 ± 4K 0 modes (in total nine points on the bottom panel). Finally, the last Figure 1(c) corresponds to j = 1, 2, 3, 4, 6, 8 and 10 (counting 15 modes in total). The evolution of three initial conditions from Figure 1 is simulated using the mKdV equation (2). The numerical snapshots taken at the same moment of time t = 160.0 are shown on Figure 2. It can be readily seen that the nonlinear stage of the MI develops much faster if the initial perturbation contains more spectral modes. For instance, the first initial condition with j = 1 unstable mode shows some signs of the MI. Namely, the cascading modes are already present in the FOURIER domain. However, their broadening is not developed yet. It is translated to the physical space by the absence of any amplification of the base wave, while on Figures 2(b,c) the MI is fully developed with amplication factors up to two times.
In order to show that the MI can be fully developed even from a single perturbation mode we present a numerical snapshot on Figure 3 corresponding to the initial condition depicted at Figure 1     but taken at the moment of time t = 620.0. A comparable amplification is achieved. However, it happens almost four times slower than for initial conditions shown in Figures 1(b,c).
This series of numerical experiments shows that a broader initial perturbation not only does not suppress the dynamical energy cascade, but also the position of all the cascading modes is not affected. The cascade shape remains exponential and it was studied in details (in the case j = 1) in [14,15]. We performed the same battery of tests for numerous pertubative modes (see Figure 1(b,c)) and with double initial amplitude (thus, double initial wave steepness). In this study we show only one snapshot corresponding to seven spectral modes and the initial amplitude a = 0.16 = 2 × 0.08. Despite the violent MI observed in this picture ( Figure 4, upper panel) the structure of the cascade still can be recognized, even if some deviation from the exponential shape can be recognized.

Double cascade
As it was demonstrated in [15] for the mKdV equation (2), a shift in the base mode in FOURIER space towards higher frequencies can lead to the formation of a double energy cascade. As in the case of the direct cascade, the experiments were performed with only one unstable mode. In this study we enrich the perturbation with up to seven spectral modes in order to see whether this double cascade can still be observed in such extreme conditions. The numerical parameters for this series of experiments are provided in Table 2. In the previous work [15] it was notices that the increase in the spectral domain yields the appearance of a second occurrence of the double energy cascade. Consequently, we studied several sizes of the computational spectral domain. Here we report the results for N = 2 × 16 384 in order to demonstrate the formation of the double cascade along with its repetition in two occurrences if the spectral domain can host them.
The corresponding results are shown on Figure 5 for the initial time t = 0 (panel (a)), t = 5.0 (panel (b)) and the same time t = 160.0 (panel (c)) as in experiments reported above.
As a conclusion, we observe that the double cascade appears almost immediately (see Figure 5(b, lower panel) corresponding to t = 5.0) and it is accompanied by a fully developed MI in the physical space (see Figure 5   results [15] shows that even the location of cascading modes appears to be the same up to the graphical resolution.
Our main goals are: • Investigate the existence of the direct and inverse energy cascades • If they exist: to compare their formation time scales to compare their structure, i.e. the number of cascading modes, their location, spacing, etc.
In order to meet these goals we will take the same numerical set-up as described above for the direct and double cascades correspondingly (see Tables 1 & 2). The initial perturbations contain one, four or seven modes as it is shown on Figure 1. Several dozens of computations have been performed for each of two equations (5), (6). Below we will show the main findings.

Direct cascade
The selected results of our numerical simulations are shown on Figure 6 for the initial condition consisting of seven perturbative modes and the base wave amplitude a = 0.08. This choice of parameters comes from the fact that we would like to highlight the differences with the mKdV equation (2). The upper pane of Figure 6(a) shows the time moment t = 160.0, while Figure 6(b) reports the system state at t = 3600.0. We can see that during the system evolution new cascading modes appear in the FOURIER spectrum tail (high frequencies). The structure and location of peaks remain the same, while their magnitudes can slightly grow with time.
However, there exist an important qualitative difference with the mKdV case reported in the previous Section. Namely, even at t = 160.0 (with all other equal parameters and the initial perturbation) the physical space of the mKdV equation already demonstrates the nonlinear stage of the MI development (see Figures 2(b,c)). On the other hand, the quartic gKdV even for t = 3600.0 (more than twenty times longer!) does not show any signs of the amplification in the physical space (see Figure 6).
In order to show that the nonlinear MI stage does not occur in the quartic gKdV equation, we performed an additional series of the experiments with higher base wave amplitude. On Figure 9(a) we show the physical and FOURIER spaces for the case of double initial amplitude a = 2 × 0.08 = 0.16, while on the lower panel Figure 9(b) we show simultaneously FOURIER spectra for the simple and double amplitudes at t = 3600.0. As in previous cases, one can see that there is no significant solution amplification in the physical space. Concerning the FOURIER space, the superposition shown on Figure 9(b) clearly shows the coincidence of the cascading modes locations, while new high frequency cascading modes appear. We note that the FOURIER spectrum shown with red dots contains roughly speaking four times more energy (since the initial amplitude is doubled). However, this extra energy is mainly redistributed towards the appearence of new peaks in the tail of the FOURIER spectrum, by conserving the energy of several first cascading modes at the long waves end of the spectrum. Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 9 July 2020 doi:10.20944/preprints202007.0171.v1

Double cascade
Numerical experiments to study the double energy cascade formation in the quartic gKdV equation were performed similarly to the study of the mKdV equation reported in the previous Section. The main conclusions are the same as for the direct cascade. Namely, in the physical space we do not witness any significant amplification of the base wave.
On the other hand, in the FOURIER space several entities of the double cascade are formed. The very interesting fact is that the choice of the spectral domain influences the shape of the double cascade for the quartic gKdV equation. In each spectral domain the total number of cascading modes remains invariant (in the particular case under consideration this number is equal to 19 modes). However, the energy distribution among these modes is different and actually depends on the size of the spectral domain. As we observed in [15], the spectral domain consisting of N = 2 × 16 384 FOURIER modes can host one entity of the double cascade, while doubling the size of the spectral domain yields logically the appearance of two entities of the cascade (with the total number of cascading modes being 14 in that computation). In order to illustrate this fact for the quartic gKdV equation as well, on Figure 8 we show the formed double cascade for different sizes of the spectral domain.
All pictures are taken at t = 3600.0, even though the location of the modes and the distribution of energy among them is already clearly visible at t = 60.0. Longer simulations are necessary in order to show that we do not observe any nonlinear development of the MI in the physical space on the longer horizon. The effect of the energy redistribution in the spectral domain deserves a separate study. It will be performed elsewhere.

Quintic gKdV equation
The numerical set-up in this Section is precisely the same as above. However, here we focus on the quintic gKdV equation (6) with parameter µ = 6.0.
The time evolution under the quintic gKdV equation dynamics (6) is totally different from what we observed for the quartic gKdV (5). Namely, in the quintic case a well-developed MI can be observed both in FOURIER and physical spaces with physical amplifications of the base wave amplitude going up to two times and higher.
On the other hand, these numerical results are qualitatively the same as for the mKdV equation (2) reported above and in [14,15]. Namely, the direct and double cascades are accompanied by the nonlinear stage of the MI in the physical space. It can be qualitatively observed on Figure 9(b). An interesting difference with the mKdV consists in the appearance of the instability treshold in the quintic gKdV equation (6). Below this supposed instability treshold the quintic gKdV behaves similarly to the quartic gKdV equation in the physical space. It is manifested by the long time dynamics of the modulated solution in the physical space without the development of the nonlinear stage of the MI (see Figure 9(a)). In contrast, when we take the doubled initial base wave amplitude a = 2 × 0.08 = 0.16, the MI in the physical space is well developed already at t = 1200.0 as it can be seen on Figure 9(b).
Concerning the FOURIER space, the number of cascading modes, the spacing between them and the shape of the energy spectrum is preserved regardless the base wave amplitude (below or above the instability treshold, see Figure 9). In contrast to the double cascade, the increase in the spectral domain does not lead to the redistribution of the energy among cascading modes. We do not provide the corresponding pictures here. They can be obtained from the authors after a simple request.
For the double cascade in the quintic gKdV we observe the same entities as in the mKdV equation (see Figure 5 and also [15]). However, what is new is the presence of the instability threshold as it was explained for the direct cascade. In order to see the difference the reader can compare

Discussion
In the present study we investigated several KdV-type equations involving higher order nonlinearities. The main focus was on the formation and properties of the dynamical energy cascade in the FOURIER space. The second goal was to describe the cases when the formation of this cascade is accompanied by the nonlinear stage of the MI development in the physical space. This question is of capital importance since the nonlinear stage implies important wave amplifications which are put forward nowadays as one of the main physical mechanisms of the freak wave formation [29]. Since, the total number of numerical experiments performed in this study exceeded a few hundreds, we included only the most typical ones.
• First of all, we extended the study of the mKdV equation presented in [14,15] for the case of a broader initial perturbation of the base wave. It was found that previous findings can be transposed to this case as well: formation of the direct and double cascades, their structure is preserved and appears to be quite stable. However, in the broad perturbation case all nonlinear processes are accelerated (i.e. they develop in shorter times). • Not only the general structure of the energy cascade is preserved, but also the number of the cascading modes (for the fixed spectral domain and the base wave amplitude) and their location in the FOURIER space. • The main difference between the direct and double cascades consists in the fact that the change of the size of the spectral domain does not induce the energy redistribution in the direct cascade, while it changes completely the number of entities, and thus the energy distribution, in the double cascade(s). • In all the mKdV cases, the formation of dynamical cascades is accompanied by the nonlinear stage of the MI development and the substantial amplification of the base wave in the physical space. • In the case of the quartic gKdV equation, all the results about the dynamical cascades formation in the spectral domain remain the same. However, no amplification of the base wave amplitude was observed in the physical space. • The quintic gKdV case is mostly similar to the mKdV equation (with positive nonlinearity) except for the presence of the amplitude threshold in the physical space, below which the amplification does not occur on the time horizons investigated in this study. However, if the amplitude is taken above this threshold, the base wave amplification occurs on longer time scales (virtually three times slower than in mKdV for the parameters considered in this study).
The main findings are summarized in a compact form in Table 3.
Strictly speaking the theory of the dynamical energy cascade [5,14] as well as the model equations considered above are valid at best only for moderate nonlinearities ε ∼ O(10 −1 ). However, Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 9 July 2020 doi:10.20944/preprints202007.0171.v1 one can ask a legitimate question: what happens if we push the nonlinearity in the wave system to the finite size ∼ O(1) as in the strong turbulence? For this purpose we performed an additional experiment for the quartic KdV equation (5) with the base wave amplitude a = 10 × 0.08 = 0.8. Together with seven perturbative modes the initial condition reaches ε ∼ O(1) announced above. Two numerical snapshots are shown on Figure 11. First of all, one can notice that there is no visible difference between t = 20.0 depicted on the upper panel and t = 100.0 shown below. It means that the fully developed quasi-stationary stage is reached on times T ∼ O(1). However, most surprisingly even for such a strong nonlinearity one can still recognise the presence of the direct dynamic energy cascade, despite the fact that we step out from all reasonnable ranges of model and theory applicability.
Recently obtained theoretically results, presented in [30,31], indicate some direction for further studies in order to be able to interpret in detail the available results of numerical simulations and get answers to at least some of the open questions. So, the modulation instability for the KORTEWEG-DE VRIES (KdV)-family of equations u t + s u p u x + u xxx = 0 , with s = ±1 and p > 0 being an arbitrary integer, has been studied in [30]. As modulation instability is defined for the NLS-type of equation, the main question was to deduce an envelope equation for the KdV-family (7). It was proven that general form of envelope equations is different for p = 2 n (symmetric case) and for p = 2 n − 1 (asymmetric case). Moreover, in the latter case the coefficients of the envelope equations depend explicitly on a numerical factor (called pedestal) which in turn is defined by the parameters of initial excitation. This is in accordance with our results demonstrating substantial difference between quartic and quintic cases. The coefficients of the envelope equations have been deduced explicitly in [31] and key characteristics of the modulation instability (that is, the instability increment, its maximum, boundary and maximal wave numbers) were written out. Again, the clear difference between symmetric and asymmetric cases has been demonstrated. In particular, the increment of instability in the latter case depends on the initial conditions in a way which assumes the existence of a lower threshold needed for the developing of instability (see [31,Equation (57)]). This is in agreement with our conclusions.
Summarizing, these novel theoretical results might be useful for planning numerical and laboratory experiments and for explaining the available data. An important issue here would be to choose small parameters and initial amplitudes facilitating these results which is not a simple problem (c.f. [31]).