Particle Trajectories in Nonlinear Schrödinger Models

The nonlinear Schrödinger equation is well known as a universal equation in the study of wave motion. In the context of wave motion at the free surface of an incompressible fluid, the equation accurately predicts the evolution of modulated wave trains with low to moderate wave steepness. While there is an abundance of studies investigating the reconstruction of the surface profile η\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document}, and the fidelity of such profiles provided by the nonlinear Schrödinger equation as predictions of real surface water waves, very few works have focused on the associated flow field in the fluid. In the current work, it is shown that the velocity potential ϕ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document} can be reconstructed in a similar way as the free surface profile. This observation opens up a range of potential applications since the nonlinear Schrödinger equation features fairly simple closed-form solutions and can be solved numerically with comparatively little effort. In particular, it is shown that particle trajectories in the fluid can be described with relative ease not only in the context of the nonlinear Schrödinger equation, but also in higher-order models such as the Dysthe equation, and in models incorporating certain types of viscous effects.


Introduction
Recent years have seen a flurry of activity aimed at understanding the motion of fluid particles in free surface flows. The problem has been studied from various point of view, including field campaigns, wave flume experiments, asymptotics, and rigorous mathematical analysis. One of the most well known results about particle trajectories in free surface flows is Stokes's finding that the particle motion associated to smallamplitude (linear) periodic waves features a net forward drift which is attributed to the decrease of the Eulerian particle velocity with increasing depth.
Recent advances in laboratory technology, in particular those relating to particle image velocimitry (PIV) and particle tracking velocimitry (PTV) have led to rapid improvements of our understanding of the motion of fluid particles in free surface flow. The influence of the boundary layers on the particle drift in a regular wave train was studied in the seminal paper [23] which was in part inspired by experiments reported in [3]. In essence, particle drift will be positive near the bottom and near the free surface but negative in an intermediate region. In recent works, experimental measurements are coupled with high-performance data analysis techniques to paint a rather complete picture of particle motion in highly nonlinear waves created in a wave flume [8,17,18]. The findings presented in these works are also related to the importance of the effect of the boundary layer both at the bed and at the free surface. Indeed, it is observed that the Stokes drift may take a very different form than what was originally found by Stokes [28]. In particular, these recent studies confirm and extend the results of [23].
However, depending on the motion of the wave maker, these results may vary. For example, the experiments reported in [9] seem to confirm the essence of Stokes's original work. Indeed, using both experiments and high-order asymptotics, a strong case is made in [9] that there is a net forward drift throughout the fluid column. These results are also in line with mathematical advances in the understanding of particle motion in free surface flow (see [16] for a review). In particular, a firm mathematical proof was given that particle trajectories in Stokes waves are not closed [10].
In the case of finite depth, even if viscosity is not taken into account, mass conservation during the creation phase may lead to zero drift when averaged over the fluid column [23,25]. On the other hand, in the current analysis, the case of (infinitely) deep water is investigated and the bottom boundary layer is ignored. In some studies, in particular in the presence of a background current in deep water, it has been observed that there is no (Lagrangian) particle drift either in the average, or even in the pointwise sense [21]. Some reports of field campaigns also point to the absence of the Stokes drift in wave trains in the open ocean [27]. With these partially conflicting pieces of evidence, it appears that there is a need for being able to describe particle paths as accurately as possible using theoretical models. Recent work has focused on the numerical description of particle trajectories based on approximating solutions of the inviscid Euler equations [27,32].
In the current contribution, we study particle trajectories due to wave motion described in the narrow-banded spectrum approximation in an idealized infinitely deep fluid. We focus on irrotational flow which is described in terms of a potential function, φ(x, z, t), satisfying Laplace's equation, φ x x + φ zz = 0, in the fluid domain.
In general terms, if (ξ(t), ζ (t)) denotes the position of a fluid particle in the (x, z)plane at a time t, then the particle motion is described by the ordinary differential equations The corresponding initial conditions are given by (ξ(0), ζ (0)) = (ξ 0 , ζ 0 ) where (ξ 0 , ζ 0 ) represents the initial coordinates of the particle in the (x, z)-plane. In this paper, all variables are dimensional with standard SI units.
As stated above, the focus of this paper is to examine particle trajectories associated with wave motion which can be approximately described by the nonlinear Schrödinger (NLS) equation and a few of its generalizations. This equation arises in the case where a carrier wave of a certain frequency is modulated slowly. In terms of the dimensional slow variables X = x and T = t, the NLS equation is written as where g = 9.8 m/s 2 is the gravitational acceleration, B describes the dimensional amplitude envelope of the oscillations of the carrier wave, k 0 > 0 represents the dimensional wave number of the carrier wave, ω 0 = √ gk 0 represents the dimensional frequency of the carrier wave, and = 2|a 0 |k 0 represents the (dimensionless) wave steepness. This equation has been well studied mathematically (see, for example, [29]) and has been shown to favorably predict experimental measurements of modulated wave trains when < 0.1 (see, for example, [22]).
Given a solution, B(X , T ), the free surface is reconstructed by whereη is a dimensional average term and B 2 is the dimensional amplitude of the first harmonic. These quantities are defined in terms of B and are examined in detail in Sect. 2.1. The abbreviation c.c. stands for complex conjugate. Additionally, in the derivation of the NLS equation, the following ansatz is used for the velocity potential whereφ =φ(X , Z , T ), A = A(X , Z , T ), and Z = z is the dimensional slow vertical variable. Substituting this expression into the Laplace equation leads to the following transport equation with imaginary speed at leading order in . In addition, at leading order in the Bernoulli condition gives Solving this boundary-value problem for A yields the first-order approximation to the velocity field everywhere in the fluid. It is then possible to compute the particle trajectories associated to any given surface profile numerically. A number of examples for particle paths that result from NLS solutions are given in the remainder of this section. Examples of particle paths that result from solutions to higher-order and dissipative generalizations of NLS are examined in the following sections. Similar considerations regarding particle paths and the potential function have been introduced for a variety of long-wave equations in [1,2]. In particular, particle paths and streamlines for waves in the KdV regime were described in detail in [4,19,20]. Somewhat different procedures have also recently been used to understand properties of particle motion in the context of the narrow-banded spectrum approximation in the presence of shear flows [12] and in the presence of point vortices [11].

