## Abstract

A highly accurate technique for calculating the parameters of the Stokes 120° gravity wave on deep water is proposed. Tanaka's nonlinear transformation of the inverse plane is used to stretch the region near the wave crest for accelerating the convergence of the classical Michell series. Several families of new irregular self-intersecting profiles with a 120° angle at the crest are presented as well. They all degenerate into the Stokes 120° profile as numerical accuracy increases. These solutions seem to represent some kind of the so-called parasitic (or ghost) solutions, which emerge due to discretization of the original continuous problem. Probable physical relevance of such irregular solutions is discussed. The validity of the Stokes theorem about a 120° angle is tested numerically. No solutions with crest angles different from 120° were found. Hence, the Stokes 120° wave seems to be the unique gravity wave with a corner at the crest.

## 1. Introduction

The discovery of giant ocean waves, which are now called freak or rogue waves (e.g. the recent review by Kharif & Pelinovsky 2003), renewed the interest of scientific community in classical water-wave problems. Proper understanding of rogue wave phenomena requires a deeper theoretical and experimental background than that available up to date. Contemporary computational techniques allow the results obtained several decades ago to be revised and some new effects to be revealed. The canonical model of hydrodynamics, where fluid is assumed to be ideal and incompressible and waves are supposed to be potential, two-dimensional and symmetric, provides a good starting point for theoretical studies of steep wave phenomena. Having been engaged in the project on extreme ocean waves, Lukomsky *et al*. (2002) calculated the new approximate wave solutions different from the well-known Stokes wave solutions (see also Lukomsky & Gandzha 2003). These solutions describe irregular flows with stagnation points inside the flow domain and discontinuous streamlines near wave crests. The existence of such approximate solutions was confirmed by Debiane & Kharif (see Gandzha *et al*. 2002) and Clamond (2003). Discontinuous streamlines of irregular flows closely resemble the vertical jets with sharp-pointed tips from standing gravity waves forced beyond the maximum height (see Longuet-Higgins 2001; Longuet-Higgins & Dommermuth 2001). However, there have been some disputes as to whether such irregular flows/waves can really exist or they are only artefacts of numerical treatment (see Constantin 2004; Lukomsky *et al*. 2004). To shed some light on this problem, we decided to calculate directly the Stokes 120° wave rather than the almost-highest gravity waves considered in the studies cited above and to search for other sharp-crested solutions with a corner at the crest. Note in this respect that the uniqueness of the Stokes 120° solution has not been established rigorously yet, although the existence of the highest Stokes wave was proved by Amick *et al*. (1982) and Plotnikov (1982) (for details see Plotnikov & Toland 2004; Fraenkel 2007).

The theory of corner waves dates back to Stokes (1880) who conjectured that finite-amplitude gravity waves (which now bear his name) are limited by a sharp-crested configuration, i.e. the steepest possible oscillatory irrotational wave forms a sharp angle at its crest. Stokes proved that this angle should be equal to 120°. Michell (1893) was the first to compute such a wave numerically. Based on the fact that the wave must have a sharp 120° crest, he developed a series expansion of the complex velocity in terms of the complex potential, with the series coefficients determined by a set of nonlinear algebraic equations. With the first three coefficients taken into account, Michell estimated the steepness of the 120° wave at *A*≈0.142 (the trough-to-crest height divided by wavelength) and the ratio of the phase speeds of the 120° wave and infinitesimal wave at *c*≈1.2. Havelock (1918) took into consideration one more term of the Michell expansion and obtained *A*≈0.1418. Nekrasov (1919) considered the expansion inverse to that of Michell and acquired similar results. Upon the appearance of electronic computers, Yamada (1957) computed 12 coefficients of the series and obtained *A*≈0.1412.

