Water-wave gap solitons: An approximate theory and accurate numerical experiments

# Water-wave gap solitons: An approximate theory and accurate numerical experiments

V. P. Ruban Landau Institute for Theoretical Physics, 2 Kosygin Street, 119334 Moscow, Russia
July 16, 2019
###### Abstract

It is demonstrated that a standard coupled-mode theory can successfully describe weakly-nonlinear gravity water waves in Bragg resonance with a periodic one-dimensional topography. Analytical solutions for gap solitons provided by this theory are in a reasonable agreement with accurate numerical simulations of exact equations of motion for ideal planar potential free-surface flows, even for strongly nonlinear waves. In numerical experiments, self-localized groups of nearly standing water waves can exist up to hundreds of wave periods. Generalizations of the model to the three-dimensional case are also derived.

###### pacs:
47.15.K-, 47.35.Bb, 47.35.Lf

## I Introduction

As we know from the nonlinear optics, specific self-localized waves can propagate in periodic nonlinear media, with a frequency inside a spectrum gap. These waves are referred to as gap solitons (alternatively called Bragg solitons; see, e.g., Refs.CM1987 (); AW1989 (); CJ1989 (); ESdSKS1996 (); PPLM1997 (); BPZ1998 (); RCT1998 (); CTA2000 (); IdS2000 (); CT2001 (); CMMNSW2008 ()). Also in the field of Bose-Einstein condensation, gap solitons (GS) have been known EC2003 (); PSK2004 (); MKTK2006 (). Recently, it has been realized that GS are also possible in water-wave systems R2008PRE-2 (). In particular, very accurate numerical experiments have shown that finite-amplitude standing waves over a periodic one-dimensional topography are subjected to a modulational instability which spontaneously produces Bragg quasisolitons — localized coherent structures existing for dozens of wave periods. However, in the cited work R2008PRE-2 (), no analytical approach was presented. As a result, many important questions about water-wave GS were not answered, concerning their shape and stability. The present work is intended to clarify this issue, at least partly. More specifically, for a given periodic bottom profile with a spatial period , we shall derive, in some approximation, coefficients for a standard model system of two coupled equations, describing evolution of the forward- and backward-propagating wave envelopes (see, e. g., Refs.AW1989 (); CJ1989 (); RCT1998 (); CT2001 ()),

 i(∂t±Vg∂x)A±=ΔA∓+(ΓS|A±|2+ΓX|A∓|2)A±, (1)

where is the time, is the horizontal coordinate in the flow plane, and are slow functions. Let at equilibrium the free surface be at . Then elevation of the surface is given by the following formula,

 η(x,t) = Re[A+eiκx−iω0(κ)t+A−e−iκx−iω0(κ)t] (2) + higher-order terms in κA±,

where is the wave number corresponding to the main Bragg resonance, is the frequency at the gap center, is the gravity acceleration, and is an effective depth of the water canal [definitely, is not a mean depth; more precisely it will be specified later by Eqs.(6) and (7)]. The coefficients in Eqs.(1) are: an effective group velocity , a half-width of the frequency gap, a nonlinear self-interaction , and a nonlinear cross-interaction .

Generally, it is assumed in derivation of the above simplified standard model that: (a) dissipative processes are negligible, (b) a periodic inhomogeneity is relatively weak (that is ), (c) the waves are weakly nonlinear, (d) original (without inhomogeneity) equations of motion, when written in terms of normal complex variables , contain nonlinearities starting from the order three:

 i˙Ak ≈ ω∗(k)Ak+12∫T(k,k2;k3,k4)A∗k2Ak3Ak4 (3) ×δ(k+k2−k3−k4)dk2dk3dk4,

where is a linear dispersion relation in the absence of periodic inhomogeneity (a weak inhomogeneity adds some small terms to the right hand side of Eq.(3); the most important effect arises from a term , where is a “small” linear non-diagonal operator). It is also required, (e) the coefficient of the four-wave nonlinear interaction should be a continuous function. In application to water waves the requirements (d) and (e) mean that: (i) all the second-order nonlinearities are assumed to be excluded by a suitable canonical transformation (the corresponding procedure is described, e.g., in Refs.K1994 (); Z1999 ()); (ii) the model system (1) can be good only in the limit of relatively deep water, since on a finite depth the function is known to contain discontinuities which disappear on the infinite depth (see, e.g., Ref.Z1999 ()). Therefore we introduce a small parameter

 ε≡exp(−2κh0)≪1, (4)

