Bound Coherent Structures Propagating on the Free Surface of Deep Water

This article presents a study of bound periodically oscillating coherent structures arising on the free surface of deep water. Such structures resemble the well known bi-soliton solution of the nonlinear Schrödinger equation. The research was carried out in the super-compact Dyachenko-Zakharov equation model for unidirectional deep water waves and the full system of nonlinear equations for potential flows of an ideal incompressible fluid written in conformal variables. The special numerical algorithm that includes a damping procedure of radiation and velocity adjusting was used for obtaining such bound structures. The results showed that in both nonlinear models for deep water waves after the damping is turned off, a periodically oscillating bound structure remains on the fluid surface and propagates stably over hundreds of thousands of characteristic wave periods without losing energy.


Introduction
Nonlinear coherent structures arising on the free surface of fluid are of significant interest to oceanology. Today, there is a large number of theoretical models for describing such structures. In the case of shallow water it is the most famous Korteweg-De Vries (KdV) equation [1] and its various modifications and generalizations [2][3][4]. An undoubted advantage of studying such nonlinear objects in shallow water (in this case solitons or breathers) is the presence of a large number of models integrable in terms of the inverse scattering problem [5,6] and availability of exact analytical solutions.
In the case of deep water waves the only integrable model is the well-known Nonlinear Schrödinger equation (NLSE) [7]. The NLSE also has various exact solutions [8][9][10][11] and generalizations [12]. However, in practice, using these models may not be sufficient to describe the dynamics of deep water waves and various nonlinear wave phenomena since these models are derived in the assumption of spectral width narrowness which significantly limits their range of applicability. In this case, it is necessary to use more accurate models such as the super-compact Dyachenko-Zakharov equation (SCDZE) [13] or the full system of nonlinear equations for an ideal potential fluid written in conformal variables [14,15]. In [16], the authors showed that the full model also has a solution in the form of a giant breather. The classical NLSE soliton was used as the initial condition for obtaining such a solution. This discovery was the starting point for further experimental observation and research of such breathers on the surface of deep water [17][18][19][20]. It was shown in [21] that breather solution can be found numerically by using the Petviashvili method [22] in the model of SCDZE.
The SCDZE can be rewritten in terms of an envelope wave packet [23]. Under the assumption that the spectral width is narrow it can be reduced to the NLSE, and the SCDZE breather would then correspond to the NLSE envelope soliton. However, despite some similarities, the behaviour of breathers and solitons can differ significantly. For example, it was shown in work [24] that there are significant differences in the interaction of a pair of breathers and a pair of solitons. In the SCDZE there is an energy exchange between breathers depending on their relative phase. Moreover, they change their amplitudes and velocities during the interaction, which is not observed for solitons in the fully integrable model of NLSE.
It should be noted that, so far, in the SCDZE, only stable solutions in the form of a single breather are known. However, an unexpected result was obtained in [25], where multiple collisions of breathers in the periodic domain were studied. It was shown that the only one coherent structure remained at the end regardless of the breathers number and initial parameters. Still, the scenarios of long-time dynamics of the system could be different. In some cases, breathers were coupled into a single periodically oscillating structure resembling the well known NLSE bi-soliton solution [8,26]. This coupled solution allowed us to assume that the more accurate models also have solutions similar to the NLSE bi-solitons, and these solutions can be found numerically or even analytically. This work is devoted to the search for such solutions, as well as to the study of their characteristics.

Deep Water Wave Equations
The fully nonlinear 2D equations describing gravity waves on the surface of ideal deep fluid are well known and can be written as: Here x and y are the horizontal and vertical coordinates, t is time, g is the free-fall acceleration, η(x, t) is the shape of the surface, φ(x, y, t) is the hydrodynamic potential inside the fluid.
One-dimensional potential flow of an ideal fluid of infinite depth in the presence of gravity is a Hamiltonian system. As was shown by Zakharov [27], the surface elevation η(x, t) and the velocity potential at the surface ψ(x, t) = φ(x, y, t)| y=η of the fluid are canonically conjugated variables and satisfy the following Hamilton's equations: Here H is the Hamiltonian, that is, total energy of the fluid: In order to obtain and further investigate the bound coherent structures, we will use two models for deep water waves in our work-the super-compact Dyachenko-Zakharov equation (SCDZE) obtained from truncated Hamiltonian when considering waves propagating in the same direction (Section 3) and the full system of nonlinear equations for an ideal potential incompressible fluid written in conformal variables (Section 4).