Olfe & Rottman (1980) ultimately formalized the Michell method to include an arbitrary number *N* of coefficients. They carried out computations up to *N*=120 and obtained *A*≈0.14106 and *c*≈1.09228. Olfe & Rottman also calculated subharmonic gravity waves with 120° crests (such waves generally have several 120° crests and several rounded crests of smaller height per period). Srokosz (1981) supplemented the paper by Olfe & Rottman with the analysis of particle trajectories. Olfe & Rottman confirmed the earlier observation of Jeffreys (1951) that the rate of convergence of the Michell series drops as the number of terms increases. This is due to the fact that the second term of the local expansion about the corner includes an irrational power of the complex potential, as was shown by Grant (1973). In view of this, Williams (1981) developed a rather laborious technique for computing sharp-crested gravity waves with Grant's asymptotic taken into account explicitly. For the highest wave on deep water, Williams acquired the estimates *A*≈0.141063 and *c*≈1.09228. One can see, however, that Williams' technique did not significantly improve the accuracy of results as compared with the classical Michell method.

Hence, we decided to start with the standard Michell method and to use Tanaka's nonlinear transformation of the complex potential (Tanaka 1983) to improve the series convergence. We show that such a modification only slightly changes the resulting equations and allows the parameters of the Stokes 120° wave to be calculated with an accuracy of 11 decimal digits (within the standard double-precision arithmetic). While analysing the solutions produced by the Michell method, we came across new self-intersecting profiles with a 120° angle at the crest, some of them being very similar to the irregular solutions found by Lukomsky *et al*. (2002). We demonstrate that these new profiles all degenerate into the Stokes 120° profile as numerical accuracy increases. These solutions seem to be an example of the so-called parasitic (or ghost) solutions to nonlinear boundary-value problems (Domokos & Holmes 2003). Parasites are extraneous solutions which appear in discretizations of continuous problems. They are usually referred to as ‘non-physical’ or ‘numerically irrelevant’ solutions, since they disappear in the continuum limit. Ghosts were introduced by Domokos & Holmes (2003) as a new class of approximate solutions to the nonlinear boundary-value problems. They can be obtained as special limits of parasitic solutions and represent observable physical states in imperfect (with respect to the original model) experiments (see also Domokos 1997). In view of this, Domokos & Holmes suggested that the concept of non-physical or numerically irrelevant solutions, introduced for discretized boundary-value problems, should be re-examined. In other words, ghosts (and parasites as their discrete analogues) can reflect some physical effects that are not taken into consideration in the original idealized model. By perturbing the original boundary-value problem, one can gain access to such more specific solutions. We make a conjecture that weak surface tension can play the role of such a physical perturbation in our case.

Finally, the Michell method is generalized for an arbitrary crest angle with the aim of searching for solutions different from the Stokes 120° solution and to test the validity of the Stokes theorem about a 120° angle.

## 2. Governing equations

Consider steady two-dimensional periodic waves on the irrotational, inviscid and incompressible fluid with unknown free surface under the action of gravity. Waves are assumed to propagate without change of their form from left to right along the horizontal axis *x* with constant speed *c* relative to the motionless fluid at infinite depth. Choose a reference frame moving with the waves at speed *c*. In this wave-related reference frame with the horizontal coordinate *θ*=*x*−*ct* equal to the wave phase, the fluid motion is independent of time *t* and is governed by the equations(2.1)(2.2)(2.3)(2.4)where *ϕ*(*θ*, *y*) is the velocity potential (the velocity is equal to ∇*ϕ*), *η*(*θ*) is the elevation of the unknown free surface, *y* is the upward vertical axis and the dash over *η* designates averaging over the wave period. Here, equation (2.1) is the Laplace equation in the flow domain, equation (2.2) is the dynamical boundary condition (the Bernoulli equation, being the Bernoulli constant), equation (2.3) is the kinematic boundary condition (no fluid crosses the surface) and equation (2.4) is the condition that fluid is motionless at infinite depth. The dimensionless variables are chosen such that length and time are normalized by the wavenumber *k* and the frequency of an infinitesimal wave, respectively, *g* being the acceleration due to gravity. In this case, the dimensionless wavelength is *λ*=2*π*.

When the wave crest forms a corner, its summit becomes a stagnation point, where the fluid velocity is zero. In this case, one finds from equation (2.2) that(2.5)

By placing the *y*-axis origin at the wave crest, we get *η*(0)=0 and =0.