and we consider in the main approximation only the principal effect of weak spatial periodicity, namely creation of a narrow frequency gap with under the main Bragg resonance conditions. Then we imply a standard procedure for obtaining approximate equations for slow wave envelopes, where the deep-water limit of is used for the coefficients and . Thus we neglect in the actual nonlinear wave interaction some relatively small terms with coefficients of order .

Of course, the functions should be sufficiently “narrow” in the Fourier space, since dispersive terms proportional to second-order derivatives are not included into the model.

After derivation of all the coefficients in section II, some known “solitonic” solutions of Eqs.(1) will be compared to numerical results for exact hydrodynamic equations, with nearly the same initial conditions as in the solitons (in section III). We shall see that very long-lived self-localized groups of standing water waves are possible. In some region of soliton parameters, water-wave GS exist up to hundreds of wave periods, until unaccounted by Egs.(1) processes change them significantly. In section IV we discuss some promising directions of further research, concerning three-dimensional generalizations of the coupled mode equations. Some auxiliary calculations are placed in two Appendices.

## Ii Coefficients of the model

We start our consideration with a short discussion of conditions when dissipation due to bottom friction, caused by water (kinematic) viscosity , is not important in wave dynamics. Obviously, a viscous sub-layer should be relatively thin in this case: . In a nearly linear regime, a width of the sub-layer can be estimated as , where . This gives us the following necessary condition for applicability of the conservative theory:

 Λ3/4g1/4ν−1/2≫1. (5)

Generally speaking, one cannot exclude a possibility that in a strongly nonlinear regime the vorticity can sometimes be advected by a wave-produced alternating velocity field far away from the rigid bottom boundary. Such vortex structures are typically generated near curved parts of the bed, and they can significantly interact with surface waves. However, we assume this is not the case; otherwise, the problem becomes too complicated. Though we do not have simple criterion to evaluate influence of the bottom-produced vorticity, with m we still hope to be correct when neglecting water viscosity, as well as compressibility and surface tension. This allows us to exploit the model of purely potential free-surface ideal fluid flows, commonly used in the water wave theory.

Since in this work we consider the case of relatively deep water, we can write , where is the frequency corresponding to the infinite depth. Later we will see that values are of the most interest.

Let us introduce conformal curvilinear coordinates determined by an analytic function , with , so that

 x+iy=B(~ζ)=~ζ−2κ−1∞∑n=1βnεnsin(2nκ~ζ), (6)

with real coefficients . Without loss of generality, we assume . The unperturbed water surface corresponds to real values of , while at the bottom we have , and

 X(b)(ζ1)+iY(b)(ζ1)=B(ζ1−ih0) (7)

is a parametric representation of the bed profile, which can be highly undulating (see, for example, Fig.1).

In these conformal coordinates, a spectrum of linear potential waves is determined through the following equation (compare to Ref.R2004PRE (), where an analogous approach but slightly different notations were used):

 (8)

with

 B′(ζ1)=1−4∞∑n=1nβnεncos(2nκζ1). (9)

Here is a linear operator which is diagonal in Fourier representation: for any function we have

 [^ktanh(h0^k)]f(ζ1)=∫ktanh(h0k)fkeikζ1dk/2π. (10)

The eigenfunction takes the following form,

 Ψμ(ζ1)=eiμζ1+∞∑n=−∞cne2inκζ1, (11)

with some coefficients . With a given , we have an infinite homogeneous linear system of equations for . Non-trivial solutions exist for some discrete values . The first gap in the spectrum is the difference between the two first eigenvalues at , that is . Approximately, for small these eigenvalues are determined by the coefficient (compare to Ref.R2004PRE ()),

 ω2(1,2)(κ)≈gκtanh(h0κ)(1∓2εβ1). (12)

It should be noted that corresponds to , while corresponds to . Thus, in the first order on , the half-width of the gap in the spectrum of linear waves is

 Δ≈ω∗(κ)εβ1≡ω∗(κ)~Δ, (13)