Bound Structures in the Super Compact Dyachenko-Zakharov Equation
Assuming that wave steepness is small, the Hamiltonian can be represented as an infinite power series of η and ψ (see [27]). We consider this series up to the fourth order term: Here, operatork means multiplication by |k| in k-space.
Equations of motion (2) with truncated Hamiltonian (4) can be used for numerical simulations of water surface dynamics [28]. However, they are not convenient for analytic study because η(x, t) and ψ(x, t) are not optimal canonical variables. Since the considered system is Hamiltonian, one can choose better Hamiltonian variables by performing a proper canonical transformation. In 2017, the most optimal canonical transformation from the pair of classical physical Hamiltonian variables η(x, t) and ψ(x, t) to normal complex variable c(x, t) was found. It not only cancels non-resonant third and fourth-order terms, but greatly simplifies the expression for the Hamiltonian when considering unidirectional waves (see [13]). This canonical transformation can be written as power series of c k up to the third order: and is described in detail in [13]. As a result of the transformation the Hamiltonian takes the following form in x-space: Here, operatorV k is multiplication by ω k k in Fourier space. Then, the equation of motion or the super compact Dyachenko-Zakharov equation (SCDZE) for unidirectional deep water waves has the following form: Here, operatorω k is a Fourier multiplier by the gk, where g is free fall acceleration. Operator ∂ + x corresponds to ikθ(k) in Fourier space, where θ(k) is Heaviside step function. Physical variables of surface η(x, t) and potential ψ(x, t) can be recovered by canonical transformation (see details in [13]) the first and the second order terms of which can be written in x -space in a compact way: Here,Ĥ is Hilbert transformation with eigenvalue isign(k) and operatorsk α act in Fourier space as multiplication by |k| α .
The SCDZE has a breather solution in the following form: Breather solution (9) is determined by two parameters: k 0 is a carrier wave number, and δ. It can be found numerically only by the Petviashvili method (see details in [13]). It was found in [21] that the SCDZE is not integrable so that the collision of a pair of breathers is not purely elastic and is accompanied by an exchange of energy between breathers and the radiation of incoherent waves (see earlier papers [21,29], and also [20,24]).
It was observed for the first time in [25] that a collision of two single breathers with very close velocities could result in the formation of a periodically oscillating structure resembling the NLSE bi-soliton. Therefore, the obvious step to obtain the bound structure is to take two single breathers and then set them at the same point in space as initial conditions.