Equation (2.3) means that the free surface is a streamline, say *ψ*=0, where *ψ* is the stream function. Hence, it is convenient to proceed from the physical plane *ζ*=*θ*+i*y* to the inverse plane of the complex potential *w*=*ϕ*+i*ψ*, with the velocity potential and stream function being new independent variables (*θ*=0 corresponds to *ϕ*=0, and *θ*=*π* corresponds to *ϕ*=−*cπ*). The complex potential *w* and the complex coordinate *ζ* are analytic functions of each other, and any function *w*(*ζ*) satisfies the Laplace equation exactly due to the Cauchy–Riemann conditions(2.6a)

(2.6b)

Hence, the problem is to determine the function *w*(*ζ*), which is analytic in the lower part of the *w*-plane and satisfies the boundary conditions(2.7)(2.8)Equation (2.7) is obtained by differentiating the Bernoulli equation with respect to *ϕ*.

## 3. Method of solution

### (a) Michell's expansion

Let us transform the lower half of the *w*-plane into a unit circle with the use of the following transformation due to Nekrasov(3.1)

The free surface *ψ*=0 corresponds to |*u*|=1, the crest *ϕ*=0 to *u*=1, the troughs *ϕ*=±*π* to and the infinite depth (*ψ*=∞) is located at *u*=0.

Stokes (1880) showed that the asymptotic *w*∼*ζ*^{3/2} satisfies the Bernoulli equation in the infinitesimal vicinity of the corner summit. The corresponding asymptotic of the complex velocity is(3.2)

Hence, Michell (1893) proposed to search for solutions to equation (2.7) in the form of the following series(3.3)

The coefficients *b*_{n} must be real for the symmetric waves. To satisfy the condition at infinite depth, one should put *b*_{0}=−*c*. Other coefficients are determined by substituting expansion (3.3) into equation (2.7). Further procedure is briefly described by Olfe & Rottman (1980). They found that the rate of convergence of the method drops as the number *N* of the coefficients taken into consideration increases. Therefore, we decided to improve the convergence of the Michell series in order to get more accurate results.

### (b) Modification of the method

We introduce an additional transformation of the *u*-plane to a new *ξ*-plane defined as(3.4)In the *ξ*-plane, infinite depth is located at *ξ*=−*γ*, and *Χ* plays the role of the complex potential. This transformation, which was originally proposed by Tanaka (1983), maps a unit circle of the *u*-plane on to a unit circle of the *ξ*-plane with stretching the crest region while contracting the trough region (figure 1). The rate of stretch and contraction is intensified as the parameter *γ* approaches unity. As a result, the series constructed in the *ξ*-plane should converge much more rapidly than the same series in the *u*-plane. Note that the transformation (3.4) is only applicable to the waves with one crest and one trough per period and is not suitable for the subharmonic waves with several crests and troughs per period. As a matter of fact, the method can be generalized to deal with the subharmonic waves as well, but such a generalization lies beyond the scope of this study.

One can easily prove that asymptotic (3.2) remains unchanged in the *ξ*-plane(3.5)

Hence, the Michell expansion in the *ξ*-plane is written as(3.6)where *ν*≡1/3 and *N* is the number of coefficients taken into account in an approximate solution. The boundary condition at infinite depth (2.8) takes the form(3.7)

Further procedure is not much more awkward than that elaborated by Michell (1893) and Olfe & Rottman (1980). Substitution of expansion (3.6) into equation (2.7) and equating the coefficients of the like exponents to zero gives the following system of nonlinear algebraic equations for *b*_{n} and *c*(3.8)where

A step-by-step derivation of equations (3.8) is given in appendix A. Equations (3.8) supplemented with condition (3.7) were solved numerically by Newton's iterations. The Jacoby matrix of the system can be found symbolically and is presented in appendix B. The numerical procedure was implemented in Fortran 95. The equations were first solved in the simplest case *N*=1 and *γ*=0 with the use of the symbolic manipulator. These solutions were used as initial guesses in calculations with higher *N*. In the same way, the solutions obtained at some *N* were used as initial guesses in calculations at still higher *N*, and so on. As a rule, we simply doubled the number of equations at each new step. In this case, no more than five iterations were ordinarily required to achieve a tolerance of 10^{−14}. The run time of one iteration on AMD Athlon 1.4 GHz processor and the maximum run memory are listed for various *N* in the first column of table 1.