where is a small dimensionless quantity.

As to the nonlinearity coefficients and , their values for the case of infinite depth (in other words, their zeroth-order approximations in ) can be easily extracted from Ref.OOS2006 ():

 ΓS≈12ω∗(κ)κ2,ΓX≈−ω∗(κ)κ2. (14)

It should be noted, the cited work OOS2006 () relies on results obtained earlier by Krasitskii K1994 (), who calculated kernels of so-called reduced integrodifferential equation for weakly nonlinear surface water waves (see also the paper by Zakharov Z1999 (), and references therein). It is important in many aspects that for deep-water waves the coefficients and have the opposite signs, and their ratio is .

With the same zeroth-order accuracy, the group velocity is

 Vg≈12ω∗(κ)κ. (15)

Now all the coefficients have been derived, and the simplified coupled-mode equations for relatively deep water waves in Bragg resonance with a periodic bottom take the following explicit form:

 i(∂tω∗+∂x2κ)a+=~Δa−+12(|a+|2−2|a−|2)a+, (16) i(∂tω∗−∂x2κ)a−=~Δa++12(|a−|2−2|a+|2)a−, (17)

where are dimensionless wave amplitudes. Analytical solutions are known for the above system (see AW1989 (); CJ1989 (); RCT1998 (); CT2001 ()), describing moving localized structures, the gap solitons. In the simplest case the velocity of GS is zero, and the solutions essentially depend on a parameter , a relative frequency inside the gap ():

 a± = √I(x)exp[−iδ~Δω∗t+iγ0±iφ(x)], (18) I(x) = 4~Δ(1−δ2)cosh[4~Δ√1−δ2κx]+δ, (19) φ(x) = arctan[√1−δ1+δtanh(2~Δ√1−δ2κx)]. (20)

These expressions correspond to purely standing, spatially localized waves with frequency (concerning their stability, see Ref.RCT1998 (), where, however, stability domains were presented for a different ratio ; there are some numerical indications that the above GS are stable in a parametric interval , where a critical value ).

It should be noted, one can hardly expect a detailed correspondence between the very simple model (16-17) and the fully nonlinear dynamics, but just a general accordance sometimes is possible. In particular, the model does not describe nonlinear processes resulting in generation of short waves which take the wave energy away from a soliton, thus influencing its dynamics. The model is also not generally good to study collisions between solitons, since wave amplitude can significantly increase in intermediate states.

## Iii Numerical experiments

In order co compare the above approximate analytical solutions to nearly exact numerical solutions, we chose the following function :

 B(~ζ)=~ζ+iD0κln(1+εCe2iκ~ζ1+εCe−2iκ~ζ), (21)

with real parameters , and . Hence, . For and both close to 1, Eq.(21) gives periodically arranged barriers (see, for example, Fig.1). The barriers are relatively thin as , and relatively high as . However, in numerical experiments with high-amplitude waves, a strong tendency was noticed towards formation of sharp wave crests over very thin barriers (say, when ), already after a few wave periods. With sharp crests, the conservative potential-flow-based model fails (it is also clear that tops of narrow barriers must generate strong vortex structures). Therefore we took in most of our computations in order to have a smooth surface for a longer time. Exact equations for ideal potential free-surface planar flows were simulated (their derivation can be found in Ref.R2004PRE (), some generalizations are made in Refs.R2005PLA (); R2008PRE ()). As in Ref.R2008PRE-2 (), we dealt with dimensionless variables corresponding to , . The dimensionless time is then related to the physical time by a factor . For instance, the period of linear deep-water waves with the length is . At , we set the horizontal free surface, while the initial distribution of the surface-value velocity potential was

 ψ0(ζ1)=2κ−3/2√I(ζ1)cos[κζ1+φ(ζ1)]≈ΨGS(ζ1), (22)

in accordance with approximate relation .

Many simulations with different parameters were performed, and a very good general agreement was found between numerical and analytical results in the weakly-nonlinear case, that is for small steepness . So, with , , , and (example I), some noticeable deviations from the purely-standing-wave regime were observed only after (see Figs.2-3). In a real-world experiment it could be several minutes with m.