Numerical Method for Obtaining Bound Structures. Experiments with Two Breathers
Numerical simulations were carried out in the periodic domain with length L. The SCDZE was solved using the pseudo-spectral Fourier method, and the fourth-order Runga-Kutta method was applied to calculate time evolution. The FFTW3 library [30] was used for the fast Fourier transform procedure. The multiplication of grid functions was carried out in x-space; the direct and inverse Fourier transform were used for calculating derivatives and non-local terms.
The algorithm used for finding a solution is very similar to the one described in [16]. The main defining parameters were wave steepness µ and amplitude ratio κ = C 1 C 2 , where C 1 and C 2 -amplitude of the corresponding breather. In the process of calculating the time evolution the initial structure begins to radiate incoherent waves. Consequently, the steepness of the resulting structure will always be less than the steepness of the initial one. Thus, it is necessary to distinguish between initial steepness and resulting steepness. We will denote them as µ ini and µ f in respectively. Here we also note that we define amplitude ratio κ only at the initial moment of time. The damping procedure was used to remove radiation and for this, a term of the following form was added to the right-hand side of SCDZE: .
Selected function f (x) and coefficient γ concentrate the damping at the edges of the domain without affecting the structure itself, which is always in the center of the domain. To keep the solution in the center, the velocity of the reference frame was constantly adjusted. Coefficient γ should not be chosen too high to prevent numerical instability. On the other hand, it should also not be too small; otherwise, the resulting radiation will not be damped. We chose γ in the range of 10 −2 -10 −3 in our experiments.
As mentioned earlier, we first set the initial conditions in the form of two single breathers found by the Petviashvili method and placed them at one point of the domain. As one can see from Figure 1, an oscillating structure is formed with time and stably exists for 10 6 seconds for any amplitude ratio κ. Initial steepness µ ini was 0.2 while resulting steepness µ f in reaches ≈0.18. We plot these profiles after the damping process is switched off, the moment in time t itself does not matter in principle since these structures are periodic in time. The profiles are plotted in approximately the same phase to be compared with each other. The damping process itself lasted 3 × 10 5 seconds, damping coefficient γ was equal to 10 −3 . The radiation level does not exceed 0.01 or 0.5% of maximum amplitude. The value of amplitude ratio κ does not significantly affect the resulting structure. A further increase in the breathers amplitudes results in a significant amplitude amplification during their interaction and wave breaking.
The full long film with the procedure of damping and further dynamics of the bound structure with µ ini = 0.2, κ = 1.5 and µ f in ≈ 0.18 up to the 10 6 seconds can be found at the following link (accessed on 11 March 2021): http://kachulin.itp.ac.ru/mdpi-2021-films/ SCEq-2br-mu018-l.mp4.

Numerical Method for Obtaining Bound Structures. Experiments with Bi-Solitons of the NLSE
Equation (7) can be re-written in terms of envelope of wave packet C(x, t). If we suppose that the Fourier spectrum of waves has a maximum at k = k 0 we can introduce envelope function C(x, t): c where ω k 0 = gk 0 is the corresponding linear frequency. Using Equation (10), which is actually canonical transformation, one can easily derive the exact equation for the envelope C that we named the "Dyachenko-Zakharov envelope equation" [23]. This equation written in the reference frame moving with the group velocity V 0 = ω k 0 2k 0 has the following form: The equation was derived without any assumptions about the spectral width of the wave packet. Moreover, Equation (11) is the exact envelope form of Equation (7) and has the same range of applicability. One can extract the nonlinear Schrödinger (NLS) equation in the assumption of narrowband spectrum:

∂C ∂t
The NLSE bi-soliton solutions ( [8,26]) are well-known examples of bound nonlinear coherent structures. Thus, based on the above approximation, another way to set the initial wave field is to use the NLSE bi-solitons the explicit expression for which has the following form: The solution is periodic in time with T = 4π k 2 0 (C 2 1 +C 2 2 ) , C 1 and C 2 -solitons amplitude. We introduce also the variable κ = C 1 C 2 which plays a role in the dynamics of the resulting structure. This parameter is also defined only at the initial moment of time.
The examples of used as initial conditions bi-solitons with µ ini = 0.1 and different κ values are shown in Figure 2. The plotted profiles correspond to Equation (13) at t = 0, and one can see that there are two simultaneous zeros of |c(x)| at such a phase. The calculation process is not much different-the bi-soliton begins to radiate, it is suppressed by damping procedure and then it turns off. Figure 3 shows |c(x)| for the bound structures resulted from the initial NLSE bi-solitons with different amplitude ratios κ and initial steepness µ ini = 0.1 (a), 0.2 (b), 0.3 (c) and 0.4 (d). As can be seen, structures exist in a wide range of defining parameters µ and κ. The bound objects stably propagate for at least hundreds of thousands of wave periods T = 2π w 0 ≈ 10s. The damping coefficient γ was 10 −3 . The duration of damping was chosen in such a way as to control the radiation level in the range of 10 −2 -10 −3 . As µ ini increases, the influence of κ becomes insignificant. In particular, in the case of µ ini = 0.4, the variable of κ does not affect the profile of the resulting structure at all. Another feature is that increasing µ ini and κ results in a greater deviation from the original bi-soliton solution. Note that in all obtained structures two zeros of |c(x)| can occur only alternately, but not simultaneously, unlike the NLSE bi-soliton solution.  The full long film with damping procedure and further dynamics up to 10 6 seconds of the bound structure resulted from the NLSE bi-soliton with µ ini = 0.3, κ = 1.5 and µ f in ≈ 0.25 can be found at the following link (accessed on 11 March 2021): http://kachulin. itp.ac.ru/mdpi-2021-films/SCEq-mu025-l.mp4. The small window at the top right shows the dynamics of the structure on a logarithmic scale to see the radiation level.
An interesting feature when using bi-soliton as an initial condition is that there is a range of parameters µ and κ at which obtaining a bound structure becomes impossible. Instead, there is a split of bi-soliton in two separate solitons (Figure 4). The phenomenon manifests itself at small amplitude ratios (κ < 1.25) depending on the steepness, but at κ = 1.1 bi-soliton splits at any value of µ within the applicability of the SCDZE. In the case of initial conditions in the form of two single breathers, no separation occurs.  The film with NLSE bi-soliton splitting can be found here (accessed on 11 March 2021): http://kachulin.itp.ac.ru/mdpi-2021-films/Bi-Sol-split.mp4.

