Qualitative Description of the Particle Trajectories for the $N$-Solitons Solution of the Korteweg-de Vries Equation

The qualitative properties of the particle trajectories of the $N$-solitons solution of the KdV equation are recovered from the first order velocity field by the introduction of the stream function. Numerical simulations show an accurate depth dependance of the particles trajectories for solitary waves. Failure of the free surface kinematic boundary condition for the first order type velocity field is highlighted.


Introduction
The Korteweg-de Vries (KdV) equation,  without losing its shape, Russell published a report of his experiments and observations in 1844 [20]. The Korteweg-de Vries equation was later derived independently by Boussinesq in 1877 [4] and by Korteweg and de Vries in 1895 [19].
The asymptotic behaviour of solutions of (1) arising from smooth and fast decaying initial disturbance obtained by the inverse scattering transform (Gardner, Greene, Kruskal and Miura, ( [13,1967], [14,1974])), confirmed the numerical experiments of Zabusky and Kruskal ( [22,1965]): a finite number of solitons are travelling to the right while an oscillating and decaying wavetrain disperses to the left.
A soliton of (1), the approximation in the KdV regime of a solitary wave, is given by where a > 0, s ∈ R. The following properties hold for (2) 1. The amplitude is given by a and is reached at x = s − √ gh 0 (1 + a 2h0 )t; 2. The traveling speed is √ gh 0 (1 + a 2h0 )t; 3. The taller the soliton, the narrower.
The behaviour of N solitons is given by the N -solitons solution, representing the KdV approximation of the interactions of N solitary waves propagating in the same direction. They behave like N distinct solitons travelling at different constant speed except when interactions occur, in which case they are left unscathed from an interaction except for a phase shift. Moreover, if two solitons with the same order of amplitude interact, they exchange their mass and speed during the interaction. The Korteweg-de Vries equation not only provides an approximation of the surface elevation in shallow water but admits as well an approximation of the underlying particle trajectories. To the first order approximation in the perturbation parameter, one obtains, with the divergence-free condition, the first order velocity field of (1), Since the KdV equation also holds to the second order, one recovers a second order approximation of the horizontal velocity field at the flat bottom ( [21]), Thus, one may study the particle trajectories, for a given velocity field, associated to solutions of KdV by solving, In the case of a soliton, the particle paths obtained with the first order velocity field (3) fail to capture one essential property of the particle trajectories observed experimentally for solitary waves: the higher in the fluid a particle is initially located, the larger its horizontal displacement [5]. It is a consequence of the independence on the y variable of the horizontal component of the velocity field. Since the soliton solution is considered as an accurate approximation of the solitary waves solution of the full governing equations ( [9]), it is crucial to seek for higher approximations of the velocity field in order to derive better approximation of the particle trajectories. Figure 2 illustrates the experimental results of [5] for different solitary waves.
In the literature, two methods have been developed to recover a higher order approximation of the velocity field for the KdV equation. Introducing a stream function for the first order approximation and using the fact that the velocity potential and the stream function are harmonic conjugate, Constantin ([7]) and Henry ([15]) have obtained a higher order velocity field for a soliton and periodic solutions respectively. Qualitatives results on the particle paths were obtained in both articles.
A careful analysis of the derivation of the KdV equation in the second order approximation allowed Borluk and Kalisch to recover from (4) a nonlaminar velocity field underneath the surface ( [3]). They obtained numerically a better path of the particles in the case of the soliton solution, the cnoidals solution and the 2-solitons solution. A similar work in the case of the Boussinesq equation had been done before in [1,2].
The main result of this article, following the work of Constantin and Henry ([7,15]), is the qualitative description of the particle trajectories for the N -solitons solution of KdV. Theorem 1.1 Let N ∈ N denotes the number of solitons and a i , 1 ≤ i ≤ N , denotes the amplitude parameter of each solitons. For sufficiently small amplitude parameters i := a i /h 0 , the particle trajectories described by the velocity field satisfy the following: (a) A particle on the flat bed move in straight line at a positive speed; (b) All the particles below the surface have a positive horizontal speed; Moreover, if a particle doesn't encounter soliton interactions, we have the following: (c) For each such particles, there exists t 1 , . . . , t 2N −1 such that the particle move upward for, and downward for, (d) The particle's position at times t 1 , t 3 , . . . , t 2N −1 follows the phase shift of the solitons.
We would like to emphasize the nonlinear nature of this result. As it will be pointed out in Section 3, the phase shifts of the solitons are the result of the nonlinearity of the KdV equation. Therefore, one cannot obtain the right behaviour of the particles for the N -solitons solution with the superposition principle used on the 1-soliton velocity field obtained in [7]. Moreover, there exists no linear approximation of solitons. This fact is of relevance when comparing the results to the periodic setting. Indeed, for the periodic travelling waves of the full governing equations, a linear approximation of the governing equations leads to nonlinear equations for the particles paths ( [12]) and the predictions of the linear theory are reflected by what goes on in the full nonlinear system ( [6,11]). Moreover, note that as the period of the periodic travelling wave of KdV increases towards infinity, one recovers the 1-soliton ( [8,17]).
The outline of the paper is the following. In Section 2, inspired from [7,10], we recall the derivation of the Korteweg-de Vries to the first order in the perturbation parameter. In Section 3, we present the Nsolitons solution and we use the method of [7,15] in order to prove Theorem 1.1. This is done by noticing that this method is not limited to travelling waves. Finally, a numerical analysis of the particle trajectories is provided in Section 4. We show that the monotonicity in the y variable of the horizontal displacement is recovered. Moreover, the numerical investigation indicate the possible lost of the free surface kinematic boundary condition for velocity fields of KdV.

