1 Introduction
In his now-classic work, Stokes had provided a supremely elegant proof that if the crest of an irrotational wave has a sharp edge, then this edge must necessarily make an angle of
$120^{\circ }$
(Stokes Reference Stokes1880b
). The argument is as follows: for an ideal fluid in deep water with complex potential
$f={\it\phi}+\text{i}{\it\psi}$
and spatial variable
$z=x+\text{i}y$
, the constant-pressure condition on the free surface,
$y={\it\eta}(x)$
, may be taken as
$$\begin{eqnarray}\displaystyle \frac{1}{2}\left|\frac{\text{d}f}{\text{d}z}\right|^{2}+g{\it\eta}=0, & & \displaystyle\end{eqnarray}$$
where
$g$
is the gravitational parameter, the corner of the wave is assumed to lie at the origin with
$z=0$
and
${\it\phi}=0$
, and the free surface is the streamline
${\it\psi}=0$
. In the neighbourhood of the corner, if we assume the complex potential takes the form
$f\sim Az^{n}$
, then it necessarily follows from (1.1) that
$n=3/2$
. Examination of the angular change in
$z$
across the crest yields the desired result.
In fact, Stokes’ demonstration was somewhat unexpected, as it corrected an earlier claim by Rankine. As Stokes explains:
In a paper published in the Philosophical Magazine, Vol. XXIX (1865), p. 25, Rankine gave an investigation which led him to the conclusion that in the steepest possible oscillatory wave of the irrotational kind, the crests become at the vertex infinitely curved in such a manner that a section of the crest by the plane of motion presents two branches of a curve which meet at a right angle. (Stokes Reference Stokes1880b , p. 225)
After correcting Rankine’s result, Stokes then wondered whether such a wave of greatest height could be obtained in practice:
This however leaves untouched the question whether the disturbance can actually be pushed to the extent of yielding crests with sharp edges …After careful consideration I feel satisfied that there is no such earlier limit, but that we may actually approach as near as we please to the form in which the curvature becomes infinite, and the vertex becomes a multiple point where the two branches with which alone we are concerned enclose an angle of
$120^{\circ }$
. (Stokes Reference Stokes1880b
, p. 227)
Stokes’ conjecture on the existence of such waves would spark a century long inquiry, until it was positively affirmed by Toland (Reference Toland1978) with the help of many other foundational works (see the review by Toland (Reference Toland1996) for a more comprehensive listing). Equally important, but distinct from studies on the theoretical existence of such waves, is research that attempts to describe the structure of the wave solutions close to the wave of greatest height. We shall focus on this line of inquiry.
For any wave with amplitude below that of the greatest, the main crest singularity moves off the free surface and lies on the imaginary axis in the complex plane. It was Grant (Reference Grant1973) who first highlighted the complexity of this process, and primarily demonstrated two key difficulties. First, higher-order corrections to Stokes’ local result for the steepest wave are extremely non-trivial and require expansions in transcendental powers (cf. Trinh, Chapman & Vanden-Broeck Reference Trinh, Chapman and Vanden-Broeck2011, appendix A for a similar discussion). Second, for wave solutions away from the steepest configuration, the singularities change in power. As Grant writes:
Now also consider the effect of increasing the amplitude. At maximum amplitude,
$z(f)$
has one singularity, of order
$2/3$
. For any lesser amplitude, it has singularities only of order
$1/2$
. The only way a continuous approach to greatest amplitude is possible is for
$z$
to have several coalescing singularities. (Grant Reference Grant1973, p. 261)
This final statement, regarding the coalescence of square-root singularities, was frustratingly left unexplored. Although Grant understood (particularly aided by the numerical results of Schwartz (Reference Schwartz1972, Reference Schwartz1974)) that the main singularity would approach the free surface, there was no evidence given of any supplemental singularities, nor any explanation of the mechanism by which the local power would change.
The search for such singularities would not be straightforward. Previously, the location of the main singularity was identified using techniques in series acceleration, such as in the work of Schwartz (Reference Schwartz1974) and Longuet-Higgins & Fox (Reference Longuet-Higgins and Fox1978). As explained in Baker & Graves-Morris (Reference Baker and Graves-Morris1996, p. 44) branch points in an analytically continued function will typically manifest as sequences of poles in a Padé approximant. Thus, numerical results on the distribution of poles within a series-accelerated approximation can often be deciphered in order to intuit the presence of a branch-point singularity. Recently, Dyachenko, Lushnikov & Korotkevich (Reference Dyachenko, Lushnikov and Korotkevich2014, Reference Dyachenko, Lushnikov and Korotkevich2015), performed extensive numerical computations of the Stokes wave with great attention to the numerical error; these results were used to construct the Padé approximants, and then to reconstruct Grant’s main singularity.
Such series-acceleration techniques can be unreliable or difficult to interpret (see e.g. p. 20 of Drennan (Reference Drennan1988) and Dallaston & McCue (Reference Dallaston and McCue2010)). Indeed, as it pertains to the Stokes wave, following the early efforts of Schwartz (Reference Schwartz1974), Tanveer (Reference Tanveer1991) had described the situation as follows:
However Schwartz’s apparent conclusion on the form of the nearest singularity …is in error since it suggested that the branch-point power changes continuously between
$1/2$
and
$1/3$
. This erroneous conclusion highlights the difficulty of extracting reliable and accurate information from a Padé approximate method (which) is practically feasible only for the nearest singularity and just to leading order. (Tanveer Reference Tanveer1991, p. 140)
Using an alternative method of analytic continuation, Tanveer (Reference Tanveer1991) calculated precise estimates of the main singularity on the imaginary axis, and based on the form of the associated free-surface equation, further conjectured the existence of a mirror singularity, which lies on the negative imaginary axis in the adjacent unphysical Riemann sheet. In § 3, we will review Tanveer’s reason for believing the existence of the mirror singularity.
Although Tanveer’s work had confirmed the presence of the main singularity and conjectured the existence of its mirrored pair, the work did not comment on the elaborate Riemann sheet structure that would result from the combined influence of these two points. During a recent conference talk, Lushnikov, Dyachenko & Korotkevich (Reference Lushnikov, Dyachenko and Korotkevich2015), indicated that an infinite number of Riemann sheets would emerge if both singularities were continuously encircled in an alternating manner; these results have now appeared in a preprint (Lushnikov Reference Lushnikov2015).
The main goal of this paper is to construct the Riemann surface that represents the extension of the water wave into the complex plane; in this way, we shall develop a more complete picture of the singularity structure. Put very simply, we want to see the singularities. We compare and contrast various numerical methods for achieving this analytic continuation, and demonstrate that in fact, the singularity structure of a finite amplitude wave in deep water is much more complicated than previously anticipated. In particular, not only do we confirm the existence of Tanveer’s (Reference Tanveer1991) mirror singularity, but we also demonstrate the numerical existence of further square-root singularities, shifted in a diagonal fashion, away from the imaginary axis and lying on subsequent Riemann sheets of the water wave.
As we explain next, the appearance of such diagonally shifted singularities is often a generic phenomenon associated with nonlinear differential equations (cf. Chapman, Trinh & Witelski Reference Chapman, Trinh and Witelski2013). Moreover, singularities that arise on an unphysical Riemann sheet have been noted in other free-surface flows as well. We highlight, for example, the case of the wave–structure interaction studied in appendix A of Trinh & Chapman (Reference Trinh and Chapman2013b ), where such singularities were crucially related to the subsequent analysis of the resultant water waves.
1.1 The importance of understanding singularities
The characterization of the singularity structure of the Stokes wave, or indeed any solution of a singular differential equation, brings with it several advantages (Tanveer Reference Tanveer1991). First, knowledge of the set of singularities can lead to the development of better analytical methods, and in rare cases, even exact solutions. In the limiting configuration of a steep wave, the curvature near the crest is governed by the effects of the nearest singularities. Thus, the corresponding asymptotic analysis will often require detailed information regarding the type and distribution of such points. In § 10, we shall discuss the implications of our work with the previous matched asymptotics procedure of Longuet-Higgins & Fox (Reference Longuet-Higgins and Fox1977, Reference Longuet-Higgins and Fox1978) for the ‘almost-highest wave’.
Second, deeper understanding of the singularity structure can lead to better numerical algorithms for computing waves. Although Stokes’ classical approach is sufficient for calculating the profiles of moderate waves (as performed by e.g. Schwartz (Reference Schwartz1974) and Cokelet (Reference Cokelet1977)), the presence of the crest singularity severely limits the radius of convergence of Stokes’ series representation. In order to compute steep waves, it becomes essential to adopt a different series expansion that accounts for the change in the power of the main crest singularity (cf. Olfe & Rottman (Reference Olfe and Rottman1980), Vanden-Broeck (Reference Vanden-Broeck1986)). In theory, if we know the detailed locations and scalings of the set of singularities, we can subtract any undesired behaviour from the solution and extend the radius of convergence.
Finally, and perhaps of widest significance, the importance of knowing the steady-state singularities extends to problems of the time-dependent type as well. As explained by Tanveer (Reference Tanveer1991), the occurrence of finite-time interfacial singularities in many free-surface problems is typically associated with singularities in the unphysical complex plane hitting the real domain. For instance, most directly connected with the Stokes wave is the study of how the approach of the nearest crest singularity to the free surface is associated with the breaking of the wave (cf. Cokelet (Reference Cokelet1977), Schwartz & Fenton (Reference Schwartz and Fenton1982), Baker & Xie (Reference Baker and Xie2011)). However, the role played by complex singularities in the development of physical instabilities is a generic one, and there are a wide range of applications, ranging from Kelvin–Helmholtz and Rayleigh–Taylor instabilities, to rupture and pinch-off phenomena (see e.g. Tanveer Reference Tanveer1991 and Chapman et al. Reference Chapman, Trinh and Witelski2013).
2 Mathematical formulation
Let us consider symmetrical two-dimensional periodic waves of wavelength
${\it\lambda}$
under the influence of gravity
$g$
at the surface of a fluid of infinite depth. The waves move from right to left with constant speed
$c$
relative to an inertial frame, so we shall take a frame of reference moving with the waves and seek steady solutions. The fluid is assumed to be inviscid and incompressible, and the motion irrotational; thus the solution may be represented as an analytic function,
$z(f)$
, of the complex spatial variable
$z=x+\text{i}y$
and complex potential
$f={\it\phi}+\text{i}{\it\psi}$
, where
${\it\phi}$
is the velocity potential and
${\it\psi}$
the streamfunction.
The reader will note the standard reversal of the independent and dependent variables. Although Stokes (Reference Stokes1847) originally developed his method in the natural way, seeking to solve
$f$
as a function of
$z$
, he published a supplement several decades later (Stokes Reference Stokes1880a
) where he detailed the much simpler inverse formulation. The elegance of the trick relies on the fact that when solving for
$z(f)$
, the boundary conditions are imposed at a known location in the
$f$
-plane. Our mathematical presentation follows closely from those formulations discussed in Wehausen & Laitone (Reference Wehausen, Laitone, Flugge and Truesdell1960), Schwartz (Reference Schwartz1974), and Vanden-Broeck (Reference Vanden-Broeck1986, Reference Vanden-Broeck2010).
The Cartesian coordinates are chosen with the
$x$
-axis at the mean water level and the
$y$
-axis is directed vertically upwards against the direction of gravity. The free surface is chosen to be the constant streamline
${\it\psi}=0$
, and
${\it\phi}=0$
at the wave crest. The phase velocity,
$c$
, is defined to be the average horizontal velocity at a constant level,
$y$
, of the fluid. That is,
$$\begin{eqnarray}\displaystyle c=\frac{1}{{\it\lambda}}\int _{0}^{{\it\lambda}}\frac{\partial {\it\phi}}{\partial x}\,\text{d}x=\frac{1}{{\it\lambda}}{\it\phi}({\it\lambda},y), & & \displaystyle\end{eqnarray}$$
and consequently, it follows from the assumed symmetry and periodicity of the wave that
${\it\phi}(n{\it\lambda}/2,y)=(cn{\it\lambda})/2$
for integer
$n$
. Thus, we shall seek to solve for the fluid within the single periodic domain
$-(c{\it\lambda}/2)\leqslant {\it\phi}\leqslant (c{\it\lambda}/2)$
and
${\it\psi}\leqslant 0$
.
The system is non-dimensionalized using
${\it\lambda}$
as the unit length and
$c$
as the unit velocity; henceforth, all variables are dimensionless. The flow in the physical and potential planes is shown in figure 1(a,c). We have shifted the vertical axis of the physical plane for a better visualization of the singularities to be studied in later sections. We write Bernoulli’s equation in complex conjugate form (cf. Wehausen & Laitone Reference Wehausen, Laitone, Flugge and Truesdell1960, § 34) as
$$\begin{eqnarray}\displaystyle \frac{1}{2}+\frac{2{\rm\pi}}{{\it\mu}}\text{Im}(z)z^{\prime }({\it\phi})\overline{z^{\prime }({\it\phi})}=Bz^{\prime }({\it\phi})\overline{z^{\prime }({\it\phi})}, & & \displaystyle\end{eqnarray}$$
for
$-(1/2)\leqslant {\it\phi}\leqslant (1/2)$
and
${\it\psi}=0$
, and where the overbar denotes complex conjugation. The introduction of the Bernoulli constant,
$B$
, differs from (1.1) due to the chosen mean level, and we have defined
$$\begin{eqnarray}\displaystyle {\it\mu}=\frac{c^{2}}{g({\it\lambda}/2{\rm\pi})} & & \displaystyle\end{eqnarray}$$
for the square of the Froude number, representing the balance between inertial and gravitational effects.
(a) Profiles of numerically computed Stokes waves for
${\it\epsilon}$
ranging from 0.01 to 0.99 (increasing in steepness). The
${\it\epsilon}=0.99$
wave is shown dotted. Notice that rather than centring the mean height of all the profiles, the vertical axis has been shifted. The nodes are singularities in the analytically continued plane. (b) Values of
${\it\epsilon}$
versus
${\it\mu}$
. (c) Singularities shown in the complex
$f$
plane.
A solution of the nonlinear water wave problem then consists of solving (2.2) for an analytic function
$z=z(f)$
in
${\it\psi}\leqslant 0$
, such that
$z$
is periodic and the velocity is uniform at infinite depth. This requires that
$z(f+1)=z(f)$
and
$z^{\prime }(f)\rightarrow 1$
as
${\it\psi}\rightarrow -\infty$
.
One approach is to posit a series expansion for the solution. For example, following Stokes (Reference Stokes1880a ), we may set
$$\begin{eqnarray}\displaystyle z^{\prime }(f)=1+\mathop{\sum }_{n=1}^{\infty }a_{n}\text{e}^{-2{\rm\pi}\text{i}nf}, & & \displaystyle\end{eqnarray}$$
and then impose the boundary conditions. This determines the coefficients
$a_{n}$
, which are assumed to be real by the imposed symmetry of the wave.
The following numerical scheme is described by Vanden-Broeck (Reference Vanden-Broeck2010). First, the infinite series (2.4) is truncated after
$N-1$
terms. Next, the
$N+1$
unknowns
$B$
,
${\it\mu}$
and
$a_{n}$
,
$n=1,\ldots ,N-1$
, are determined by applying Bernoulli’s equation (2.2) at
$N$
equally spaced points in a single periodic domain. A final equation fixes the amplitude of the wave, and the resultant nonlinear system of
$N+1$
equations is solved using e.g. Newton’s method.
In regards to the amplitude condition, in this paper, we enforce a fixed value of
${\it\epsilon}$
defined by
$$\begin{eqnarray}\displaystyle {\it\epsilon}^{2}=1-\left|z^{\prime }(0)z^{\prime }({\textstyle \frac{1}{2}})\right|^{-2}, & & \displaystyle\end{eqnarray}$$
the same convention as in Vanden-Broeck (Reference Vanden-Broeck1986). In the limit
${\it\epsilon}\rightarrow 0$
, the wave height decreases to zero and
${\it\epsilon}=1$
corresponds to Stokes’ greatest wave. As noted by Vanden-Broeck, the series expansion of (2.4) provides an accurate approximation to the wave solutions for
${\it\epsilon}^{2}$
less than roughly
$0.6$
. We will mainly limit ourselves to waves that can be accurately calculated using (2.4).
It was shown by Chen & Saffman (Reference Chen and Saffman1980) that gravity waves of infinite depth are not unique, and the classical Stokes wave can bifurcate into new families of solutions at large values of
${\it\epsilon}$
. A similar result for finite depth was confirmed by Vanden-Broeck (Reference Vanden-Broeck1983). These irregular waves are distinguished by having more than one crest per wavelength. In this paper, we focus on the classic Stokes wave (and thus one crest per wavelength).
Numerical computations of the Stokes waves using (2.4) are presented in figure 1(a) for steepness values in the range
$0.01\leqslant {\it\epsilon}^{2}\leqslant 0.80$
. The table of
${\it\epsilon}$
versus
${\it\mu}$
values in figure 1(b) confirms that we have achieved at least five decimal accuracy in comparison with numbers reported in Cokelet (Reference Cokelet1977) and Vanden-Broeck (Reference Vanden-Broeck1986). The dashed curve in figure 1(a) corresponds to the wave with
${\it\epsilon}^{2}=0.99$
, which is computed using the generalized Havelock method detailed in Vanden-Broeck (Reference Vanden-Broeck1986). The crest singularities shown in the figure will be explained in the forthcoming sections.
3 Preliminary theory of the singularities
The case of
${\it\epsilon}=1$
corresponds to Stokes’ greatest wave and for this configuration, the singularity is located at the corner
$f=0$
of the
$120^{\circ }$
crest. For any
${\it\epsilon}<1$
, the singularity moves into the upper half-complex-plane, say, at
$f=f_{A}$
; we will continue to refer to this point as the main crest singularity. Let us first review Grant’s (Reference Grant1973) approach to studying
$f_{A}$
, and in particular, the argument that if such a singularity exists, then it must be of square-root type – here, by square-root type, we mean that the local power of a solution component (such as
$z$
or
$z^{\prime }$
) is
$m/2$
for some integer
$m$
.
3.1 On the crest singularity and its mirror image
Our first objective is to explain how the free surface, where
$f={\it\phi}\in \mathbb{R}$
, can be analytically continued into the upper half-plane. By the assumed symmetry of the free surface about the crest,
${\it\phi}=0$
, we have the following relation between
$z$
and its conjugate:
$$\begin{eqnarray}\displaystyle \overline{z({\it\phi})}=-z(-{\it\phi}), & & \displaystyle\end{eqnarray}$$
valid for real values of
${\it\phi}$
.
Next we apply the symmetry condition (3.1) to Bernoulli’s equation (2.2) and reduce the conjugate relations. The resultant equation can then be used as a prescription for continuing into the rest of the complex
${\it\phi}$
-plane. We write the free-surface potential as
${\it\phi}={\it\phi}_{r}+\text{i}{\it\phi}_{c}$
and, abusing notation somewhat, we relabel
${\it\phi}\mapsto f\in \mathbb{C}$
. After simplification, (2.2) yields
$$\begin{eqnarray}\displaystyle {\mathcal{Q}}(f){\mathcal{P}}(f)z^{\prime }(f)={\it\mu}, & & \displaystyle\end{eqnarray}$$
where we have introduced the functions
$$\begin{eqnarray}\displaystyle {\mathcal{P}}(f)\equiv {\it\mu}B+\text{i}{\rm\pi}[z(-f)+z(f)]\quad \text{and}\quad {\mathcal{Q}}(f)\equiv 2z^{\prime }(-f). & & \displaystyle\end{eqnarray}$$
When
$f$
is in the upper half-plane, then those components in
${\mathcal{P}}$
and
${\mathcal{Q}}$
that require values of
$z$
and
$z^{\prime }$
at
$-f$
are known from the convergent in-fluid series (2.4). Thus (3.2) forms a well-defined initial value problem for
$z$
in the upper half-plane. In § 4, we shall discuss how (3.2) can be used to explore the topology of the Riemann surface of
$z$
.
The following dominant balance argument can be verified a posteriori. Let us suppose that the singularity lies at
$f_{A}$
in the upper half-plane. Based on the imposed symmetry of the fluid and the form of (3.2), we would then expect
$f_{A}$
to be purely imaginary. Near this point, we assume that the solution can be represented in the form of
$$\begin{eqnarray}\displaystyle z\sim a_{0}+a_{1}(f-f_{A})^{{\it\alpha}}, & & \displaystyle\end{eqnarray}$$
for constants
$a_{0}$
and
$a_{1}$
, and
$0<{\it\alpha}<1$
. The restriction on
${\it\alpha}$
ensures that
$z^{\prime }$
is sufficiently singular. Since
${\mathcal{Q}}$
is an analytic function in the upper half-
$f$
-plane, then at leading-order in (3.2), we have
$$\begin{eqnarray}\displaystyle {\mathcal{P}}(f_{A})={\it\mu}B+{\rm\pi}\text{i}[z(-f_{A})+z(f_{A})]=0, & & \displaystyle\end{eqnarray}$$
which yields a value for
$a_{0}$
in (3.4). At next order, we may then use the fact that
${\mathcal{Q}}(f_{0})$
is bounded and non-zero to conclude that
$$\begin{eqnarray}\displaystyle (f-f_{A})^{2{\it\alpha}-1}=O(1)\Rightarrow {\it\alpha}=1/2, & & \displaystyle\end{eqnarray}$$
thus establishing the square-root behaviour as remarked by Grant (Reference Grant1973).
In fact, it is a consequence of the periodicity of the fluid that there will be further copies of the crest singularity, situated at
$f_{A}+\mathbb{Z}$
, with each copy separated by the unit length of the periodic interval. The general properties of the singularities that arise due to the underlying periodicity of (2.4) are non-trivial, and we explain these details in § 8. For the moment, let us focus on those singularities that are generated within the main interval.
Like Tanveer (Reference Tanveer1991, p. 149), we may develop the following argument to propose that a mirror singularity might be expected at
$-f_{A}$
. Suppose that using (3.2), the solution is analytically continued around the branch cut from
$f_{A}$
, and re-enters the lower half-plane; this part of the plane corresponds to an unphysical Riemann sheet, rather than the in-fluid domain. Since
${\mathcal{Q}}$
is singular at
$-f_{A}$
, then based on the condition (3.2), it is likely that
$z$
will also be singular. Although Tanveer did not elaborate beyond this point, we may go further and apply the same dominant balance arguments as above. It may be verified that, unlike the main crest singularity associated with property (3.5),
${\mathcal{P}}$
tends to a non-zero constant as
$f$
tends to
$-f_{A}$
. Using (3.6) for
${\mathcal{Q}}$
, an asymptotic balance for (3.2) yields
$$\begin{eqnarray}\displaystyle z^{\prime }=O(f+f_{A})^{1/2}\Rightarrow z\sim b_{0}+b_{1}(f+f_{A})^{3/2}, & & \displaystyle\end{eqnarray}$$
for constants
$b_{0}$
and
$b_{1}$
. The weaker nature of the singularity at the mirror point does not appear to have been remarked upon by others.
3.2 On general singularities
Later in § 6, not only do we confirm the existence of the main crest singularity and its mirror image, but in § 7, we will also demonstrate the existence of further singularities: a countably infinite set, diagonally shifted from the positions of
$f_{A}$
and
$-f_{A}$
. In fact, the dominant balance arguments that we have applied to conclude (3.6) and (3.7) can be repeated for a general singularity, say, at
$f=f^{\ast }$
.
We propose that there are two types of singularities, whose classification depends on the behaviour of
${\mathcal{P}}$
as
$f$
approaches
$f^{\ast }$
. The two types are described by
$$\begin{eqnarray}\displaystyle {\mathcal{P}}(f)\rightarrow \left\{\begin{array}{@{}ll@{}}0 & \text{with }z\sim a_{0}+a_{1}(f-f^{\ast })^{1/2},\\ \text{const.}\neq 0 & \text{with }z\sim b_{0}+b_{1}(f-f^{\ast })^{3/2}.\end{array}\right. & & \displaystyle\end{eqnarray}$$
Recall the argument applied to justify the existence of the mirror singularity; that is, a singularity at
$f_{A}$
must produce a non-analyticity in
${\mathcal{P}}$
or
${\mathcal{Q}}$
at the reflected point
$-f_{A}$
. This argument can be similarly applied to the case of a general singularity. Suppose for instance that
$f^{\ast }$
is a singularity of
$1/2$
type in (3.8), so with
${\mathcal{P}}\rightarrow 0$
. Studying the functional form of
${\mathcal{P}}$
in (3.3), it would seem unlikely for it, too, also to tend to zero as
$f$
tends to
$-f^{\ast }$
after crossing the lone branch cut from
$f^{\ast }$
; after all, there is no reason for the analytic continuation to be symmetric in this fashion. Therefore, it is likely that the mirror singularity
$-f^{\ast }$
switches to
$3/2$
type. Based on this very informal argument, we intuit that the chance existence of a singularity due to
${\mathcal{P}}=0$
can cause the emergence of a further mirror singularity where
${\mathcal{P}}\neq 0$
.
We will apply a variety of techniques to compute the locations of the singularities. The most crude approach is to observe the values of
$z$
and
$z^{\prime }$
along particular curves, surface meshes or contour plots. Branch cuts appear as differences in levels or kinks within such plots. Alternatively, we can observe numerical unboundedness as the singularity in
$z^{\prime }$
or
$1/z^{\prime }$
is approached. Generally, very approximate locations of the singularities are not difficult to obtain.
However, a more elegant approach can be designed based on contour integration. Near a general singularity, we have, after inversion,
$$\begin{eqnarray}\displaystyle z^{\prime }(f)\sim d(f-f^{\ast })^{{\it\alpha}}\Rightarrow f(z^{\prime })\sim f^{\ast }+(z^{\prime }/d)^{1/{\it\alpha}}, & & \displaystyle\end{eqnarray}$$
for constant
$d$
. We divide the second relation by
$z^{\prime }$
and integrate the result along a closed counterclockwise contour,
$C^{\prime }$
, containing
$z^{\prime }=0$
. The second term on the right-hand side of (3.9) integrates to zero, and this yields
$$\begin{eqnarray}\displaystyle f^{\ast }=\frac{1}{2{\rm\pi}\text{i}}\oint _{C^{\prime }}\frac{f}{z^{\prime }}\,\text{d}z^{\prime }=\pm \frac{1}{2{\rm\pi}\text{i}}\oint _{C}f\frac{z^{\prime \prime }}{z^{\prime }}\,\text{d}f. & & \displaystyle\end{eqnarray}$$
The last equality in (3.10) follows from the change of variables with
$\text{d}z^{\prime }=z^{\prime \prime }\text{d}f$
, and subsequent integration of the map of the contour,
$C$
, in the
$f$
-plane. The
$\pm$
sign follows from the ambiguity of not knowing whether a positive orientation in the
$z^{\prime }$
-plane is associated with a positive orientation in the
$f$
-plane. However, once this orientation has been established (e.g. for the main crest singularity), the conformal property of complex functions assures us that the sign will not change for other singularities. In our numerical computations to follow, we integrate (3.10) using the trapezoid rule, where the values of
$f$
and
$z^{\prime }$
are found from the scheme of § 4 and
$z^{\prime \prime }$
is calculated using centred differences.
Finally, note that based on the dominant balance arguments above, only singularities of square-root type are admissible in the analytic continuation of the Stokes wave. Locally, the Riemann surface near a square-root branch point contains only two sheets (or two branches), and thus positive or negative rotations are equivalent. However, due to the uniqueness of the analytic continuation process, this local property must also be true for the global surface. For instance, consider a point
$f$
on the Riemann surface for
$z(f)$
which is reached through a positive (counterclockwise) rotation about a given branch point. If we were to replace the positive rotation by the equivalent negative (clockwise) rotation, this process would yield the same value for
$z(f)$
.
4 Analytic continuation of the Stokes wave
Once the coefficients of the series representation have been computed using the free-surface condition (2.2), the in-fluid values, with
$\text{Im}(f)={\it\psi}<0$
are obtained directly from (2.4). However, the analytic continuation of the free surface to the upper half
$f$
-plane is severely limited by the main singularity identifiable with the crest of the steepest wave. In the past, authors had attempted to use techniques in series acceleration (primarily Padé approximants) in order to extend the radius of convergence of (2.4). However, such techniques can be unreliable and produce spurious singularities. We now implement two methods for constructing the analytic continuation of the Stokes wave.
4.1 Method 1: Reflection about the origin
The first method of analytic continuation, which we introduced in § 3.1, was originally proposed by Grant (Reference Grant1973) and then extended by Tanveer (Reference Tanveer1991). Here, we improve the method to allow for exploration onto further branches of the solution. From (3.2), we have the equation
$$\begin{eqnarray}\displaystyle z^{\prime }(f)=\frac{1}{2z^{\prime }(-f)\left\{B-\displaystyle \frac{{\rm\pi}}{{\it\mu}\text{i}}[z(f)+z(-f)]\right\}}. & & \displaystyle\end{eqnarray}$$
The simplest way to solve (4.1) is to begin on the known free surface calculated from (2.4), and integrate upwards along a vertical path in the
$f$
-plane. For example, consider the path labelled
${\mathcal{A}}{\mathcal{B}}$
in figure 2(a). The path is reflected about the origin to
$\underline{{\mathcal{A}}{\mathcal{B}}}$
. Since the values of
$z(-f)$
and
$z^{\prime }(-f)$
are known from the series expansion (2.4), then (4.1) provides an initial value problem that can be solved using any standard numerical stepping scheme (in our case, an adaptive Runge–Kutta solver).
Projection of the Riemann surface into
$(\text{Re}\,f,\text{Im}\,f,\text{Re}\,z)$
-space for the solution with
${\it\epsilon}=0.1$
. The surface is meshed by numerically integrating (4.1) along vertical paths beginning from the free-surface solution along the
$\text{Re}\,f$
axis.
We use this method to mesh the
${\it\epsilon}=0.1$
surface shown in figure 2(b). In the region
$\text{Im}(f)<0$
, this yields the physical fluid, which can be verified to match the convergent series values. In the upper half-plane, we observe a clear branch cut and branch point, labelled A, which is the main singularity as conjectured by Grant (Reference Grant1973). We shall say more about the specifics of the surface and singularity in § 6.
Now suppose we wish to obtain values on other branches of the solution. Consider the rectangular path labelled
${\mathcal{A}}{\mathcal{B}}{\mathcal{C}}{\mathcal{D}}$
in figure 2(a). The values of the solution along this rectangular path only depend on the values along
$\underline{{\mathcal{A}}{\mathcal{B}}{\mathcal{C}}{\mathcal{D}}}$
, and can thus be determined in an analogous manner, by solving (4.1) beginning from
${\mathcal{A}}$
.
Once the values in the upper half-plane are known, the path of continuation can proceed into the lower half-plane as well. If the prior journey had not circled a branch point, then we remain on the same Riemann sheet, and the lower half-plane will be the same as the in-fluid region calculated from (2.4). Indeed, this is the case in figure 2(b) where we observe the rectangular contour lying on the same surface mesh as determined through vertical integration. On the other hand, if a critical point was encircled, the continuation has taken us onto an unphysical Riemann sheet lying ‘adjacent’ to the physical fluid. This continuation process can then be iterated with the values in each lower or upper half-plane reflected in order to provide the values for the next half-plane.
The strength of the reflection scheme for continuation is that it is accurate; for a good approximation to the initial values on the free surface, the only error in the scheme occurs from the integration of (4.1), which can be controlled. However, the approach is limited by the requirement that the chosen path for continuation must be symmetric about the origin, and thus it does not allow us to derive the values
$z(f)$
for a general trajectory. In our numerical results for § 6, we implement the reflection method for rectangular contours centred about the origin.
4.2 Method 2: Boundary integral continuation
We propose an alternative approach that uses the boundary integral formulation of the Stokes wave problem and is not limited by symmetry. Boundary integral schemes provide an alternative treatment of potential free-surface problems, particularly in cases where the problem geometry (here fluid of infinite depth) is more complicated. An extensive review of such techniques can be found in the book by Vanden-Broeck (Reference Vanden-Broeck2010).
To begin, the complex velocity is written as
$\text{d}f/\text{d}z=\exp ({\it\tau}-\text{i}{\it\theta})$
for the fluid speed
$\text{e}^{{\it\tau}}$
and streamline angle by
${\it\theta}$
. The logarithmic hodograph variable,
${\it\Omega}(f)$
, is then introduced by
$$\begin{eqnarray}\displaystyle {\it\Omega}(f)\equiv {\it\tau}-\text{i}{\it\theta}=\log [1/z^{\prime }(f)]. & & \displaystyle\end{eqnarray}$$
Since
${\it\Omega}(f)$
is an analytic function in the fluid region, Cauchy’s integral formula applied to a rectangular contour bordering
${\it\phi}\in [-0.5,0.5]$
and
${\it\psi}\in (-\infty ,0]$
provides a boundary integral relationship between the real and complex parts of
${\it\Omega}$
. For values along the free surface,
$f={\it\phi}+\text{i}0$
, the boundary integral is given by
$$\begin{eqnarray}\displaystyle {\it\tau}({\it\phi}+\text{i}0)=\unicode[STIX]{x2A0D}_{-0.5}^{0.5}{\it\theta}_{free}({\it\varphi})\cot [{\rm\pi}({\it\varphi}-{\it\phi})]\,\text{d}{\it\varphi}, & & \displaystyle\end{eqnarray}$$
where the integral is of principal value.
At this point, the integral (4.3), combined with Bernoulli’s equation (2.2), provides a closed system for the determination of
${\it\tau}$
and
${\it\theta}$
along the axis, and indeed this forms the basis for such numerical boundary integral schemes. Instead, we will be interested in analytically continuing the integral (4.3) off the axis. The analytic continuation of the governing system gives two equations for
${\it\tau}={\it\tau}(f)$
and
${\it\theta}={\it\theta}(f)$
, valid now for
$f\in \mathbb{C}$
,
$$\begin{eqnarray}\displaystyle \text{e}^{3{\it\tau}}\frac{\text{d}{\it\tau}}{\text{d}f}+\frac{2{\rm\pi}}{{\it\mu}}\sin {\it\theta}=0\quad \text{and}\quad {\it\tau}+a\text{i}{\it\theta}={\mathcal{H}}(f), & & \displaystyle\end{eqnarray}$$
with
$a=\pm 1$
for
$\text{Im}(f)\gtrless 0$
, and where we have defined
${\mathcal{H}}$
according to
$$\begin{eqnarray}\displaystyle {\mathcal{H}}(f)\equiv \int _{-0.5}^{0.5}{\it\theta}_{free}({\it\varphi})\cot [{\rm\pi}({\it\varphi}-f)]\,\text{d}{\it\varphi}. & & \displaystyle\end{eqnarray}$$
The first equation in (4.4) is the differentiated and analytically continued form of Bernoulli’s equation (cf. equation (6) in Vanden-Broeck (Reference Vanden-Broeck1986)), written in terms of
${\it\tau}$
and
${\it\theta}$
. The second equation is the analytically continued boundary integral, where
$a=1$
for continuation into the upper half-plane and
$a=-1$
for the lower half-plane. Notice that upon taking
$\text{Im}(f)\rightarrow 0$
in the equation, we obtain a principal value integral and residue contribution, and subsequently return to (4.3), which indeed verifies that the continuation is done properly.
Consider a continuation into the upper half-plane that returns down to the real axis. Upon crossing into the lower half-plane, we have
$$\begin{eqnarray}\displaystyle {\it\tau}+\text{i}{\it\theta}={\mathcal{H}}(f)+2\text{i}{\it\theta}_{free}(f), & & \displaystyle\end{eqnarray}$$
where
${\it\theta}_{free}$
is the analytic continuation of the free-surface angle into the lower half-plane. Now if
${\it\theta}$
on the left-hand side of (4.6) is still on the same Riemann sheet as
${\it\theta}_{free}$
, that is to say, a branch point or pole was not encircled during the prior continuation, then
${\it\theta}(f)={\it\theta}_{free}(f)$
, and we recover (4.4b
) with
$a=-1$
. This would yield the values within the physical fluid region. However, if our analytic continuation took us onto a different Riemann sheet, then (4.6), with additional work to find the value of
${\it\theta}_{free}$
, is required in order to obtain the values in the lower half-plane.
The analogous rule applies to the case where we wish to continue from (4.6) back into the upper half-plane. For this, we must now include an additional residue,
$$\begin{eqnarray}\displaystyle {\it\tau}+\text{i}{\it\theta}={\mathcal{H}}(f)+2\text{i}{\it\theta}_{free}(f)-2\text{i}\hat{{\it\theta}}_{free}(f). & & \displaystyle\end{eqnarray}$$
If it was the case that no critical point was encircled while in the lower half-plane, we would remain on the same Riemann sheet as in (4.6),
${\it\theta}_{free}=\hat{{\it\theta}}_{free}$
, and we thus recover (4.4b
) with
$a=1$
. However, if we have traversed onto the another Riemann sheet, then in general
${\it\theta}_{free}$
will possess different values from
$\hat{{\it\theta}}_{free}$
, the latter of which is the analytic continuation of the free surface directly into the adjacent upper half-plane.
Let us now describe the numerical procedure for continuation along an arbitrary path in the
$f$
-plane, written as
${\it\gamma}$
and shown in figure 3. We first break
${\it\gamma}$
into segments each time the axis is crossed, and write
$$\begin{eqnarray}\displaystyle {\it\gamma}={\it\gamma}_{1}(s)\cup {\it\gamma}_{2}(s)\cup \cdots {\it\gamma}_{n}(s), & & \displaystyle\end{eqnarray}$$
with each subsegment parameterized by
$s\in [0,1]$
. The first contour,
${\it\gamma}_{1}$
is in the upper half-plane, and the initial point,
${\it\gamma}_{1}(0)$
lies on the free surface. By our notation,
${\it\gamma}_{i}$
lies in the upper half-plane if
$i$
is odd and the lower half-plane if
$i$
is even. Let
${\it\tau}_{i}$
,
${\it\theta}_{i}$
and
${\it\Omega}_{i}$
, be the analytically continued values along the subsegments,
${\it\gamma}_{i}$
. In addition to the subscript notation, we shall use an additional superscript, such as
${\it\Omega}_{i}^{(j)}$
, for individualized Riemann sheets, to be explained below.
Projection of the Riemann surface into
$(\text{Re}\,f,\text{Im}\,f,\text{Re}\,{\it\Omega})$
-space for the solution with
${\it\epsilon}=0.1$
. The surface is meshed by using the boundary integral continuation.
First, we solve for
${\it\Omega}_{1}={\it\Omega}_{1}^{(1)}$
in the upper half-plane using (4.4b
). Once
${\it\gamma}_{1}$
reaches its end along the real axis, we continue onto
${\it\gamma}_{2}$
in the lower half-plane and introduce the residue contribution (4.6). We call this newest contribution,
${\it\Omega}_{2}^{(2)}$
, that is to say, from (4.2) it is composed of
${\it\tau}_{2}^{(2)}-\text{i}{\it\theta}_{2}^{(2)}$
. The quantities
${\it\tau}_{2}^{(2)}$
and
${\it\theta}_{2}^{(2)}$
correspond to physical fluid quantities and can be calculated by returning to (4.4b
) and solving with
$a=-1$
; more simply, it consists of the in-fluid values directly obtained from the series representation (2.4).
Once we arrive at the end of
${\it\gamma}_{2}$
and continue onto
${\it\gamma}_{3}$
, then again, we must include a separate residue contribution and this involves the introduction of
${\it\Omega}_{3}^{(3)}$
. Thus, the desired values of
${\it\Omega}$
on each section of the plane will depend, in a recursive manner, on the values of the preceding Riemann sheets. Consequently,
$$\begin{eqnarray}\displaystyle {\it\Omega}_{n}^{(j)}\text{depends on }{\it\Omega}_{n}^{(n)},{\it\Omega}_{n}^{(n-1)},\ldots {\it\Omega}_{n}^{(j+1)}. & & \displaystyle\end{eqnarray}$$
Then the analytic continuation,
${\it\tau}_{n}^{(j)}$
and
${\it\theta}_{n}^{(j)}$
, along the curve
${\it\gamma}_{i}$
will involve solving the
$n$
th-order system of equations,
$$\begin{eqnarray}\displaystyle \frac{\text{d}{\it\tau}_{n}^{(j)}}{\text{d}f}=-\left[\frac{2{\rm\pi}}{{\it\mu}}\sin {\it\theta}_{i}^{(j)}\right]\text{e}^{-3{\it\tau}_{n}^{(j)}}\quad \text{for }1\leqslant j\leqslant n, & & \displaystyle\end{eqnarray}$$
with the initial conditions,
$$\begin{eqnarray}\displaystyle \left.{\it\tau}_{n}^{(j)}(f)\right|_{f={\it\gamma}_{n}(0)}=\left\{\begin{array}{@{}ll@{}}\!\displaystyle {\it\tau}_{free}({\it\phi})\Bigg|_{{\it\phi}={\it\gamma}_{n}(0)} & \text{for }j=n,\\ \!\displaystyle {\it\tau}_{(n-1)}^{(j)}(f)\Bigg|_{f={\it\gamma}_{n-1}(1)} & \text{for }1\leqslant j\leqslant n-1.\end{array}\right. & & \displaystyle\end{eqnarray}$$
The initial condition in (4.10b
) imposes that the newest residue contribution begins from the physical free surface,
${\it\tau}_{free}$
, calculated from (2.4) and (4.2). The initial conditions (4.10b
) impose that the previously generated contributions all continue from the end of the preceding path.
The system of differential equations is then closed with the system of recursively defined algebraic equations for
${\it\theta}$
given by
$$\begin{eqnarray}\displaystyle (-1)^{n+1}\text{i}{\it\theta}_{n}^{(n)}={\mathcal{H}}(f)-{\it\tau}_{n}^{(n)},\quad \text{for }n\geqslant 1, & & \displaystyle\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle (-1)^{j+1}\text{i}{\it\theta}_{n}^{(j)}={\mathcal{H}}(f)-{\it\tau}_{n}^{(j)}+2\text{i}\mathop{\sum }_{k=j+1}^{n}(-1)^{k}{\it\theta}_{n}^{(k)},\quad \text{for }1\leqslant j\leqslant n-1\text{ and }n\geqslant 2.\qquad & & \displaystyle\end{eqnarray}$$
Given a path,
${\it\gamma}={\it\gamma}_{1}\cup \cdots \cup {\it\gamma}_{n}$
, the system (4.10) can be solved using any standard ordinary differential equations scheme (in our case, an adaptive Runge–Kutta solver).
Projection of the Riemann surface into
$(\text{Re}\,f,\text{Im}\,f,\text{Re}\,{\it\Omega})$
-space for the solution with
${\it\epsilon}=0.1$
. The surface is meshed by using the boundary integral continuation.
The result of using the boundary integral continuation with the path shown in figure 3 is displayed in figure 4 for the Stokes wave with
${\it\epsilon}=0.1$
. The centreline of the surface mesh in figure 4 is the exact computed value, and for visualization purposes, we have extended the centreline into a surface mesh in the following way: at each point, we construct a planar patch using a vector,
$\boldsymbol{v}_{1}$
, along the arclength direction, and its normal,
$\boldsymbol{v}_{2}$
, assumed to satisfy conformality of the function
${\it\Omega}(f)$
. Unscaled, these are written
$$\begin{eqnarray}\displaystyle \boldsymbol{v}_{1}=[\text{Re}({\it\gamma}^{\prime }),\text{Im}({\it\gamma}^{\prime }),\text{Re}({\it\Omega}_{s})]\quad \text{and}\quad \boldsymbol{v}_{2}=[\text{Re}(\text{i}{\it\gamma}^{\prime }),\text{Im}(\text{i}{\it\gamma}^{\prime }),\text{Re}(\text{i}{\it\Omega}_{s})],\qquad & & \displaystyle\end{eqnarray}$$
where
${\it\Omega}_{s}$
is the derivative of
${\it\Omega}$
with respect to the parameterization variable
$s$
. The real operator in the third component of both vectors can be changed if a different three-dimensional representation is desired. In this way, the surface mesh shown in figure 4 is an accurate representation of the local topology of
${\it\Omega}$
.
Notice that in contrast to figure 2, in figure 4, we prefer to plot
${\it\Omega}$
since it is readily available from the calculations. However, the
$z$
values can be retrieved through (4.2). In this paper, we often switch between different solution measures and thus different visualizations of the Riemann surface. For example, we may prefer to plot the real or imaginary parts of the profile
$z$
, or the logarithmic hodograph variable,
${\it\Omega}=\log (1/z^{\prime })$
or the complex velocity,
$\exp ({\it\Omega})$
. Because the solution must always be of square-root type in
$z$
, then continuation properties between
$z$
and
$(1/z^{\prime })$
will be analogous. Depending on whether the real or imaginary part is taken, visualizations of
${\it\Omega}$
may display a logarithmic branch point. Generally, our strategy is to choose a measure where the relevant properties (e.g. branch cuts or singularities) are easily and unambiguously identified.
5 Notation for the singularities and Riemann sheets
In this section, we develop a notation to describe the singularity structure and collection of Riemann sheets for the complexified Stokes wave. For this, we must anticipate the results presented in §§ 6–8; there, we shall demonstrate that when viewed in the
$f$
-plane, a set of complex singularities of the Stokes wave will form two symmetrical V-shaped arrays, as illustrated in figure 5(a). Singularities are initially marked with an upper-case alphabetical letter, e.g.
$A$
,
$B$
,
$C$
and so forth. Those in the lower plane are distinguished by an underline, such as
$\underline{A}$
, while those in the left plane are distinguished by a superscripted
$\ell$
, as in
$B^{\ell }$
.
(a) Notation for the singularities; the sequence illustrated in (b) to (d) shows the continuation paths (arrow) to reach the Riemann sheet
$A\underline{A}\underline{B}$
(e). Branch cuts are taken upwards/downwards for the upper/lower half-planes.
In fact, the two main crest singularities,
$A$
and
$\underline{A}$
, will also appear in a periodic fashion, and such copies will be marked by a subscripted index, such as
$A_{\pm 1}$
,
$A_{\pm 2}$
and so forth, connected with the singularity
$A$
in the main interval. We have devoted a special section in § 8 to explain the periodic properties of the singularities; however, our principal challenge in this paper will be first to understand the singularities within the main periodic interval.
It is important to remark that singularities in the Stokes wave may not necessarily – and indeed they do not – exist on all Riemann sheets. The illustration in figure 5(a) is effectively a projection of all such singularities onto the two-dimensional
$f$
-plane. Thus, depending on the path of continuation that is used to reach a given point in the
$f$
-plane, the singularity may or may not exist at that given location.
In order to complete the description of the Riemann surface, we must thus describe how each individual Riemann sheet is reached using a path of analytic continuation. To this end, we introduce a sequence of letters for each individual sheet. The letters correspond to the singularities introduced above, with the additional symbol ‘
$\varnothing$
’ to refer to the main Riemann sheet on which the physical fluid is found, and from which the analytic continuation begins. Consider, for example,
$$\begin{eqnarray}\displaystyle \text{Riemann sheet }A\underline{A}\underline{B}. & & \displaystyle\end{eqnarray}$$
Each of the three letters in the sequence corresponds to crossing a branch cut from the associated singularity in the positive (counterclockwise) sense. Thus, Riemann sheet
$A\underline{A}\underline{B}$
is reached by beginning on
$\varnothing$
, and then encircling the singularities at
$A$
,
$\underline{A}$
and
$\underline{B}$
, one after the other, as illustrated in figure 5(b–e). As established in § 3.2, crossing a branch cut in the positive direction is equivalent to crossing in the negative direction; however, for precision of our numerical results, we will always consider the above notation as corresponding to a positive crossing.
6 The crest singularity and its mirror image
The purpose of this section is to apply the two methods for analytic continuation developed in § 4 in order to develop a visualization of the Riemann surface for the Stokes wave and to verify the presence of the proposed singularities. We will generally depend on the boundary integral scheme in § 4.2, which is unrestricted in the path of analytic continuation. However, the reflection scheme in § 4.1 was used throughout to verify the accuracy of paths that were sufficiently simple (notably paths in § 6.1).
6.1 The main Riemann sheet
$\varnothing$
and the crest singularity
$A$
Previously, during the explanation of the reflection scheme of § 4.1, we presented clear visual evidence of the main crest singularity, as it appears in the surface plot of figure 2. This singularity, which we refer to as ‘
$A$
’, lies on the main Riemann sheet
$\varnothing$
.
Let us confirm that
$A$
is indeed a square-root branch point. Using the boundary integral scheme, we analytically continue along a path that begins from the origin, tends upwards to
$f_{1}=0.05\text{i}$
, and then performs two rotations around an approximated singularity location for
$A$
at
$f=f_{A}$
. For the purpose of the visualization, we use a perturbed circle parameterized by
$$\begin{eqnarray}\displaystyle f({\it\theta})=|f_{1}-f_{A}|\left[1+{\it\delta}\sin \left\{{\textstyle \frac{1}{2}}({\it\theta}_{0}+{\it\theta})\right\}\right]\text{e}^{\text{i}({\it\theta}_{0}+{\it\theta})}, & & \displaystyle\end{eqnarray}$$
where
${\it\theta}_{0}=\text{Arg}(f_{1}-f_{A})$
and
$0\leqslant {\it\theta}\leqslant 4{\rm\pi}$
.
(a)
$\text{Re}(z)$
as a function of the arclength,
$s$
of the path of continuation; (b) the Riemann surface in
$(\text{Re}\,f,\text{Im}\,f,\text{Re}\,{\it\Omega})$
-space showing that two rotations around the
$\underline{A}$
singularity brings the path back to itself. The surface is computed using the boundary integral method and corresponds to
${\it\epsilon}=0.3$
Stokes wave. The singularity is located at
$f\approx 0.176\text{i}$
.
The numerical result, shown in figure 6, corresponds to the Stokes wave with steepness parameter
${\it\epsilon}=0.3$
. The three-dimensional representation in figure 6(b) shows the square-root behaviour as projected into
$(\text{Re}\,f,\text{Im}\,f,\text{Re}\,{\it\Omega})$
-space. In particular, note that the point,
$f_{1}=0.05\text{i}$
, which is shared by the two circular orbits, takes two distinct values at the beginning of the first and second orbits, before finally returning to its original value at the end. This behaviour is to be expected for a path that has traversed both Riemann sheets near a square-root branch point.
Once the numerical values along the path are known, we may then apply the integral equation (3.10) to calculate a more accurate location for the singularity. For example, for
${\it\epsilon}=0.3$
, this computation yields
$f_{A}\approx 0.1756\text{i}$
. As shown in figure 1(c), in the limit
${\it\epsilon}\rightarrow 0$
, the singularity tends to
$\text{i}\infty$
, while in the limit
${\it\epsilon}\rightarrow 1$
, the singularity tends to the origin. The general dependence of the singularity location on
${\it\epsilon}$
will be presented later in § 8.
Not only do these numerical results allow us to confirm the square-root nature of the solution, but we may also confirm the precise multiple of the square-root scaling. In § 3, we argued (from Grant Reference Grant1973), that
$$\begin{eqnarray}\displaystyle z\sim a_{0}+a_{1}(f-f_{A})^{1/2}, & & \displaystyle\end{eqnarray}$$
near the singularity. We now examine figure 6(a), which plots the real part of
$z$
as a function of the arclength along the path. As the path undergoes the two circular orbits,
$\text{Re}(z)$
undergoes one periodic cycle. This confirms the
$1/2$
power in (6.2), as opposed to some other
$(f-f_{A})^{m/2}$
power. This same scaling will not be true of the mirror singularity in the lower half-plane.
6.2 The adjacent sheet and its mirror singularity
$\underline{A}$
In § 3.1, we explained why a singularity might be expected at the mirror point
$-f_{A}$
on the Riemann sheet adjacent to
$\varnothing$
. The conjecture, which was first proposed by Tanveer (Reference Tanveer1991), is based on the reflection property of (3.2) and the non-analyticity of
${\mathcal{P}}$
and
${\mathcal{Q}}$
due to
$f_{A}$
. Let us call this mirror singularity ‘
$\underline{A}$
’, with the underline to remind us of its placement in the lower half-plane.
(a)
$\text{Re}(z)$
as a function of the arclength,
$s$
of the path of continuation. (b) The Riemann surface in
$(\text{Re}\,f,\text{Im}\,f,\text{Re}\,{\it\Omega})$
-space showing that two rotations around the
$\underline{A}$
singularity brings the path back to itself. The surface is computed using the boundary integral method and corresponds to the
${\it\epsilon}=0.3$
Stokes wave. The two singularities are located at
$f\approx \pm 0.176\text{i}$
.
Consider the case of
${\it\epsilon}=0.3$
. Using the boundary integral scheme, we perform the following two motions: first, beginning from
$f=0$
on the main sheet
$\varnothing$
, we perform a positive (counterclockwise) orbit around
$A$
, and then return to the origin. This brings us to the adjacent Riemann sheet. Next, we analytically continue downwards to
$f_{1}=-0.05\text{i}$
, and then perform two positive rotations centred on
$-f_{A}$
using the analogous parameterization to (6.1). The result is shown in figure 7(b) and confirms that two rotations brings the function back to its original value. Calculation of the integral (3.10) confirms that the singularity is indeed located at the reflected point.
In § 3.1, we demonstrated that if such the mirror singularity was found, then the solution must of the form
$$\begin{eqnarray}\displaystyle z\sim b_{0}+b_{1}(f+f_{A})^{3/2}, & & \displaystyle\end{eqnarray}$$
as opposed to the
$1/2$
scaling of (6.2). This is confirmed with figure 7(a) where
$\text{Re}[z(f)]$
is observed to undergo three nearly periodic cycles instead of the single cycle in figure 6(a).
One might ask: what happens if the path does not encircle a singularity? For this, a comparison can be drawn to the analytic continuation shown in figure 2 where the rectangular path (shown solid) does not encircle the
$A$
point, and therefore the continuation into the lower half-plane returns to the physical fluid (where it is analytic everywhere and single valued). Had there been no singularity in the lower half plane, we would expect all paths of continuation to be single valued along closed orbits.
7 The existence of diagonal singularities
By now, we have confirmed the existence of the two main singularities. Now we will show that further singularities emerge when proceeding deeper into the surface. These new singularities are accessed through paths that traverse branch cuts from
$A$
and
$\underline{A}$
, and are positioned diagonally from
$\underline{A}$
– more specifically, they tend to a constant argument as
${\it\epsilon}\rightarrow 0$
. As explained in § 3.2, the initial emergence of these diagonal singularities is somewhat out of chance; their existence depends on initial conditions of the Stokes wave and they are associated with locations where
${\mathcal{P}}=0$
. However, once they appear, they will also cause the appearance of further mirror singularities, where
${\mathcal{P}}\neq 0$
.
In fact, after submission of this work, we became aware that within a preprint, Lushnikov (Reference Lushnikov2015) had independently conjectured the existence of such diagonal singularities. The conjecture was based on an observation that the sudden transition of powers (cf. p. 257) from two singularities on the imaginary axis to the single highest-wave singularity could be achieved using an infinite sequence of nested square roots; the resultant function then admits off-axis singularities.
7.1 Diagonal singularities
$\underline{B}$
and
$\underline{B}^{\ell }$
As an example, consider the case of
${\it\epsilon}=0.5$
and the Riemann sheet
$A\underline{A}$
. In order to reach this sheet, we have taken the continuation path shown in the lower two-dimensional plane of figure 8(b); it begins on the positive
$\text{Re}(f)$
axis and forms an elliptical orbit around
$A$
and
$\underline{A}$
. We have chosen to illustrate the Riemann surface using its projection into
$(\text{Re}\,f,\text{Im}\,f,\text{Im}\,1/z^{\prime })$
-space, as this provides the most clear visualization. Unlike the previous visualizations, where we have only meshed along the continuation path, here, we provide a more complete surface mesh for the Riemann sheet. The surface mesh is created by first encircling
$A$
and
$\underline{A}$
, and then returning to the real
$f$
-axis. From the axis, the continuation moves in straight vertical lines (similar to the scheme to create figure 2). This process has the effect of imposing vertical branch cuts from the singularities.
These two pictures confirm the existence of the diagonal singularities
$\{\underline{B},\underline{B}^{\ell }\}$
for the Stokes wave with
${\it\epsilon}=0.5$
. Continuation follows a path that encircles
$A$
and
$\underline{A}$
, and then the purported singularity at
$\underline{B}$
. (a) The profile of
$\text{Re}(z)$
along the path confirms the singularity scaling. (b) The surface mesh for
$A\underline{A}$
and continuation mesh showing the double rotation around
$\underline{B}$
.
On the Riemann sheet shown in figure 8(b), we observe the usual
$A$
singularities and its periodic copies. Only
$A_{-1}$
and
$A_{1}$
are seen, but other copies are shifted a unit length away from one another. Notice that periodic copies of
$\underline{A}$
do not appear; we shall explain the periodicity structure in more detail in § 8. Most importantly, observe the presence of the two branch cuts clearly visible in the lower half-plane. These branch cuts emerge from two points located symmetrically about the imaginary axis:
$$\begin{eqnarray}\displaystyle f_{\underline{B}}\approx 0.557-0.396\text{i}\quad \text{and}\quad f_{\underline{B}^{\ell }}\approx -0.557-0.396\text{i}. & & \displaystyle\end{eqnarray}$$
The doubly rotated path in figure 8(b) confirms the square-root nature of the new diagonal singularity,
$\underline{B}$
, and examination of figure 8(a) verifies that the singularity is of
$1/2$
type, as anticipated by (3.8).
A plot of the locations of the
$B$
-type singularities is shown in figure 9(a), along with the dependence of the singularity magnitudes and arguments as a function of
${\it\epsilon}$
, seen in subplots (b,c). It is confirmed that as
${\it\epsilon}\rightarrow 0$
,
$|f_{B}|\rightarrow \infty$
, while
$\text{Arg}\,f_{B}$
tends to a constant. The numerical computations are restricted to
${\it\epsilon}\leqslant 0.6$
, since above this value it becomes difficult to preserve accuracy of the singular integral. Based on figure 9, however, we would expect that
$|f_{B}|\rightarrow 0$
as
${\it\epsilon}\rightarrow 1$
, but the limiting angle is unclear.
Location of the singularities
$\{B,B^{\ell },\underline{B}^{\ell },\underline{B}\}$
(from first to fourth quadrants) calculated using the integral (3.10). Singularities
$\underline{B}$
and
$\underline{B}^{\ell }$
are found on Riemann sheet
$A\underline{A}$
. Singularity
$B^{\ell }$
is found on
$A\underline{A}\underline{B}$
and singularity
$B$
is found on
$A\underline{A}\underline{B}^{\ell }$
.
7.2 Mirror diagonal singularities
$B$
and
$B^{\ell }$
We may confirm that as
$\underline{B}$
and
$\underline{B}^{\ell }$
are approached,
${\mathcal{P}}\rightarrow 0$
, which follows our argument leading to (3.8). Thus, from the same argument, it is expected that the singularities (7.1) will cause the emergence of mirror singularities (reflected about the origin) due to the non-analyticity of
${\mathcal{P}}$
and
${\mathcal{Q}}$
. In fact, these mirror singularities can be more simply argued based on the reflection scheme of § 4.1.
Consider a rectangular path, similar to figure 2, that orbits both
$A$
and
$\underline{A}$
. In the first quadrant (top right) of sheet
$A\underline{A}A$
,
$z(f)$
is derived based on values of
$z(-f)$
from the third quadrant (bottom left). Consequently,
$z(-f)$
lies on Riemann sheet
$A\underline{A}$
, which contains the
$\underline{B}^{\ell }$
singularity with
${\mathcal{P}}\rightarrow 0$
. Thus, by (4.1), it follows that there is a singularity at the reflected point in the first quadrant. These mirror singularities,
$B$
and
$B^{\ell }$
are confirmed in the later figure 10(f), and we may also verify that they are of the
$3/2$
type in (3.8).
Riemann sheets for the
${\it\epsilon}=0.5$
Stokes wave. The white nodes mark the crest singularity and its two periodic copies, the mirror singularity and the first diagonal singularities. Black nodes mark singularities arising from diagonals. The contours correspond to
$\text{Re}\,{\it\Omega}=\text{Re}({\it\tau}-\text{i}{\it\theta})$
.
7.3 Further diagonal singularities
The diagonal structure does not end there. In fact, it can be verified that there are two further singularities of
$C$
-type:
$\{C,C^{\ell },\underline{C},\underline{C}^{\ell }\}$
. For example, if we perform two elliptical orbits instead of the single one of figure 8(b), we would arrive at Riemann sheet
$A\underline{A}A\underline{A}$
. On this sheet, we discover a similar topology to the one for
$A\underline{A}$
, but with the diagonal set
$\{\underline{B},\underline{B}^{\ell },\underline{C},\underline{C}^{\ell }\}$
. These
$C$
-type singularities are diagonally shifted further than for the previous case of the
$B$
-type singularities, and they too will induce the mirror set
$\{C,C^{\ell }\}$
.
If we strictly confine ourselves to climbing the
$A\underline{A}$
structure, then the pattern of singularities is clear: each additional tour adds a further set of diagonal points. Thus for example, the three-tour Riemann sheet
$A\underline{A}A\underline{A}A\underline{A}$
contains
$\{\underline{B},\underline{B}^{\ell },\underline{C},\underline{C}^{\ell },\underline{D},\underline{D}^{\ell }\}$
. The more difficult challenge is to characterize the surface once we proceed onto Riemann sheets accessed through the diagonal points. This will chiefly be our task in § 8.
8 The global structure of the Riemann surface
In this section, our goal is to provide a basic list of rules and observations that characterize the global structure of the Riemann surface for the Stokes wave. As it turns out, a full description of the surface will prove to be quite difficult, but we hope to give the reader a taste of the daunting complexity of the singularity structure.
Earlier in figure 5(a), we had presented our notation describing the main crest singularity
$A$
, its mirror image
$\underline{A}$
and also the main diagonal singularities that form the two
$V$
-shaped arrangements,
${\mathcal{V}}=\{B,B^{\ell },\underline{B},\underline{B}^{\ell },C,C^{\ell },\underline{C},\underline{C}^{\ell },\ldots \}$
. Our numerical results have indicated that beyond these ‘first diagonals’, further singularities (diagonals of diagonals) arise once the path of continuation encircles the former set. Newer diagonals do not share the same symmetry properties as
${\mathcal{V}}$
.
8.1 Symmetry-breaking paths of continuation
Our discussion will centre on numerical results for the Stokes wave with
${\it\epsilon}=0.5$
. In figure 10(a), we plot the locations of the singularities discussed earlier in §§ 6 and 7 using a white node. Contour plots for
$\text{Re}\,{\it\Omega}=\text{Re}({\it\tau}-\text{i}{\it\theta})$
are shown in subfigures (b,j). We have chosen this solution measure as it is one for which the singularities are apparent from the shape of the contours. The white nodes in the smaller subfigures correspond exactly to those in figure 10(a), and the black nodes mark additional singularities that arise once symmetry is broken in the path of analytic continuation.
This notion of symmetry breaking is a central reason for the difficulty in establishing the full Riemann surface. Earlier we presented a method of analytic continuation that was based on formula (4.1). Here, the reflective property of the equation indicates that if there exists a singularity at
$f=f^{\ast }$
, then there is likely a singularity at
$-f^{\ast }$
as well. However, this criterion only applies if, by the time the continuation path reaches
$-f^{\ast }$
, the former
$f^{\ast }$
is still singular on the same Riemann sheet. This is certainly true for the transition between Riemann sheets
$\varnothing$
and
$A$
, where singularity
$A$
is located on the imaginary axis on both sheets, and hence induces the creation of the mirror singularity
$\underline{A}$
. However, in general, we must expect that paths of continuation that encircle a diagonal singularity to break the symmetry of the Riemann sheet; thus, arguments that depend on reflection will no longer directly apply.
8.2 Rules and observations
By now, the first two observations are well known:
-
(i) On the main sheet
$\varnothing$
, there is only the single crest singularity
$A$
and its periodic copies,
$A_{\pm i}$
,
$i=1,2,3,\ldots$
-
(ii) Two rotations about any singularity returns the path to its original value; clockwise and counterclockwise rotations are equivalent.
The first observation is verified in figure 10(b), and in numerous points throughout this work. The second observation follows from the fact that, as discussed in § 3, the only possible dominant balance in (3.2) must necessarily yield singularities of square-root type. We have also made mention of the periodic copies of
$A$
, of which
$A_{1}$
and
$A_{-1}$
can be seen in the figure. This leads us to:
-
(iii) Periodic copies of
$A$
and
$\underline{A}$
can appear on certain Riemann sheets; when this occurs, they appear at the points
$\pm f_{A}+\mathbb{Z}$
.
Note that while the physical water wave is required to be periodic, there is no such limitation on the analytic continuation. However it is a consequence of the two methods of continuation explained in § 4 that certain periodic properties may still be preserved. As shown in figure 10(b), on
$\varnothing$
, there will be periodic copies of the main singularity,
$A$
, separated by unit lengths. In fact, this can be more easily seen by using the circular map given by
$t=\exp (-2{\rm\pi}\text{i}f)$
. This maps the lower half-
$f$
-plane to the interior of a unit circle in the
$t$
-plane, and all singularities
$f_{A}+\mathbb{Z}$
to a single point
$t_{A}$
. By similar logic, the mirror singularity
$f_{\underline{A}}$
will also be accompanied by periodic copies on certain Riemann sheets, as in figure 10(c).
-
(iv) By solely encircling the combination
$A\underline{A}$
, diagonal singularities will appear in pairs, from
$\{\underline{B},\underline{B}^{\ell },\underline{C},\underline{C}^{\ell },\underline{D},\underline{D}^{\ell },\ldots \}$
, with a new type for each orbit of
$A\underline{A}$
. From the chain of
$A\underline{A}$
orbits, encircling
$A$
once further will cause the upper diagonal singularities
$\{B,B^{\ell },C,C^{\ell },D,D^{\ell },\ldots \}$
to appear in an analogous fashion.
Thus for example, on sheet
$A\underline{A}$
, singularities
$\{\underline{B},\underline{B}^{\ell }\}$
appear. On sheet
$A\underline{A}A\underline{A}$
, singularities
$\{\underline{B},\underline{B}^{\ell },\underline{C},\underline{C}^{\ell }\}$
appear, and so forth. This statement was the focus of § 7, and we also see the behaviour in figure 10(d,g). Adding the additional circling of
$A$
produces the mirror singularities; thus
$A\underline{A}A$
has singularities
$\{B,B^{\ell }\}$
, as shown in figure 10(f).
Indeed there are a plethora of other rules that may derive with the help of the continuation schemes discussed in § 4. However, we will end here with one final observation – the most vague and open of the ones we have proposed.
-
(v) Circling a diagonal singularity causes the generation of further singularities on the new sheet. The locations of these singularities are unpredictable and not subject to symmetry arguments in the
$f$
-plane.
The observation was discussed in § 8.1, and can be seen in figure 10(e,g,h). In such figures we see that upon circling a diagonal singularity, new singularities may arise (marked as a dark node). Such singularities arise in a non-symmetric fashion.
9 Conclusion
Using numerical computations of the analytic continuation of the Stokes wave, we have confirmed the existence of the main crest singularity, as shown by Grant (Reference Grant1973), Schwartz (Reference Schwartz1974) and many others. We have also confirmed the existence of its mirror image, as first conjectured by Tanveer (Reference Tanveer1991). Further exploring deeper layers of the Riemann surface reveals a countably infinite number of diagonal singularities.
10 Discussion
Throughout his numerous correspondences with William Thompson (Lord Kelvin), Stokes made frequent note of his hesitation that the highest possible wave would indeed exhibit a crest angle of
$120^{\circ }$
. However by September 1880, he seemed to have convinced himself of the theory’s veracity:
You ask if I have done anything more about the greatest possible wave. I cannot say that I have, at least anything to mention mathematically. For it is not a very mathematical process taking off my shoes and stockings, tucking up my trousers as high as I co[u]ld, and wading out into the sea to get in line with the crests of some small waves that were breaking on a sandy beach (
$\ldots$
) it did seem to me that the waves began to break while their sides still made only a blunt angle, a good deal less than
$90^{\circ }$
. I feel pretty well satisfied that the limiting form is one presenting an edge of
$120^{\circ }$
. (Stokes Reference Stokes and Wilson1880c
, p. 498)
While Stokes may have been able to reassure himself through direct observations of physical waves, for us, such interactions are more difficult – not least because the singularities we study are confined to the complex plane. However, a connection with the physicality of the problem need not be out of reach: through our construction and exploration of the Riemann surfaces, we have sought to offer the same sense of discovery that Stokes had found that day, wading into the sea.
It is furthermore remarkable that, perhaps the simplest possible problem involving the study of a nonlinear water wave, a topic of study that dates back over a century and a half, continues to confound us in the present day with new challenges. While we have attempted to highlight the complexity of the underlying issues, there are naturally three questions that follow from this work.
(i) What is the complete mathematical description of the Riemann surface for the Stokes wave? An exhaustive search of the individual Riemann sheets of the surface is yet to be undertaken. Here, we have sought to describe the patterns and connections between the branch-cut transitions for the first few Riemann sheets adjacent to the physical fluid. While this has led to the particular example shown in figure 10, and our proposed rules and observations of § 8, a more complete description of the full surface topology will be the subject of forthcoming work.
(ii) What are the asymptotic properties of the singularities in the two limits of
${\it\epsilon}\rightarrow 0$
(small waves) and
${\it\epsilon}\rightarrow 1$
(steep waves)? Though the case
${\it\epsilon}\rightarrow 0$
, may seem relatively benign, we still do not understand the various properties of the singularities in this limit. For example, it seems from figure 9 that the diagonal singularities tend to infinity along a constant angle as
${\it\epsilon}\rightarrow 0$
. Can these angles, as well as other properties of the Riemann surface be analytically predicted? There are strong connections between these observations and the occurrence of periodic arrays of singularities in other studies on nonlinear differential equations (see e.g. Costin & Costin Reference Costin and Costin2001, Chapman et al.
Reference Chapman, Trinh and Witelski2013).
Clearly, in the limit
${\it\epsilon}\rightarrow 1$
, the cancellation of the leading square-root power cannot be described using something as trivial as
$$\begin{eqnarray}\displaystyle z\sim \left[\vphantom{\tilde{b}_{0}}a_{0}+a_{1}(f-f_{A})^{1/2}+\cdots \,\right]+\left[\tilde{b}_{0}-a_{1}(f+f_{A})^{1/2}+\cdots \,\right], & & \displaystyle\end{eqnarray}$$
between the main crest singularity and its mirror image, and with
$f_{A}\rightarrow 0$
. For one reason, we have established that the mirror singularity is of order
$3/2$
. For another, encircling both singularities and climbing the
$A\underline{A}$
structure does not return the continuation to its original value – the topology of the Riemann surface is infinitely sheeted along this path. In contrast, two rotations in (10.1) will return the solution to its original value. Thus, any asymptotic description of the matching region between outer (away from the singular crest) and inner solutions cannot be composed of any finite number of square-root sums, as in (10.1). Indeed this was the observation from a preprint by Lushnikov (Reference Lushnikov2015), who conjectured that an infinite product of increasingly nested square roots might explain the recombination of the two main crest singularities to Stokes’ single crest singularity. As noted by the author, such a functional form would necessarily have additional off-axis singularities. It is intriguing that we have found such singularities here as well.
The work by Longuet-Higgins & Fox (Reference Longuet-Higgins and Fox1977, Reference Longuet-Higgins and Fox1978) studying the ‘almost-highest wave’ is most relevant, and we hope to undertake a comprehensive study reviewing the connections between the previous work and our present paper. Problems with coalescing singularities often reveal challenging asymptotic structures for study (Chapman et al. Reference Chapman, Trinh and Witelski2013; Trinh & Chapman Reference Trinh and Chapman2015).
(iii) What are the physical and practical consequences of the newfound singularity structure? Finally, the question most relevant to the practitioner (numerical, theoretical and experimental) is what we may take away from this study, as it relates to real life consequences of Stokes wave theory. We had already discussed, in § 1.1, many of the issues where better understanding of the singularity structure would illuminate. As many theories on time-dependent instabilities rely upon the interaction of the nearest singularity to the free surface, the question of how to account for the coalescence of the (presumed) infinity of singularities as
${\it\epsilon}\rightarrow 1$
seems to be of paramount importance. In addition, we also highlight intriguing questions on the role of Stokes wave singularities in the development of capillary ripples riding on steep gravity waves (cf. Longuet-Higgins Reference Longuet-Higgins1995). As shown in Trinh & Chapman (Reference Trinh and Chapman2013a
,Reference Trinh and Chapman
b
), it is expected that the singular limit of small surface tension will present a challenging route for further study.
Acknowledgements
We are indebted to Professors S. J. Chapman and H. Woudhuysen (Oxford) for their insightful suggestions and for providing close support during the preparation of this work. P.H.T. thanks Miss C. Goodwin for providing editorial input. Both authors express gratitude to Lincoln College for financial support.