NLS Plane-Wave Solutions
The plane-wave solutions of NLS are given by where and B 0 and k are real constants. Closed-form expressions for the corresponding surface displacement and velocity potential are included in Appendix A.1. The surface displacement corresponding to this solution has a (temporal) period of For demonstrative purposes, we select = 0.1, B 0 = 1, k 0 = 1 and k = 0 as parameter values for the plane-wave solution. For these parameter values, the period of the surface is given by t * NLS = 100π √ 5/357=1.968, the crest height is 0.223, the trough height is −0.183, and the wave profile has mean zero. The crest and trough heights are not symmetric about z = 0 even though the mean term is zero because η(x, t) is comprised of six (complex) Fourier modes. See Fig. 1 for time-series plots of one period of this NLS plane-wave solution and one period of the corresponding surface displacement.
The paths of NLS plane-wave particles are determined using (7) to determine asymptotic expressions for φ and η that are valid up to O( 3 ) (the details of this process are described in detail in Sect. 2) and then numerically integrating the system of ODEs given in Eq. (1). Figure 2 includes a waterfall plot showing the path of a particle that starts on the surface. Figure 3a shows the fluid surface and the position of a NLS plane-wave particle that starts on the surface at four different t values. Figure  3b shows the paths of three NLS plane-wave particles corresponding to different ζ 0 . The top particle starts on the fluid surface while the other two start (and remain) inside the fluid. Let t NLS (ζ 0 ) represent the period of the vertical motion of a particle with initial z coordinate ζ 0 . Values for t NLS (ζ 0 ) were computed by numerically solving Eq. (1) and determining the period of the resulting solution for a range of ζ 0 values. A plot showing how t NLS (ζ 0 ) and ζ 0 are related is included in Fig. 4a. Note that The horizontal motion of the particles is not periodic because ξ(t NLS (ζ 0 )) > ξ 0 regardless of the initial position of the particle. (Note that this does not mean the particle always moves to the right.) The difference between the final and initial horizontal positions during one period of vertical motion, ξ(t NLS (ζ 0 )) − ξ 0 , is known as the horizontal Lagrangian drift. Figure 4b contains a plot relating the average horizontal Lagrangian velocity, and ζ 0 . As expected from Eq. (18d), both u NLS (ζ 0 ) and the horizontal Lagrangian drift limit to zero as ζ 0 → −∞.
Additional simulations (not shown) establish that increasing B 0 increases t NLS , u NLS and the horizontal Lagrangian drift in roughly exponential manners. Increasing the wave number of the solution, k, increases t NLS , and decreases u NLS and the drift in approximately linear manners. Equation (18c), the kinematic boundary condition, implies that a particle that starts on the surface stays on the surface. This equation is only approximately satisfied by NLS since NLS is an asymptotic approximation to the Euler equations. Therefore NLS   particles that start on the surface of the fluid do not necessarily remain on the surface. A plot of the difference between the particle's vertical position and the fluid surface, is included in Fig. 5. As an additional check on our work, we computed for a variety of values of ∈ (0, 1 4 ]. Table 1 contains a summary of these results that establish E( ) ∼ O( 4 ) for plane-wave solutions of NLS (as expected because NLS is an O( 3 ) approximation).