Derivation of the Korteweg-de Vries equation
The physical assumptions to derive the Korteweg-de Vries equation are the following. Consider approximately two-dimensional waves, that is, the wave profile is the approximately the same for all cross sections taken with respect to the crest line. Assume irrotational waves in a homogeneous, non-viscous and incompressible fluid. Moreover, assume the impermeability of the bottom (bottom kinematic boundary condition), that the same particles form the surface (free surface kinematic boundary condition) and that the motion of the air and that of water is decoupled. Then, the full governing equations are given by Let us first set the non-dimensional full governing equations. Consider the perturbation of the pressure, which express the change of the pressure as the waves move, and the change of variables, where λ is the typical, or average, wavelength and where x → λx is to be understood as x replaced by the new variable λx. Introducing the shallowness parameter, the nondimensional full governing equations reads The amplitude of the waves plays a fundamental role in the formulation of the water-wave problem. Thus, let a denotes the typical, or average, amplitude of a wave and express the wave profile by η → aη. Then, defining the amplitude parameter by = a/h 0 , the nondimensional full governing equations become, Since p and v scale like , we rescale them as well as u for consistancy, Then, The regime = O(δ 2 ) arises by requiring that the above (x, t)-variables remain of the same order. One removes the shallowness parameter and symmetrizes the equations by considering the change of variables The fourth equation of (6) allows us to define a velocity potential φ(x, y, t) such that By restricting ourselves to waves travelling from left to right and by considering a slow time scale, we obtain the following for φ, Consider the following perturbative expansions in the amplitude parameter, The zero order approximation yields, The first three equations implies, φ 0 y ≡ 0, 0 ≤ y ≤ 1, while the boundary condition on y = 1 allow us to deduce, Using the first and second order of (7), one obtains that η 0 satisfies the Korteweg-de Vries equation The region where the KdV balance occurs is when ξ = O(1), T 0 = O(1), that is, in the physical variables, when Coming back to the physical variables, one obtains the Korteweg-de Vries equation (1) and a laminar flow for solutions η of (1) 3 A higher velocity field for the N -soliton solution 3.

The N -solitons solution
Let us recall the N -solitons solution as introduced by Hirota in [16] (see also [21]) for the Korteweg-de Vries (1). Let where F is given by the power serie expansion, and Then, by replacing (10) in (1), one obtains that only the first N order term in (11) are nonzero. Moreover, one may express F as , is the product over these n indexes, provided k < l when specified. It will be convenient to denote Let us explain why this solution is called the N -solitons solution. Consider 0 < a N < ... < a 1 and s N < ... < s 1 . From (10), η can be recasted as First, suppose that we are in the region where f 1 1 and f k 1 for all k ≥ 2. The behaviour of η is given by, that is, a soliton of height a 1 and phase s 1 . Suppose now, for 2 ≤ k ≤ N , that we are in the region where f k 1, f i 1 for 1 ≤ i < k and f i 1 for k < i ≤ N . Then, that is, a soliton with a phase shift of 1 2β k ln . From the last computations, one sees that the N -solitons solution behaves like N distinct solitons when they are far apart from each other. Moreover, interactions between solitons must occur since the solitons are travelling at different speed. Indeed, a direct computation shows that a i = a j , otherwise the solution becomes a (N − k)-solitons solution, k being the number of times a i = a j , i < j. Thus, up to a change of indexes in the last computations, one sees that they are left unchanged in shape after interactions, the only notable effect of the interaction being the phase shift. Figure 3 illustrates the phase shift resulting from the interaction of 2 solitons.