What is interesting, even for larger , up to , GS can exist for dozens of wave periods. A numerical example for such a relatively high-amplitude water-wave GS is presented in Figs.4-5, where , , , and (example II). In this simulation, there were 45 oscillations before sharp crests formation (see Figs.4 and 6). As to a further evolution of such GS, only in a real-world experiment it will be possible to get reliable knowledge about it, since various dissipative processes come into play.

Concerning water-wave GS with negative , their behavior for was found stable, while for the dynamics was unstable, and partial disintegration of GS was observed after a few tens of wave periods (not shown). However, in some numerical experiments, the life time of GS was limited by the above mentioned process of sharp crest formation rather than by their own instability in frame of the model (16-17), at least with (not shown).

Finally, we would like to present an example of interaction of two GS (example III). The bed parameters are , , . Both solitons initially had and they were separated by a distance . At we set the horizontal free surface and

 ψ0(ζ1)≈ΨGS(ζ1−33Λ)+ΨGS(ζ1+33Λ). (23)

This numerical experiment also describes interaction of a single GS with a vertical wall at . Surface profiles for several time moments are shown in Figs.7-8. We see that in this example the interaction between GS is attractive. They collide and produce a highly nonlinear and short wave group near .

## Iv 3D generalizations and discussion

In this work, coefficients of the standard model (1) were derived for water-wave GS in the approximation of relatively deep water. The frequency gap in this case is small (of order ) despite strong bed undulations. It seems that a more general situation of intermediate depth cannot be described by this basic model, since an interaction of the main wave with a long-scale flow (“zeroth harmonics”) is then essential and should be included into equations. At the formal level, this corresponds to the mentioned discontinuities of the four-wave matrix element on a finite depth. Actually, in a finite-depth dynamics, three-wave interactions are more essential, and therefore they cannot be removed efficiently by a weakly-nonlinear transformation. This is the main difference between the present third-order theory and previously developed second-order theories (see, for example, Ref.HM1987 ()).

So far we considered purely two-dimensional flows, with the single horizontal coordinate . Let us now introduce two important generalizations for three-dimensional flows. Below we only derive equations, but their detailed analysis will be a subject of future work.

In the first case, the bottom topography is still one-dimensional, but we take into account weak variations of the wave field along the second horizontal coordinate , simply by adding dispersive terms, proportional to , to the coupled-mode system, as written below:

 (i∂tω∗±i∂x2κ+∂2q4κ2)a±=~Δa∓+12(|a±|2−2|a∓|2)a±. (24)

In this system, a near-band-edge approximation for the upper branch of the linear spectrum gives a 2D focusing nonlinear Schroedinger equation (NLSE). Thus, in a long-scale limit, the system (24) exhibits a tendency towards wave collapse which is known as a typical feature of 2D NLSE dynamics.

In the second case, a periodic bottom profile is essentially two-dimensional, and in the horizontal Fourier-plane there are several pairs of Bragg-resonant wave vectors. For simplicity, we present below equations for the case when has the symmetry of a square lattice, with equal periods in both horizontal directions and :

 χ=∑n1n2αn1n2[cos(2n1κx+2n2κq)+cos(2n1κx−2n2κq)], (25)

where coefficients possess the symmetry . Let us consider interaction of two wave pairs having slow complex amplitudes and , with the first pair corresponding to wave vectors , and with the second pair corresponding to . It is important that the absolute values are equal to each other: . Again we will assume . It is convenient to use new horizontal coordinates:

 x1=q+x√2,x2=q−x√2. (26)

Elevation of the free surface is then given by the formula

 ϰη = Re{e−iΩ0t[a+eiϰx1+a−e−iϰx1 (27) +b+eiϰx2+b−e−iϰx2]}+…,

where , with and , and the dots correspond to higher-order terms (again, we should note that generally ). Approximate equations of motion for the amplitudes have the following form:

 i(∂tΩ∗±∂x12ϰ)a±=ϵ1a∓+ϵ2[b++b−]+∂Hnl∂a∗±, (28) i(∂tΩ∗±∂x22ϰ)b±=ϵ1b∓+ϵ2[a++a−]+∂Hnl∂b∗±, (29)