Cnoidal-Wave Solutions of NLS
The cnoidal-wave solutions of NLS are given by where B 0 and κ ∈ [0, 1) are real parameters. Here cn(·, κ) is a Jacobi elliptic function with elliptic modulus κ and period 4K (κ) where K is the complete elliptic integral of the first kind [5]. Formulas for the corresponding surface displacement and velocity potential are included in Appendix A.2. There are two situations in which B is periodic in T . First, if κ = 1/ √ 2, then B is real for all X and T and has a T -period of 2K (1/ π k 0 , then the solution corresponds to a periodic surface displacement with t-period Plots of the NLS cnoidal-wave solution and the corresponding surface displacement Fig. 6. Plots of the paths of three particles are included in Figs. 7 and 8. Unlike plane-wave solutions, cnoidal-wave solutions do not have constant magnitude. This means that the vertical motion of the particles (except in rare cases) is quasiperiodic rather than periodic. For example, the particle that starts at a peak on the surface has an initial elevation of ζ 0 = 0.1248. The vertical motion of this particle has a quasiperiod of t = 20.2047. After one quasiperiod, the particle reaches an elevation of ζ(20.2047) = 0.1246, slightly lower than its initial height. This is caused by the fact that the particle experienced a horizontal drift of 0.4140 during this interval. Due to this quasiperiodicity, the Lagrangian drift and average Lagrangian velocity are not well defined.
The second situation in which a T -periodic solution is obtained is if B 0 and κ = 1 In order for this solution to correspond to a t-periodic surface displacement, the period of this solution must align with the period of the carrier wave. Enforcing this restriction leads to large-amplitude surface displacements that are well outside of the regime where NLS is expected to be valid. Therefore, for demonstrative purposes, we select = 0.1, k 0 = 1, κ = 0.999 and B 0 = 1 (parameters that lead to a non-periodic surface displacement in the NLS regime). Plots of this solution and the corresponding surface displacement are included in Fig. 9. Figures 10 and 11 show how three of these NLS cnoidal-wave particles move in t. Just as in the previous case, the vertical motion of the particles is quasiperiodic even though the solution is periodic.
For the cnoidal-wave solutions of NLS, increasing B 0 increases the horizontal drift for particles that start on the surface.

Solitary-Wave Solutions of NLS
The solitary-wave solutions of NLS, 3. This solution is not periodic in X or T and therefore does not lead to a periodic surface displacement. We select the parameter values = 0.1, k 0 = 1, and B 0 = 1 for demonstrative purposes. Figure 12 contains plots of this solution and the corresponding surface displacement. Figures 13 and 14 show the paths of three NLS solitary-wave particles on the interval t ∈ [0, 50]. During t ∈ [0, 50], the horizontal drift for the particle that starts on the surface is approximately 0.6170. In the next section, we present the full derivation of the velocity potential in the more general case of the higher-order NLS (Dysthe) equation and the viscous models. Then in Sect. 3, we describe particle paths for a number of examples.

Construction of the Velocity Potential
We now give the full details of the reconstruction of the velocity potential in terms of the unknown B(X , T ) of the nonlinear Schrödinger equations. In fact, we will take a slightly more general view by including viscous effects, and higher-order terms. To explain how to find the potential φ, it is convenient to first review the derivation of the evolution equation.

Derivation of the Viscous Dysthe System
Wu et al. [31] proposed the following system for a two-dimensional, infinitely deep, weakly dissipative fluid (18d) Here φ = φ(x, z, t) represents the velocity potential of the fluid, η = η(x, t) represents the free surface displacement, g represents the acceleration due to gravity, and α > 0 represents dissipation from all sources. The classical Euler equations (see, for example, Debnath [13]) are recovered from this system by settingᾱ = 0. Following the work of Zakharov [32] and Dysthe [15], we make the following modulated wave train ansatz η(x, t) = 3η + Be iω 0 t−ik 0 x + 2 B 2 e 2(iω 0 t−ik 0 x) + · · · + c.c., (19a) φ(x, z, t) = 2φ + A 1 e k 0 z+iω 0 t−ik 0 x + 2 A 2 e 2(k 0 z+iω 0 t−ik 0 x) + · · · + c.c., where ω 0 represents the frequency of the carrier wave, k 0 > 0 represents the wave number of the carrier wave, = 2|a 0 |k 0 represents the (dimensionless) wave steepness, a 0 represents a typical amplitude, and c.c. represents complex conjugate. Further assumē The slow space and time variables are defined by X = x, Z = z, and T = t. Assuming thatᾱ = 2 α and carrying out the perturbation analysis through O( 4 ) leads to (see [6,7] for details) where ω 2 0 = gk 0 and B * represents the complex conjugate of B. If the Hilbert transform, H, is defined by where the Fourier transform and its inverse are defined bŷ then the system (21) can be rewritten as We refer to this system as the viscous Dysthe (vDysthe) system. The NLS equation is obtained from this system by setting α = 0 and disregarding terms of order 2 and higher. The Dysthe system, also known as the modified NLS system, is obtained from Eq. (22) by setting α = 0. The dissipative NLS equation (dNLS) is obtained by setting 2 = 0. Given a solution to Eq. (22), obtaining the corresponding surface displacement is straightforward but tedious. Determining the velocity potential is more complicated (and more tedious) because partial differential equations (PDEs) must be solved to determine each A jk . To this order, the surface displacement and velocity potential are given by (see Sect. 2.2 for more detail) whereB = B(X + i Z, T ),B * = B * (X + i Z, T ), E = e iω 0 t−ik 0 x , and c.c. represents the complex conjugate. Unfortunately, a simple, general formula forφ 1 does not exist. However, we include formulas forφ 1 for the particular solutions we examine below.

Derivation of the Velocity Potential
To determine the Z dependence of each of the terms in the expansion for the velocity potential, PDEs must be solved at each order. At O( ) the PDE and surface boundary condition are At O( 2 ) the PDEs and surface boundary conditions are At O( 3 ) the PDEs and surface boundary conditions are Here is the two-dimensional Laplacian operator, F = F X X + F Z Z . The first of the PDEs listed at each order can be solved as an advection equation with a complex velocity. The remaining PDEs are complex advection equations with nonhomogeneous terms that can be solved using variation of parameters. For brevity, the details have been omitted, but the solutions to these systems are included in the Eq. (23).

Particle Paths for Generalizations of NLS
We now present particle paths associated with three generalizations of the NLS equation. We start with the higher-order non-viscous model, the Dysthe equation. Then we continue with the dissipative nonlinear Schrödinger equation (dNLS) and finally with the viscous Dysthe equation.

The Dysthe System
The Dysthe system is obtained from Eq. (22) by setting α = 0. This system has been shown to provide accurate predictions for the evolution of modulated wave trains for a wider range of values than NLS (see, for example, Lo and Mei [22]). Its plane-wave solutions are given by where and B 0 , k are real constants. The formulas for the corresponding surface displacement and velocity potential are included in Appendix A.4. For demonstrative purposes, we select = 0.1, B 0 = 1, k 0 = 1 and k = 0. These parameters lead to a surface displacement with crest height of 0.2247, trough height of −0.1813, and a t-period of This particular period is the same as in the NLS example examined above because we selected k = 0 in both cases. Figure 15 contains time-series plots of one period of the real and imaginary parts of the Dysthe plane-wave solution and one period of  Figure 16 shows the paths followed by three Dysthe plane-wave particles. One starts on the surface and two start inside the fluid. Figure 17 contains plots demonstrating how the period and mean horizontal velocity depend on ζ 0 for this Dysthe plane-wave solution. For a given surface amplitude and ζ 0 , the periods, horizontal drifts, and average horizontal velocities of the Dysthe particles are all greater than the corresponding NLS quantities. The differences between the NLS and Dysthe results are small due to the orders of the equations, and these differences go to zero as z → ∞. Finally, for these Dysthe plane-wave solutions, the difference between the particle's vertical position and the fluid surface, as defined in Eq. (13), decreases like O( 5 ).

The dNLS Equation
The dissipative nonlinear Schrödinger equation is obtained from (22) by setting 2 = 0. The dNLS equation was first derived as a model of water waves by Dias et al. [14]. However, Lo and Mei [22] and Segur et al. [26] had previously studied it and had shown that it accurately models the evolution of modulated wave trains. The plane-wave solutions of dNLS are given by where and B 0 and k are real constants. These solutions were chosen so that they limit to the plane-wave solutions of NLS as α → 0. Formulas for the corresponding surface displacement and velocity potential are included in Appendix A.4. The plane-wave solutions of dNLS are not periodic in T , so the corresponding surface displacements are not periodic in t. For demonstrative purposes, we selected = 0.10, B 0 = 1, k 0 = 1, k = 0, and α = 4. This value of α is roughly an order of magnitude larger than the experimentally determined parameters used by Segur et al. [26] and Carter and Govan [6]. We used a relatively large value for α so that viscous effects could be seen on relatively short time scales. Figures 18 and 19 show the paths of three dNLS plane-wave particles on the interval t ∈ [0, 75]. The top particle starts on the surface and the other two start and stay inside the fluid. Although the motion of the particles is not periodic, a number of comments can still be made. Due to dissipative effects, the motion, both horizontal and vertical, of all particles decreases as time increases. Each particle eventually spirals in toward a fixed point. The "period", the average horizontal Lagrangian velocity, and the horizontal drift all decrease as t increases. This is consistent with the NLS results where smaller-amplitude solutions lead to smaller periods, velocities, and drifts. Figure 20 contains a plot of D(t) versus t for this plane-wave solution of dNLS. Note that D(t) does not limit to zero as t → ∞. This is because the particle that starts on the surface ends up inside the fluid. This is due to two facts: (i) dNLS is an asymptotic model and (ii) the weakly viscous Euler equations are only an approximation to the true viscous system. However, dNLS remains valid in the small viscosity,ᾱ → 0 limit. Both D(t) and E( ) limit to zero as → 0. Plots of a ξ(t) and b ζ(t) for the three dNLS plane-wave particles in Fig. 18

The Viscous Dysthe Equation
The viscous Dysthe system is given in Eq. (22). The plane-wave solutions of this system are given by where  where B 0 and k are real constants. These solutions limit to the plane-wave solutions of the Dysthe system as α → 0. The plane-wave solutions of vDysthe are not periodic in T . For demonstrative purposes, we selected = 0.1, B 0 = 1, k 0 = 1, k = 0, and α = 4. This value of α is roughly an order of magnitude larger than the experimentally determined parameters used by Segur et al. [26] and Carter and Govan [6]. We used a relatively large value for α so that viscous effects could be seen on relatively short time scales. Figures 21 and 22 show the paths of three vDysthe plane-wave particles on the interval t ∈ [0, 75]. There are many similarities with the dNLS case. As t increases, the "period", the average horizontal Lagrangian velocity, and the horizontal drift all decrease. The horizontal drift and average horizontal velocity limit to zero as t → ∞. The particles spiral in toward a fixed point. Figure 23 shows that the error term D(t) is not smaller in the viscous Dysthe context than in the dNLS context. This observation further demonstrates the limitation that Eq. (18c) is only valid in the small α limit. Using a higher-order approximation to Eq. (18) does not provide a more accurate approximation to Eq. (18c) because Eq. (18c) with viscosity is only an approximation to the true viscous system.

Conclusion
In this paper, we used the classical derivation of the NLS equation to derive formulas for the velocity potential throughout the fluid. Using these formulas, we numerically computed and examined the paths of fluid particles corresponding to plane-, cnoidal, and solitary-wave solutions to the NLS equation. Following a similar procedure, we examined the paths of particles subject to the motion of plane-wave solutions to the Dysthe, dissipative NLS, and viscous Dysthe equations. We showed that dissipative/viscous effects decrease particle speed and displacement. Finally, we showed that the boundary conditions of the full water wave problem are only asymptotically satisfied by the solutions to these equations.

A.2 NLS Cnoidal-Wave Solutions
The surface displacement and velocity potential for the NLS cnoidal-wave solution given in Eq. (14) are where Note that the mean terms are nonzero, and that the velocity potential for this solution has poles at z = O( −1 ) because the Jacobi elliptic function is evaluated with a complex argument. These poles occur well below the surface of the fluid. Therefore, we do not consider them. The particle paths we consider in Sect. 1.2 are sufficiently close to the fluid surface (and are far from the poles).
Note that the mean terms are nonzero, and that the velocity potential for this solution has poles at z = O( −1 ) because the hyperbolic secant is evaluated with a complex argument. These poles occur well below the surface of the fluid. The particle paths we consider in Sect. 1.3 are sufficiently close to the fluid surface, and are far from the poles.

A.5 dNLS Plane-Wave Solutions
where ω r and ω i are given in Eq. (32). Note that the mean terms are nonzero. Note in addition that the nonlinear term in dNLS causes some portions of the velocity potential to decay faster than others.

A.6 Viscous Dysthe Plane-Wave Solutions
where ω r and ω i are given in Equation (34). Just as with all the other plane-wave solutions, the mean terms here are zero.