### (c) Physical quantities

Once the coefficients *b*_{n} and the phase speed *c* are found from equations (3.8), a variety of other wave characteristics can be evaluated. The wave steepness *A* is calculated as(3.9)

The wave profile *y*=*η*(*θ*) is found in parametric form from the following integrals(3.10a)(3.10b)where *Χ* ranges from 0 (crest) to *π* (trough) and(3.11)

The wave slope *η*_{θ}, convexity *η*_{θθ} and higher derivatives are calculated by numerical differentiation of the function *y*=*η*(*θ*).

The mean kinetic energy per unit horizontal area is computed as(3.12)where *I* is the mean momentum per unit area. The mean potential energy per unit area is calculated as(3.13)

The integrals were all evaluated numerically by the Simpson method. To resolve more accurately the singularity at the wave crest, an additional nonlinear transformation of the integration variable *Χ* was performed(3.14)

This transformation, which is due to Chen & Saffman (1980), acts in the same way as Tanaka's transformation (3.4). It stretches the crest region while contracting the trough region. The parameter *α* was chosen as close to unity as possible, namely, *α*=1–10^{−14}. This allowed us to achieve an integration accuracy of 10^{−14} in all our computations. The number of integration points ranged from 2000 for *N*=1000 to 16 000 for *N*=6400. The maximum run time of one integration with Fortran 95 and AMD Athlon 1.4 GHz processor did not exceed several seconds.

## 4. Results of computations

### (a) Stokes 120° wave

Table 1 lists some parameters of the Stokes 120° wave calculated at various *N* and *γ*. To test the overall accuracy of calculations, we substituted truncated series (A 3) and (A 4) into the Bernoulli equation (2.7) and computed the maximum deviation *Er* of the left-hand side of the equation from zero. The error decreases with increasing *N*, the magnitude of the last coefficient *b*_{N} decreasing. The case *γ*=0 corresponds to the standard Michell method. One can see that Tanaka's transformation improves the accuracy of the method by three orders of magnitude. However, note that the parameter *γ* cannot be chosen very close to unity, since the contraction of the trough region becomes so appreciable that numerical accuracy deteriorates. The greater the truncation number *N*, the larger the values of *γ* that can be chosen.

It can be seen that the results obtained by the modified Michell method at *N*=3200 and *γ*=0.99 are accurate to 11 decimal digits. Further improvement of the accuracy is possible only with the use of quadruple precision, since the coefficients *b*_{N} reach a double-precision accuracy limit of 10^{−14} at *N*=3200. Our estimates of the phase speed, steepness, potential energy and kinetic energy of the Stokes 120° wave coincide with the predictions made by Maklakov (2002) from computations of the limiting properties of almost-highest gravity waves (his values are also accurate to 11 decimal digits).

The behaviour of the Stokes 120° profile and its derivatives is of special importance in mathematical studies of the uniqueness of the Stokes wave solution (Plotnikov & Toland 2004). Hence, we calculated the slope *η*_{θ} and convexity *η*_{θθ} of the profile obtained at *N*=6400 and *γ*=0.95 (figures 2 and 3). One can see that at *θ*→0±, and the slope is almost linear on the rest of the wave period. The convexity *η*_{θθ} tends to zero at *θ*→0±, and it has a maximum at *θ*≈±0.597. The third and higher derivatives are infinite at *θ*→0±.

In view of the mathematical studies on Stokes waves, we also decided to test the validity of the Stokes theorem about a 120° angle and to check the behaviour of Michell's expansion in the case of other crest angles. The results of this investigation are presented in appendix C.

### (b) Irregular self-intersecting profiles