Nonlaminar velocity field
The velocity potential obtained in the previous section is independent of the variable y. Therefore, even if the soliton solution of the KdV equation provides a good approximation of the solitary wave, one cannot describe accurately the underlying structure of the particles paths.
Let us explain how one may recover a nonparallel displacement of the particles. Consider that φ is the velocity potential of (5) satisfying, φ x = u, φ y = v.
We introduce the stream function ψ of (5) satisfying For irrotational flows, the velocity potential and stream function are harmonic conjugates. Consider the complex velocity potential f , with t acting as a parameter, such that f (z, t) := φ(x, y, t) + iψ(x, y, t), z ∈ C, with f (z, t) = u(x, y, t) − iv(x, y, t).
Instead of using (8) to obtain the expression an expression of the velocity potential independant of the y variable, we rather solve the Laplacian equation of the velocity potential and of the stream function in the upper-half plane {z = x + iy ∈ C|y > 0} with the velocity field (3) as a boundary condition on y = 0.
If the first boundary condition is the natural non penetration condition, we use the independance of the y coordinate of the horizontal component of the velocity field to obtain the second boundary condition ( [7,15]).
The boundary condition implies By the reflection principle, the analytic function f (z) has an unique extension across y = 0. Therefore, f is defined up to a constant. Using Alembert's formula, we obtain the following expressions of u and v u(x, y, t) = 1 2 In the case where η is the N -solitons solution, an explicit expression of the velocity is obtained. Indeed, let, where, Then, one has, u(x, y, t) = 1 2 On the other hand, v(x, y, t) = i 2 the last line coming from the expression of the derivatives of F c and F s , 2β ij y .

Proof of Theorem 1.1
Before proving Theorem 1.1, let us introduce some notations to simplify some of the expressions. We denote , m , for m sums, and m , for m products. We shorten the expression of α(i 1 , ..., i n ) by α in and define α in,jm := α in α jm and so on for product over more than two terms. Finally, we use the slightly modified Einstein notation (β in ) for a sum from β i1 to β in , while β ij denotes a single coefficient. Since no confusion is possible, (β in )(β jm ) denotes the product over the two sums.
Let us turn ourselves to the proof of Theorem 1.1 Proof: (a) From (16)- (17), we recover the original velocity field when y = 0.
The coefficients of the first sum with respect to β i are all positive. Moreover, the coefficients appearing in the other sums are non-negative. Indeed, consider for example the coefficient in front of cos((2β in )y) cos((2β jm )y) sin((2β kp )y) sin((2β lq )y) : Therefore, since in order to prove (b), it is sufficient to impose, (c) Let us give an equivalent form of v. Consider first f cos(2(β in )y) sin(2(β jm )y).
For sake of simplicity, assume that 0 < a n < . . . < a 1 . From the previous expression, the behaviour of v when x → ∞ is governed by f N , which is positive in the neighborhood of x = s 1 + c 1 t if and only if x − c 1 t > s 1 . The neglected terms in the previous analysis don't affect the nature of the change of sign but rather, and slightly, its position.
In the case where f k 1, f i 1, i < k and f i 1, i > k, the leading terms of the numerator are given by We obtain that this last expression is positive in the neighbourhood of x = s k + c k t if and only if x − c k t > s k + 1 2β k ln(α(1, ..., k)/α(1, ..., k − 1)). The phase shift here reproducing the phase shift in the N -solitons solution shown in (14). Here again, the neglected terms may induce a small difference between the phase shift in this analysis and the one corresponding to (14) but the nature of the change of sign of v is preserved.
From the previous analysis, we obtain c). Indeed, v > 0 when x → ∞ and v is still positive from x → ∞ to a small neighborhood of x = s 1 + c 1 t, representing the position of the crest of the first soliton. From the neighbourhood of x = s 1 + c 1 t to the neighborhood of x = s 2 + c 2 t + 1 2β2 ln(α(1, 2)), only (19) and (20) play an important role in the sign of v. The functions f i being monotonic, only one change of sign of v occurs in this region and v therefore become positive for values of x greater than the small neighborhood of x = s 2 + c 2 t + 1 2β2 ln(α (1, 2)). The analysis of (20) shows that there a change of sign in the neighborhood of the crest of the second soliton. The same analysis can be carried through all the crests of the N -solitons solution until the final crest, that is, at the left of the last soliton where v is always negative.

Numerical analysis of the particles trajectory
This section is devoted to the numerical analysis of the particle trajectories given by the higher approximation of the velocity field (16)- (17). First, we compare the results in the case of a single soliton with the numerical approximation of the particle trajectories of the first order approximation field (3). We use the experimental results in [5] as the reference trajectory. The simulation were done using a fourth order Runge-Kutta scheme. Figure 4 represents the numerical approximation of the particles trajectory with both velocity field with the parameters of experiment (c) of [5]. Table 5 compares the total displacement in the x variable and the maximal displacement in the y variable of the velocity fields with the experimental results of [5]. From the numerical results, one recovers an important feature of the particle trajectories of solitary waves : the higher the particles are initially located, the greater the horizontal displacement. However, one notices for both velocity fields that the vertical displacement of the particle located at the surface is greater that the height of the soliton. It comes from the fact that the free surface kinematic boundary conditions assuring that the particles stay underneath the water surface is not preserved in (3) and therefore in (16)- (17). Indeed, this boundary condition is a higher order term in the approximation parameter as one can see from the third line of (7). Numerical simulations show that the overshoot is proportional to . The numerical approximations of the particles trajectory are therefore relevant for small values of and may be interpreted qualitatively for larger values of . Figure 6 illustrates the numerical approximation of the particles trajectory associated to the higher velocity field in the case of the 2-solitons in the case where