Bound Structures in the Full System of Nonlinear Equations
Following work [14] we perform conformal mapping of the free surface liquid domain z = x + iy confined by a free boundary y = η(x, t) onto a half-plane of new complex variable w = u + iv having a fixed boundary v 0 --see Figure 5. This transformation can be written in the following form: y (x,t) Then, we use the following variables suggested by Dyachenko [15]: We define functions U and B using the Hilbert operator and the projection operator P = 1 2 (1 + i H) as: Here and below the asterisks stands for complex conjugation. In new variables (15) the fully nonlinear (RV) equations have the following form: with the boundary conditions: Figure 6 shows the structure obtained from the original bi-soliton with µ ini = 0.2 and κ = 1.5. The maximum steepness of the obtained structure µ f in reaches ≈ 0.175. The structure propagated stably for 2 × 10 6 seconds without losing energy in the SCDZE model (a). After completing the calculation in the SCDZE we used canonical transformations (8) from variable c(x) to variables η and ψ. Further, variables R and V were obtained using (15), and we continued the simulating process in the model of RV equations (b). The structure also propagated stably for 1 × 10 6 seconds within this model.  The bound structure can also be obtained in the exact RV model by specifying two breathers or the NLSE bi-solitons as the initial conditions. Figure 7 shows the example of bound structure with µ f in ≈ 0.25. It was obtained using the NLSE bi-soliton with µ ini = 0.3 and κ = 1.5 as the initial condition. Panel (a) shows the surface profile η(x) while panel (b) shows a squared surface η 2 (x) in logarithmic scale to monitor the radiation level. As in other cases, the structure exists stably and does not lose energy. If initial steepness µ ini exceeds 0.5, it will always result in the splitting of the bound state mentioned earlier regardless of the initial conditions and the choice of κ.  The full long film with the procedure of damping and further dynamics of this structure up to the 10 6 seconds in the model of RV equations can be found at the following link: http://kachulin.itp.ac.ru/mdpi-2021-films/RV-mu025-l.mp4.
The short film with detailed dynamics of bound structure during one period in the model of RV equations can be found here: http://kachulin.itp.ac.ru/mdpi-2021-films/RV-mu025-s.mp4.

Conclusions
We investigated bound coherent structures in the model of SCDZE for unidirectional deep water waves and the full system of nonlinear equations for potential flows of an ideal incompressible fluid written in conformal variables. Two single breathers placed at one point of the domain, as well as well-known NLSE bi-soliton solutions were used as initial conditions in numerical experiments. In the simulations the initial structure radiated incoherent waves that were suppressed by the damping mechanism. It was shown in both variants of initial conditions that after the damping procedure is turned off, the remaining bound structure propagates stably on the surface of the fluid and does not lose energy for hundreds of thousands of characteristic wave periods. Such structures were obtained in a wide range of determining parameters-steepness and amplitude ratios. It was also revealed that there is a range of parameters at which it is impossible to obtain the bound structure. Instead, there is a split of initial structure into two separate ones. The resulting bound structures may turn out to be new solutions to the equations used and we hope that further research will help understand their nature in more detail.