Consider algebraic equations (3.8) at *γ*=0 (standard Michell's method) in the simplest case *N*=1. The substitution allows a single equation to be obtained after elimination of *c*^{2}(4.1)

This equation has three real roots

The first solution yields the Stokes 120° wave. The second solution corresponds to the subharmonic wave of class 2, which was first computed by Olfe & Rottman (1980) (it has one 120° crest and one rounded crest per period). The third solution was obtained neither by Olfe & Rottman nor by Michell. Olfe & Rottman missed this solution since they derived, for some reason unknown to us, the fourth-degree equation instead of the fifth-degree equation (4.1), although their original set of equations was identical to our equations (3.8) at *γ*=0.

At *N*=2, the third solution is transformed into

By increasing *N* further and using the solutions at smaller *N* as initial guesses for solutions at higher *N*, we traced this new solution up to *N*=1000. It yields a self-intersecting profile whose crest attains a ‘heart-shaped’ form at *N*>5 (figure 4). The heart-shaped bubble at the crest shrinks in both dimensions as *N* increases, and it should turn into a single point at *N*→∞, since pure gravity waves cannot have multi-valued profiles (Spielvogel 1970). To calculate the profile more accurately, we used Tanaka's transformation and gradually increased *γ* up to *γ*=0.98. The resulting profile became much closer to the Stokes 120° profile than the profile at *γ*=0. In line with this observation, table 2 demonstrates that the parameters of the heart-shaped solution slowly converge to the Stokes wave parameters as the numerical accuracy increases. Hence, it becomes clear that the solution will wholly degenerate into the Stokes solution at *N*→∞.

The heart-shaped profile shown in figure 4 is only the first of a plethora of self-intersecting profiles produced by the Michell method. Another typical profile is the ‘two-loop’ profile demonstrated in figure 5. This two-loop solution also degenerates into the Stokes 120° solution as numerical accuracy increases (table 3). Two examples presented in figures 4 and 5 represent two general kinds of self-intersecting profiles we found, namely, one-bubble profiles with an inverted 120° angle (like the heart-shaped profile) and two-bubble profiles with an included 120° angle (like the two-loop profile). The profiles of the same kind differ by the shape of their bubbles. The greater the truncation number *N*, the larger number of different solutions with self-intersecting profiles is possible. These solutions all converge to the Stokes 120° solution as numerical accuracy is improved. Hence, the Stokes 120° wave seems to be the unique gravity wave with a corner at the crest (see also appendix C).

### (c) Discussion

The two-bubble solutions produced by the Michell method are very similar to the irregular wave solutions found by Lukomsky *et al*. (2002). The irregular wave profiles by Lukomsky *et al*. (2002) do not form loops since the Fourier expansions used to calculate those profiles do not admit multi-valued profiles. Irregular waves seem to approximate the two-bubble 120° profiles in the same way as almost-highest waves approximate the Stokes 120° wave. The irregular wave solutions are expected to converge to the two-bubble 120° solutions while these two-bubble solutions converge to the Stokes 120° solution at *N*→∞. This resolves the numerical nature of the irregular flows by Lukomsky *et al*. (2002). These flows should all degenerate into the Stokes 120° flow at *N*→∞. Therefore, there actually exists only one family of the gravity waves known as the Stokes wave family.

Nevertheless, the irregular solutions with multi-valued profiles may have physical relevance. This point of view is strongly supported by the recent study by Domokos & Holmes (2003). They investigated approximate extraneous solutions which appear in discretized formulations of nonlinear boundary-value problems. Such solutions are also called parasites, since they disappear in the continuum limit. In view of this, they were assumed to be non-physical or ‘mathematically irrelevant’ for a long time (see a wide list of references in Domokos & Holmes 2003). In fact, the irregular solutions with self-intersecting profiles presented in this study are examples of such parasitic solutions. The main message of the paper by Domokos & Holmes (2003) is that parasitic solutions can be physically meaningful (see also Domokos 1997).

As an example, Domokos & Holmes (2003) considered the idealized Euler buckling problem. The problem is reduced to the following second-order ordinary differential equation with boundary conditions(4.2)where ( ) ′ denotes d/d*s*, *s* is the arc length and *α* is the slope of the elastic rod. The parameter *λ* ranges from zero to infinity, and *λ*=∞ corresponds to the singular limit when the rod forms a symmetric loop. The discretized version of equation (4.2) yields a plethora of parasitic solutions, which strongly depend on discretization number *N* (Domokos & Holmes 1993). These parasites were shown to vanish and converge to regular solutions in the classical limit, when, firstly, *N*→∞ at finite *λ* and, secondly, *λ*→∞. However, Domokos & Holmes (2003) proved that at certain conditions, when the limits *N*→∞ and *λ*→∞ are taken simultaneously, parasites can converge to a new class of approximate solutions bounded away from regular solutions. Such new solutions were called ghosts. There is experimental evidence (Domokos 1997; Domokos & Holmes 2003) that ghosts represent observable physical states of the buckled rod. They can be observed in imperfect (with respect to the original idealized model) physical experiments, when friction, noise and other perturbations of a physical nature manifest themselves. Domokos (1997) explained the appearance of parasites by a high degree of symmetry (degeneracy, structural instability) of the Euler boundary-value problem. By perturbing this idealized problem, one can obtain more generic cases, which in turn, might prove to be vastly more complex than the original, classical problem. Discretization is a special kind of such a perturbation.

Although the concept of ghosts was strictly introduced only for ordinary differential equations, there are no doubts that such solutions also exist in the case of partial differential equations. We did not intend to construct ghosts in our problem, but the analogy with the buckling rod is evident. Indeed, consider the parameterwhich appears in the Nekrasov integral equation (Plotnikov & Toland 2004). The case *μ*=3 corresponds to infinitesimal waves and the case *μ*=∞ corresponds to the highest Stokes wave with a 120° corner at the crest (here *q*(0) is the water velocity at the crest summit). The limit *μ*→∞ corresponds to the limit *λ*→∞ in the buckling problem, and the loop-shaped profile of the rod is the analogue of the Stokes 120° wave. Our irregular self-intersecting profiles at *μ*=∞ are the analogues of the discrete approximations to ghost solutions constructed by Domokos and Holmes at *λ*→∞. Such solutions were revealed because the degree of symmetry of the canonical model was reduced by discretization, the finite truncation number *N* being the corresponding mathematical perturbation. Since multi-valued profiles are typical for steep capillary and gravity-capillary waves (e.g. Debiane & Kharif 1997), we expect that weak surface tension would be the best physical perturbation to yield access to these ghost solutions and to give them real physical sense. The mathematical perturbation at large *N* is expected to correspond to the physical perturbation produced by a weak surface tension (longer waves), and the case of small *N* seems to correspond to the case of strong surface tension (shorter waves). When the wave motion is governed by only one factor (gravity field), the Stokes 120° wave is the only solution observed (the highest degree of symmetry/degeneracy). The solutions with multi-valued profiles (ghosts) can manifest themselves when the Stokes wave degeneration is removed by some external perturbation, e.g. by surface tension. Nevertheless, the validation of this assumption requires a deeper investigation of the surface tension effect on water waves.

## 5. Conclusions

The rate of convergence of Michell's (1893) series expansion was increased by several orders of magnitude with the use of Tanaka's (1983) nonlinear transformation. This allowed us to compute the parameters of the Stokes 120° wave with an accuracy of 11 decimal digits. A variety of new corner wave solutions with self-intersecting profiles were found. Some of these solutions closely resemble the irregular wave solutions by Lukomsky *et al*. (2002). They all degenerate into the Stokes 120° wave solution as numerical accuracy increases. An analogy between these new approximate irregular solutions and ghost solutions introduced by Domokos & Holmes (2003) for nonlinear boundary-value problems was drawn, and the probable physical relevance of such solutions was discussed. The validity of the Stokes theorem about a 120° angle was tested numerically, and the uniqueness of the Stokes 120° wave was confirmed.

## Acknowledgments

We express thanks and appreciation to Prof. V. I. Shrira for his help which allowed us to attend the EGU General Assembly (2005) held in Vienna, where some preliminary results of this study were reported, and to Dr D. Clamond for his valuable collaboration.

## Footnotes

- Received December 26, 2006.
- Accepted February 27, 2007.

- © 2007 The Royal Society