where small constants and depend on a given bed profile, and mean the complex conjugate quantities, and the function corresponds to nonlinear interactions. Using an explicit expression from Ref.Z1999 () for the deep-water four-wave resonant interaction , we have

 Hnl=14{|a+|4+|a−|4+|b+|4+|b−|4} −|a+|2|a−|2−|b+|2|b−|2 +τ⊥{|a+|2|b+|2+|a+|2|b−|2+|a−|2|b+|2+|a−|2|b−|2} −34[a+a−b∗+b∗−+a∗+a∗−b+b−], (30)

where is a normalized value of the matrix element for two perpendicular wave vectors of equal length (see Fig.9). Since , we actually may neglect in the terms proportional to .

Unfortunately, it is hardly possible to find some analytical space-dependent solutions for the nonlinear system (28)-(29), but it can be investigated by approximate methods.

The parameters , , and can in principle be calculated from solution of a linearized problem for water waves over a periodic 2D bed. An exact linearized equation for a surface value of the velocity potential can be written in the following form:

 (ω2/g−[^ktanh(h^k)]−^N)Ψω(r)=0, (31)

where is a radius-vector in the horizontal plane, , while is a self-conjugate linear operator corresponding to a bottom inhomogeneity. However, in three dimensions there is no compact form for , valid with any bottom profile. At the moment, there are only approximate expressions , obtained by expansion (up to a finite order ) of a vertical velocity at the level in powers of (relatively small) bottom deviation from a constant level . The linear self-conjugate operators have a general structure

 ^Nj=[cosh(h^k)]−1^Sj[cosh(h^k)]−1, (32)

with

 ^S1=(∇χ∇), (33)
 ^S2=−(∇χ∇)[tanh(h^k)^k](∇χ∇), (34)
 ^S3=(∇χ∇)[tanh(h^k)^k](∇χ∇)[tanh(h^k)^k](∇χ∇) +[12(∇χ2∇)(∇χ∇)−16(∇χ3∇)∇2], (35)

and so on, where is the horizontal gradient (see Appendix A).

Now we are going to calculate , , and . Let us note that with the four independent eigenfunctions in Bragg resonance are: , , , and . Accordingly, we have for the eigenfrequencies

 ω2cc/g≈ϰtanh(hϰ)+⟨Ψcc^NΨcc⟩/⟨Ψ2cc⟩, (36)

(where mean the average value in the horizontal plane), and analogously for and . Let us introduce short notations for small quantities: and similarly for and . Then we have approximate equalities,

 ωcc ≈ Ω∗(1−ϵ+νcc), (37) ωss ≈ Ω∗(1−ϵ+νss), (38) ωcs=ωsc ≈ Ω∗(1−ϵ+νcs). (39)

These frequencies should be identified with the eigenfrequencies of the linear part of system (28)-(29), for space-independent solutions:

 ω(1,1,1,1) = Ω∗[(1−ϵ0)+ϵ1+2ϵ2]=ωcc, (40) ω(1,1,−1,−1) = Ω∗[(1−ϵ0)+ϵ1−2ϵ2]=ωss, (41) ω(1,−1,1,−1) = Ω∗[(1−ϵ0)−ϵ1]=ωcs, (42) ω(1,−1,−1,1) = Ω∗[(1−ϵ0)−ϵ1]=ωsc. (43)

As the result, we obtain the required formulas for the model parameters:

 ϵ0 ≈ ϵ−14(νcc+νss+2νcs), (44) ϵ1 ≈ 14(νcc+νss−2νcs), (45) ϵ2 ≈ 14(νcc−νss). (46)

With Eqs.(32)-(35), calculation of , , and is straightforward if the function contains a finite number of Fourier harmonics, for example

 χ=α1κ(cos2κx+cos2κq)+α2κcos2κxcos2κq. (47)

Moreover, since , it is possible to simplify the operators by writing there instead of . By doing so and taking into account only and , for the bottom profile (47) we obtain approximately

 ϵ2 ≈ −2ϵ√5α1α2, (48) ϵ1 ≈ ϵ√2α2, (49) ϵ0 ≈ ϵ[1+8√5α21+(1+12√5)α22]. (50)

Thus, the expansion is certainly useful for analysis of the case , but it can hardly be valid for a strongly undulating bed. It should be noted that a global representation of the velocity potential in the form (54) (see Appendix A) is questionable in the general case. Derivation of for arbitrary , assuming , is an interesting open problem.

It is worth noting that an explicit (though approximate) form of operator allows us to derive weakly nonlinear equations of motion for water waves over a nonuniform 2D bottom. For example, the Hamiltonian functional (it is the kinetic energy plus the potential energy ) up to the 4th order in terms of the canonically conjugate variables and is written below:

 H ≈ 12∫{ψ^Kψ+gη2+η[(∇ψ)2−(^Kψ)2]}d2r (51) +12∫[ψ^Kη^Kη^Kψ+η2(^Kψ)∇2ψ]d2r,

where (see Appendix B). It is interesting to note that the bottom inhomogeneity comes into the Hamiltonian through the definition of operator only. For , it coincides with the previously known fourth-order Hamiltonian for water waves on a uniform depth (see, e.g., Ref.Z1999 (), and references therein). It is also clear that coupled-mode system (28)-(29) corresponds to the case , when the difference is neglected in the third- and fourth-order parts of the Hamiltonian, but it is kept in the second-order part. The functional determines canonical equations of motion,

 ηt = δHδψ≈^Kψ−(∇η∇)ψ−^Kη^Kψ (52) + ^Kη^Kη^Kψ+12[^Kη2∇2ψ+∇2η2^Kψ], −ψt = δHδη≈gη+12[(∇ψ)2−(^Kψ)2] (53) + (^Kψ)^K(η^Kψ)+η(^Kψ)∇2ψ.

Numerical simulation of these cubically nonlinear equations, with , will be an important subject of future research.

Further analytical and computational work is also needed to investigate formation of vortex structures near the bottom boundary and to evaluate their influence on the free surface dynamics. In any case, the present results, based on the 2D purely potential theory, are deserving much attention. Moreover, the author hopes that in a future real-world experiment all the mentioned dissipative processes will not be able to destroy water-wave GS for a sufficiently long time. Instead, with vortices and breaking wave crests, the predicted phenomenon of standing self-localized water waves over a periodic bed will be found even more rich, interesting, and beautiful.

Acknowledgments. These investigations were supported by RFBR grant No. 06-01-00665, by RFBR-CNRS grant No. 07-01-92165, by the “Leading Scientific Schools of Russia” grant No. 4887.2008.2, and by the Program “Fundamental Problems of Nonlinear Dynamics” from the RAS Presidium.

## Appendix A. Expansion of operator ^N

The expansion of in powers of is easily obtained from the integral representation of the velocity potential

 Φ(r,y)=∫[ϕkcoshk(y+h)coshkh+fksinhkyk]eikrd2k(2π)2, (54)

where is the Fourier transform of the velocity potential at , and is the Fourier transform of an unknown function which should be determined by substitution of Eq.(54) into the bottom boundary condition

 [∂Φ/∂y−∇χ⋅∇Φ]∣∣y=−h+χ=0. (55)

The resulting integral equation can be represented as follows,

 −∇⋅∫fkikcosh[k(χ(r)−h)]k2exp(ikr)d2k(2π)2−∇⋅∫ϕkiksinh[kχ(r)]kcosh[kh]exp(ikr)d2k(2π)2=0. (56)

It can be formally solved for by expanding Eq.(56) in powers of and assuming . For instance, equation (56) with the third-order accuracy is written below:

 {1+(∇χ∇)[tanh(h^k)^k]−(∇χ22∇)}[cosh(h^k)]f=[(∇χ∇)−(∇χ36∇)∇2][cosh(h^k)]−1ϕ. (57)

As the result, be obtain an approximate solution , where the operators , , and are given by Eqs.(32)-(35). A linearized system describing the free-surface dynamics is

 −ψt=gη,ηt=[^ktanh(h^k)]ψ+f, (58)

where is a surface value of the velocity potential (in the linear approximation ). It gives us the equation (31) for eigenfunctions corresponding to some fixed frequency .

## Appendix B. Hamiltonian of water waves up to 5th order

An approximate Hamiltonian of water waves can be easily derived by writing the kinetic energy of potential three-dimensional motion of an ideal fluid in the following form,

 K = 12∫d2rη(r)∫−h+χ(r)[(∂Φ/∂y)2+(∇Φ)2]dy=12∫ψ[∂Φ/∂y−∇η⋅∇Φ]∣∣y=ηd2r (59) = 12∫∇ψ⋅∫[ϕkiksinh[k(η(r)+h)]kcosh[kh]+fkikcosh[kη(r)]k2]exp(ikr)d2k(2π)2d2r = 12∫[ψ^Kϕ+η∇ψ⋅∇ϕ+η22∇ψ⋅∇^Kϕ−η36∇ψ⋅∇(∇2ϕ)+…]d2r,

with subsequent substitution

 ϕ≈ψ−η^Kψ+η^Kη^Kψ+η22∇2ψ−η^K(η^Kη^Kψ+η22∇2ψ)−η22∇2η^Kψ+η36∇2^Kψ. (60)

The approximate equality (60) follows from an expansion of Eq.(54): Thus,

 K = 12∫{ψ^K[ψ−η^Kψ+η^Kη^Kψ+η22∇2ψ−η^K(η^Kη^Kψ+η22∇2ψ)−η22∇2η^Kψ+η36∇2^Kψ] (61) +η∇ψ⋅∇[ψ−η^Kψ+η^Kη^Kψ+η22∇2ψ]+η22∇ψ⋅∇^K(ψ−η^Kψ)−η36∇ψ⋅∇(∇2ψ)+…}d2r.

After simplifying, we obtain , where

 H(4) = 12∫[ψ^Kη^Kη^Kψ+η2(^Kψ)∇2ψ]d2r, (62) H(5) = 12∫[η36(^Kψ)∇2^Kψ−ψ^Kη^Kη^Kη^Kψ−η33(∇2ψ)2−η2(^Kη^Kψ)∇2ψ−η22(^Kψ)∇2(η^Kψ)]d2r. (63)

In the same manner, it is also possible to derive the Hamiltonian with a higher-order accuracy.

## References

• (1) W. Chen and D.L. Mills, Phys. Rev. Lett. 58, 160 (1987).
• (2) A. B. Aceves and S. Wabnitz, Phys. Lett. A 141, 37 (1989).
• (3) B. J. Eggleton et al., Phys. Rev. Lett. 76, 1627 (1996).
• (4) D.N. Christodoulides and R.I. Joseph, Phys. Rev. Lett. 62, 1746 (1989).
• (5) T. Peschel, U. Peschel, F. Lederer, and B. A. Malomed, Phys. Rev. E 55, 4730 (1997).
• (6) I.V. Barashenkov, D.E. Pelinovsky, and E.V. Zemlyanaya, Phys. Rev. Lett. 80, 5117 (1998).
• (7) A. de Rossi, C. Conti, and S. Trillo, Phys. Rev. Lett. 81, 85 (1998).
• (8) C. Conti, S. Trillo, and G. Assanto, Phys. Rev. Lett. 85, 2502 (2000).
• (9) T. Iizuka and C. Martijn de Sterke, Phys. Rev. E 61, 4491 (2000).
• (10) C. Conti and S. Trillo, Phys. Rev. E 64, 036617 (2001).
• (11) K.W. Chow et al., Phys. Rev. E 77, 026602 (2008).
• (12) N. Efremidis and D. N. Christodoulides, Phys. Rev. A 67, 063608 (2003).
• (13) D. E. Pelinovsky, A. A. Sukhorukov, and Yu. S. Kivshar, Phys. Rev. E 70, 036618 (2004).
• (14) M. Matuszewski et al., Phys. Rev. A 73, 063621 (2006).
• (15) V. P. Ruban, Phys. Rev. E 77, 055307(R) (2008).
• (16) V. P. Krasitskii, J. Fluid Mech. 272, 1 (1994).
• (17) V. Zakharov, Eur. J. Mech. B/Fluids 18, 327 (1999).
• (18) V. P. Ruban, Phys. Rev. E 70, 066302 (2004).
• (19) M. Onorato, A. R. Osborne, and M. Serio, Phys. Rev. Lett. 96, 014503 (2006).
• (20) V. P. Ruban, Phys. Lett. A 340, 194 (2005).
• (21) V. P. Ruban, Phys. Rev. E 77, 037302 (2008).
• (22) T. Hara and C. C. Mei, J. Fluid Mech. 178, 221 (1987).