On the Limiting Stokes' Wave of Extreme Height in Arbitrary Water Depth

As mentioned by Schwartz (1974) and Cokelet (1977), it was failed to gain convergent results of limiting Stokes' waves in extremely shallow water by means of perturbation methods even with the aid of extrapolation techniques such as Pad\'{e} approximant. Especially, it is extremely difficult for traditional analytic/numerical approaches to present the wave profile of limiting waves with a sharp crest of $120^\circ$ included angle first mentioned by Stokes in 1880s. Thus, traditionally, different wave models are used for waves in different water depths. In this paper, by means of the homotopy analysis method (HAM), an analytic approximation method for highly nonlinear equations, we successfully gain convergent results (and especially the wave profiles) of the limiting Stokes' waves with this kind of sharp crest in arbitrary water depth, even including solitary waves of extreme form in extremely shallow water, without using any extrapolation techniques. Therefore, in the frame of the HAM, the Stokes' wave can be used as a unified theory for all kinds of waves, including periodic waves in deep and intermediate depth, cnoidal waves in shallow water and solitary waves in extremely shallow water.


Introduction
The two-dimensional steady progressive gravity wave is one of the most classic problems in fluid mechanics, which can be tracked back to Stokes (1847Stokes ( , 1880 and was widely studied by lots of researchers (Michell 1893;Nekrasov 1920;Yamada 1957;Yamada & Shiotani 1968;Schwartz 1972;Byatt-Smith & Longuet-Higgins 1976;Vanden-Broeck & Schwartz 1979;Chen & Saffman 1980;Olfe & Rottman 1980;Schwartz & Fenton 1982;Sulem et al. 1983;Hunter & Vanden-Broeck 1983;Vanden-Broeck 1986;Klopman 1990;Fenton 1990;Karabut 1998;Dallaston & Mccue 2010;Lushnikov 2016;Lushnikov et al. 2017). Among analytic approaches for this problem, perturbation methods are used most frequently. Stokes (1847Stokes ( , 1880 proposed a perturbation approach using the first Fourier coefficient, a 1 , as the perturbation quantity, and then showed that the highest free-surface wave (i.e. limiting wave, or extreme wave) in deep water would have a sharply pointed crest, enclosing a 120 • angle. Schwartz (1974) carried out this expansion for deep-water wave to the order 70 and found that, as the wave height H increases, the first Fourier coefficient a 1 first increases until it reaches a peak value, and then decreases. In other words, a single a 1 corresponds to two different wave heights for large enough wave height H. Thus, Stokes' expansion for deep-water waves is invalid for the limiting/extreme wave height.
Then Schwartz (1974) used a new expansion parameter ǫ = H/2 in his perturbation approach, and carried out the perturbation expansion to the 117th order in deep water and to the 48th order in general water depths, respectively. Utilizing the Padé approximants and the Shanks's iterated e 1 transformations (Shanks 1954), Schwartz (1974) successfully obtained convergent results for the ratio of water depth to wavelength d/λ > 0.05. However, his method has to rely on the extrapolation to obtain the dispersion relation for very high wave, since his perturbation series for the square of phase velocity, c 2 , only converges well for wave height shorter than 97% of the maximum. In addition, Schwartz (1974) found that accurate wave profile cannot be obtained for very high waves even with the aid of Padé approximant, so he added some standard terms to the crest to account for the remainder of the profile. Note that Schwartz (1974) ascribed the failure of his perturbation method in shallow water to round-off error.
A new perturbation quantity where v crest and c 0 are the fluid speed at the crest in the reference frame moving with the wave and the speed of waves of infinitesimal amplitude, respectively, was considered by Longuet-Higgins & Fenton (1974). Using this perturbation quantity, Longuet-Higgins & Fenton (1974) found that the series under the use of Padéapproximants converges better than that using ǫ = H/2. Further, another expansion parameter where v trough and c denote the fluid speed at the trough and the phase speed in the inertial frame, respectively, was used by Longuet-Higgins (1975). The computational efficiency was drastically improved by using this perturbation quantity. In addition, Longuet-Higgins (1975) proposed an alternative expansion parameter (1.3) Using (1.3) as the perturbation quantity, Cokelet (1977) carried out the expansion to the 120th order, and obtained convergent results for Stokes waves with d/λ > 0.0168. However, Cokelet (1977) pointed out that his method cannot give accurate wave profiles even in case of d/λ < 0.11. Furthermore, Dallaston & Mccue (2010) reconsidered both Schwartz's (Schwartz 1974) and Cokelet's (Cokelet 1977) schemes, but with exact calculations so as to void any round-off error. However, they found that both the series expansions of Schwartz (1974) and Cokelet (1977) actually cannot provide precise estimates of the limiting wave properties in extremely shallow water. Besides perturbation methods (Schwartz 1974;Cokelet 1977), a variety of numerical methods were proposed for the limiting Stokes' wave. One common numerical method is to minimize the mean-squared error in the kinematic and dynamic free surface boundary conditions. Chappelear (1961) expanded the velocity and the profile equation in Fourier series, and then used the method of least squares to determine the Fourier coefficients. Dean (1965) employed an analytical stream function expression with a series of unknown coefficients to describe the waves, and then used a numerical perturbation method to determine the unknown coefficients. Similarly, the numerical method was used by Williams (1981) to minimize the error in the surface boundary conditions over a series of evenly spaced points. However, a new crest term was supplemented to the integral equation in Williams' numerical method. Williams (1981) found that introducing this new term can greatly accelerate the convergence, i.e., the same level of accuracy can be reached by less terms of Fourier coefficients. This method (Williams 1981) was a significant progress in numerics for Stokes' wave, since it is one of few methods that are both free from extrapolation and can give accurate results. Unfortunately, this method (Williams 1981) still fails for the extremely shallow water d/λ < 0.0168. Rienecker & Fenton (1981) used a finite Fourier series to reduce the free surface conditions to a set of nonlinear algebraic equations, and then used Newton's iteration method to solve these nonlinear equations. By means of this method, the equations are satisfied identically at a number of points on the surface. This method was further simplified by Fenton (1988) who numerically solved all the necessary partial derivatives. However, Fenton (1988) found that it is sometimes still necessary first to solve a sequence of lower waves and then to extrapolate forward in height steps to reach the desired height. Vanden-Broeck & Schwartz (1979) proposed an efficient numerical scheme to solve the steep gravity wave. They first formulated the steep gravity waves as a system of integrodifferential equations, and then used the Newton's iteration technique to solve the coupled equations. Using this numerical method, accurate results can be obtained even in the case of d/λ = 0.008. In addition, Vanden-Broeck & Miloh (1995) employed series truncation methods, which use a refinement of Davies-Tulin's approximation (Davies 1951;Tulin 1983), to solve the steep gravity waves. By means of these methods, accurate numerical results can be obtained in the cases of d/λ 0.0168. It should be emphasized that these schemes are easier to implement than the boundary integral equation methods (Hunter & Vanden-Broeck 1983).
Besides the property of limiting Stokes wave with a included 120 • angle in the crest, the non-monotonicity of the speed and energy near the limiting wave height, first found by Longuet-Higgins (1975), also received wide attention from researchers. Longuet-Higgins & Fox (1978) proposed a matching technique for gravity waves of almost extreme form, and then successfully confirmed the existence of branch-points of order 1/2, as predicted by Grant (1973), and of turning-points in the phase velocity as a function of wave height. In addition, the asymptotic solution of Longuet-Higgins & Fox (1978) indicates that there is an infinite number of turning-points in the dispersion relation, momentum and energy for the wave height very close to maximum height. However, many methods generally only capture one or two of these turning points (Chandler & Graham 1993), but three turning points are found by Dallaston & Mccue (2010).
Note that understanding the characterization of the singularity structure of the Stokes' wave, such as the locations and scalings of the singularities, is of great help in theory (Tanveer 1991;Crew & Trinh 2016). Dyachenko et al. (2014Dyachenko et al. ( , 2016 analyzed the distance d c from the lowest singularity in the upper half-plane (i.e., the square-root branch point) to the real line which corresponds to the fluid free surface, and then suggested a power law scaling d c ∝ (H max − H) 3/2 . Using this power law scaling, Dyachenko et al. (2014Dyachenko et al. ( , 2016 presented an estimate H max /λ ≈ 0.1410633 in deep water. Moreover, a square-root branch point is found by Lushnikov (2016) to be the only singularity in the physical (first) sheet of Riemann surface for non-limiting Stokes wave. Then infinite number of square root singularities are found in the infinite number of non-physical sheets of Riemann surface after crossing a branch cut of a square root into the second and subsequently higher sheets of the Riemann surface. Furthermore, Lushnikov (2016) conjectured that non-limiting Stokes wave at the leading order consists of the infinite product of nested square root singularities, and that as increasing the steepness of the Stokes wave to the extreme form, these nested square root singularities will simultaneously approach to the real line from different sheets of Riemann surface and finally form together 2/3 power law singularity of the limiting Stokes wave. This conjecture was well supported by high precision simulations. In addition, the slow decay of the Fourier coefficients is a challenging problem for numerical methods due to the existence of the singularities for limiting/approximately-limiting Stokes' wave. In order to move all complex singularities away from the the free surface, Lushnikov et al. (2017) introduced a free parameter into an auxiliary conformal mapping so as to allow finer resolution near the crest of the wave. They found that the numerical convergence rate is dramatically improved by adapting the numerical grid near singularities.
Up to now, there are only few methods (Schwartz 1974;Cokelet 1977;Williams 1981;Vanden-Broeck & Schwartz 1979) capable to solve the two-dimensional limiting (extreme) steady progressive wave in very shallow water. Besides, almost all analytic/numerical methods fail to give accurate results (especially the wave profile) for limiting waves in extremely shallow water, such as d/λ < 0.005. In addition, most analytic/numerical methods rely on the extrapolation techniques, such as the Padé approximant, so as to accelerate the convergence and remove singularities that limit a series radius of convergence. So, an approach that can yield accurate results for the two-dimensional limiting (extreme) progressive gravity wave in arbitrary water depth without using any kind of extrapolation techniques is of great value. This is the motivation of this paper.
In this paper, the limiting Stokes' wave in arbitrary water depth is successfully solved by an analytic approximation method, namely the homotopy analysis method (HAM) (Liao 1992(Liao , 1999(Liao , 2003(Liao , 2009(Liao , 2010(Liao , 2012Van Gorder & Vajravelu 2008;Mastroberardino 2011;Kimiaeifar et al. 2011;Vajravelu & Van Gorder 2012;Sardanyés et al. 2015;Liao et al. 2016;Zhong & Liao 2017, 2018Liu et al. 2018a,b). Unlike perturbation methods (Schwartz 1974;Cokelet 1977), the HAM is independent of any small/large physical parameter. Especially, different from all analytic approximations, there is a so-called "convergence-control parameter" in the frame of the HAM, which has no physical meaning but provides us a convenient way to guarantee the convergence of series solutions. For example, perturbation techniques are invalid for large deformation of Von Kármán plate, a famous classic problem in solid mechanics. However, Zhong & Liao (2017, 2018 successfully applied the HAM to gain convergent series solution even for extremely large deformation of Von Kármán plate. It is worthwhile mentioning that some mathematical theorems of convergence have been rigorously proved in the frame of the HAM (Liao 2012). For instance, it has been proved by Liao (2012) that the power series given by the HAM converges to 1/(1 + t) in the intervals: respectively. So, the power series (1.4) converges to 1/(1 + t) either in the interval (−1, +∞) if letting < 0 impend 0, or in the interval (−∞, −1) if letting > 0 tend to 0, respectively. In other words, the power series (1.4) given by the HAM can converge to 1/(1 + t) in its entire definition interval (except the singularity t = −1) by properly choosing the so-called "convergence-control parameter" . This is a good example to illustrate that the HAM can dramatically improve the convergence of a series by means of the so-called convergence control parameter. By contrast, the traditional power series: only converges in the interval (−1, 1). Thus, the HAM can indeed greatly enlarge the convergence interval of solution series by means of properly choosing the so-called "convergence-control parameter" . Note that perturbation methods for many problems have been found to be a special case of the HAM when = −1, as illustrated by Zhong & Liao (2017, 2018, and this well explains why perturbation results are often invalid in case of high nonlinearity. In addition, the HAM provides us great freedom to choose initial approximation so that iteration can be easily used to accelerate convergence in the frame of the HAM. Besides, a better initial guess can also modify the convergence of iteration, too. More importantly, something completely new/different have been successfully obtained by means of the HAM: the steady-state resonant waves were first predicted by the HAM in theory (Liao 2011;Xu et al. 2012;Liu & Liao 2014) and then confirmed experimentally in a laboratory (Liu et al. 2015): all of these illustrate the novelty and potential of the HAM for highly nonlinear problems.
There exist two challenges for the traditional perturbation methods: (1) the series solutions diverge either when the water depth is rather small or when the wave height approaches to the peak value; (2) the computational efficiency is pretty low when the terms of Fourier coefficients are large. For the first challenge, it is found that the convergence of series solutions of the limiting Stokes wave can be guaranteed by means of choosing a proper convergence-control parameter in the frame of the HAM. Note that, since there is a singularity exactly locating at the crest for Stokes' wave of extreme form, considering many enough terms of Fourier series is inevitable if results of high precision are required. For the second challenge, we used an iteration HAM approach to greatly accelerate convergence of all unknown Fourier coefficients. Thus, by means of an iteration HAM approach with properly choosing convergence-control parameter , accurate results in arbitrary water depth can be obtained efficiently. More importantly, since all Fourier coefficients obtained by the HAM are convergent, accurate wave profiles in very shallow water can be presented without using any kinds of extrapolation techniques. Compared with the perturbation methods (Schwartz 1974;Cokelet 1977), our HAM approach is simpler, more easy-to-use and valid in almost the whole range of physical parameters, as mentioned later in this paper. All of these demonstrate the superiority of the HAM over the perturbation method for this famous problem in fluid mechanics.
This paper is organized as follows. Fundamental equations are given in § 2. Procedures of the HAM for the limiting Stokes wave problem are presented in § 3. The limiting (extreme) Stokes wave in infinite depth is considered in § 4. The limiting Stokes waves in finite depth are investigated in § 5, with comparison to previous results (Laitone 1960;Fenton 1972;Schwartz 1974;Cokelet 1977;Williams 1981). Concluding remarks are given in § 6.

Mathematical description of limiting Stokes' wave
Consider symmetrical, two-dimensional, periodic gravity waves propagating from right to left in a fluid with a horizontal bottom. The propagation of waves is only under the influence of gravity. Wave speed, c, is constant relative to an inertial frame. Assume that the fluid is inviscid and incompressible and that the motion is irrotational. Consider another reference frame moving with a wave crest. With respect to this frame, the motion is steady.
As shown in figure 1 (a), λ, H, g represent the wavelength, the wave height and the gravity acceleration, respectively. Locate the x axis at a distance d above the bottom. Let the stream function Ψ = 0 on the free surface, and Ψ = −c d on the horizontal bottom. The Bernoulli condition on the free surface reads vv + 2gy = K, where velocity v = v x − iv y , the bar denotes the complex conjugation, and K is an unknown constant, respectively. As shown in figure 1, we map the interior of fluid motion "ABODEA" in the physical "z" plane into an annulus "ABODEA" in the "ζ" plane according to the transformation: where ζ = Re iθ ; R, r 0 and θ represent the radius, inner radius and argument respectively; a 1 , a 2 , · · · , a j , · · · are unknown constant coefficients to be computed. The horizontal bottom Ψ = −c d and the free surface Ψ = 0 are then mapped onto the circles R = r 0 = e −d and R = 1, respectively. Note that r 0 = 0 and r 0 = 1 correspond to the cases of infinite depth and infinite shallow water, respectively. The complex velocity potential w can be expressed as where Φ represents velocity potential. According to (2.2), we have (2.4) So, we have the wavelength = 2π, (2.5) and the wave steepness Theoretically, a 1 , a 2 , · · · , a j , · · · need to be all reserved to identically satisfy the equation (2.10). However, we can only consider limited terms in practice. Thus, let us consider here the first r Fourier coefficients a 1 , a 2 , · · · , a r , i.e., f is approximated by Substituting (2.11) into (2.10) and then equating the coefficients of cos(kθ), where k = 0, 1, 2, · · · , r, we obtain the following (r + 1) algebraic equations †: 12) and N k [a 1 , a 2 , · · · , a r ] where N k (k = 1, 2, · · · , r) denotes a nonlinear operator, with the following definitions: a n1 a n−n1 when 1 n r, a n1 a n−n1 when r < n 2r. (2.14) Then, the next step is to solve the nonlinear algebraic equations (2.13) for the r unknown constant Fourier coefficients a 1 , a 2 , · · · , a r . Thereafter, the wave speed c can be directly given by (2.12).

The mathematical approach based on the HAM
Let a j,0 denote the initial guess of a j (j = 1, 2, · · · , r), a non-zero auxiliary parameter (called the convergence-control parameter), and q ∈ [0, 1] the embedding parameter for a homotopy, respectively. First of all, we construct a family of equations , Ω 2 (q), · · · , Ω r (q) , k = 1, 2, · · · , r, (3.1) where the nonlinear operators N 1 , N 2 , · · · , N r are defined by (2.13), and the unknown functions Ω 1 (q), Ω 2 (q), · · · , Ω r (q) correspond to the unknown constant Fourier coeffi-cients a 1 , a 2 , · · · , a r , respectively, and a k,0 is the initial guess of a k . Note that, in the frame of the HAM, we have great freedom to choose the initial guess a k,0 so that an iteration approach can be proposed based on this kind of freedom to greatly accelerate convergence, as mentioned later. More importantly, the so-called convergence-control parameter can provide us a simple way to guarantee the convergence of solution series, as shown below.
Assuming that is properly chosen so that the Maclaurin series (3.4) exists and converges at q = 1, then according to (3.3), we have the so-called homotopy-series solutions a n = +∞ k=0 a n,k , n = 1, 2, · · · , r. (3.7) Substituting (3.4) into the zeroth-order deformation equations (3.1) and then equating the like-power of q, we have the so-called mth-order deformation equations where Note that, in the frame of the HAM, we have great freedom to choose the initial guesses a 1,0 , a 2,0 , · · · , a r,0 . So, we can simply choose a k,0 = 1 k , k = 1, 2 · · · , r. (3.11) Then a 1,k , a 2,k , · · · , a r,k can be obtained by (3.9) step by step, starting from k = 1. The nth-order homotopy approximations of a 1 , a 2 , · · · , a r read Ω i,n = n k=0 a i,k , i = 1, 2 · · · , r. (3.12) Once a 1 , a 2 , · · · , a r are determined, the wave speed c can be given by (2.12).
In order to characterize the global error of our HAM approximation, we define the following squared residual error  Table 1. The squared residual error E , wave steepness H/λ and wave speed parameter (gλ)/ 2πc 2 versus iteration times in the case of r0 = 0, given by the first-order HAM-based iteration approach using c0 = −0.2, r = 100, and the initial guess (3.11).
smaller the E, the more accurate the HAM approximation (3.12). Besides, it has been proved (Liao 2003(Liao , 2012 in general that a homotopy-series converges to solution of original equations as long as all squared residual errors tend to zero. So, it is enough to check the squared residual error (3.13) only.

The limiting Stokes' wave in infinite depth
To show the validity of our HAM approach mentioned above, we first of all give convergent series solution of the limiting (extreme) Stokes' wave in infinite depth.
According to Liao (2003), the convergence of the homotopy-series solutions can be greatly accelerated by introducing the iteration technique, which uses the nth-order homotopy-approximationΩ 1,n ,Ω 2,n , · · · ,Ω r,n as new initial guesses a 1,0 , a 2,0 , · · · , a r,0 for the next iteration, say, a 1,0 =Ω 1,n , a 2,0 =Ω 2,n , · · · , a r,0 =Ω r,n . This provides us the nth-order iteration of the HAM. According to our computation, both the HAM approach without iteration and the HAM-based iteration approach can yield convergent results, but the efficiency of the HAM-based iteration approach is much higher. In particular, the first-order HAM-based iteration approach has the highest efficiency. So, we use the firstorder HAM-based iteration approach in all cases of this paper, if not specially mentioned. Table 1 presents the results in the case of r 0 = 0, corresponding to infinite depth of water, given by the convergence-control parameter = −0.2, r = 100 (i.e. one hundred truncated terms of Fourier series) and the initial guess (3.11). Note that the squared residual error E defined by (3.13) quickly decreases to the tiny level 10 −17 . This illustrates that all Fourier coefficients a 1 , a 2 , · · · , a 100 , given by our HAM approach, are convergent. Figure 2 shows the homotopy-approximation of the first Fourier coefficient, a 1 , versus iteration times in the case of r 0 = 0, given by r = 100 and the convergence-control parameter = −0.4, −0.25, −0.1, respectively. Note that = −0.4 leads to divergence of iteration, = −0.25 corresponds to a quickly convergent iteration, but = −0.1 a slowly convergent iteration, respectively. Obviously, the optimal value of corresponds to the fastest convergence, as pointed out by Liao (2010). It is found that convergent results can be obtained by our iteration HAM approach with arbitrary values of ∈ [−0.27, 0). So, the convergence-control parameter indeed provides us a simple way to guarantee convergence and to accelerate convergence. This clearly illustrates the important role of the convergence-control parameter in the frame of the HAM. Table 2 presents the convergent results in the case of r 0 = 0, given by different r. Note that the steepness of the limiting wave in infinite water depth tends to a fixed value H/λ = 0.14108 when r is large enough, say, r > 5000. This is reasonable, since the precision of our results is controlled by r, i.e. the truncated number of the Fourier  Table 2. Wave steepness H/λ and wave speed parameter (gλ)/ 2πc 2 versus truncated terms r in the case of r0 = 0 (in infinite depth), given by the first-order HAM-based iteration approach using c0 = −0.2. series (2.11). Note that Schwartz (1974) gave H max /λ = 0.14118 but Dyachenko et al. (2016) gave H max /λ = 0.141058 for limiting wave in deep water, with 0.071% and 0.016% relative errors compared to our results, respectively. It should be emphasized that, by means of the HAM, convergent results of all Fourier coefficients a j can be obtained. This distinguishes the HAM from other methods. It is found that the high-order Fourier coefficients a j drop rather slowly (e.g., a 1 = 0.29223, a 100 = 0.01576, a 500 = 0.005415, a 1000 = 0.003739, a 3000 = 0.002905, a 5000 = 0.002759). We have H/λ = 0.14085 even by means of r = 500, and have the more accurate result H/λ = 0.14108 by r > 5000. All of these indicate that f defined by (2.8) converges pretty slowly indeed. However, it should be emphasized that, whatever r we choose, convergent values of all Fourier coefficients a j can be directly obtained by our iteration HAM approach without using any extrapolation and Padé approximant techniques. Figure 3 shows the comparison of limiting wave profiles given by Schwartz (1974), Dyachenko et al. (2016) and our HAM approach. The agreement between them is satisfactory. This indicates the validity of our HAM-based approach. Our limiting wave profile has a sharply pointed crest with an enclosing angle 119.3 • , which is very close to the theoretical value 120 • . However, Schwartz (1974) mentioned that "while the method of Páde fractions yields accurate profiles for wave heights somewhat short of the maximum, it is insufficient for the description of very high waves", and that "Páde fractions do not converge well in the immediate neighbourhood of branch-points; moreover, only the first few coefficients a j , can be determined with acceptable accuracy". Thus, Schwartz (1974) had to use the so-called "series completion method" to gain a satisfactory wave profile. By contrast, using the HAM, convergent results of all Fourier coefficients a j (j = 1, 2, 3, · · · , r) and the convergent wave profile for the limiting wave can be obtained without using any extrapolation techniques such as Padé technique, the series completion method and so on. This illustrates that the HAM-based approach is superior to perturbation methods (Schwartz 1974;Cokelet 1977). This is mainly because, unlike perturbation methods, the HAM provides us a convenient way (through the socalled convergence-control parameter ) to guarantee the convergence of solution series.

The limiting Stokes' wave in finite depth
Note that r 0 = 0 and r 0 = 1 correspond to the case of infinite depth and infinite shallow water, respectively. Without loss of generality, let us first consider the case of r 0 = 0.05. In the frame of the HAM, we have great freedom to choose the initial guesses of a 1 , a 2 , · · · , a r . Considering the continuous variation of a 1 , a 2 , · · · , a r as r 0 increases from 0 to 1, the convergent values of a 1 , a 2 , · · · , a 5000 in the case of r 0 = 0, obviously, are much better than (3.11) as the initial guess for the case of r 0 = 0.05. In other words, if we have obtained the convergent results of a 1 , a 2 , · · · , a 5000 in the case of r 0 = 0, then  Table 3. Wave steepness H/λ and wave speed parameter (gλ)/ 2πc 2 versus iteration times, m, in the case of r0 = 0.05, given by the first-order HAM-based iteration approach using the convergence-control parameter c0 = −0.2, the truncated terms r = 5500 and the initial guess (3.11).  Table 4. Wave steepness H/λ and wave speed parameter (gλ)/ 2πc 2 versus iteration times, m, in the case of r0 = 0.05, given by the first-order HAM-based iteration approach using the convergence-control parameter c0 = −1.2, the truncated terms r = 5500 and the initial guesses (5.1).
it is better to take a k,0 = a k when 1 k 5000, a 5000 when k > 5000, (5.1) as the initial guesses of a 1 , a 2 , · · · , a r in the case of r 0 = 0.05. It is found that, in the case of r 0 = 0.05, the optimal convergence-control parameter is about −0.2 if the initial guess (3.11) for r 0 = 0 is taken, and 400 times iteration is required to gain convergent results H/λ = 0.14026 and (gλ)/(2πc 2 ) = 0.8421, as shown in Table 3. However, if we take the initial guess (5.1), the optimal convergent-control parameter becomes −1.2, and we obtain the same convergent results H/λ = 0.14026 and (gλ)/(2πc 2 ) = 0.8421 by just thirty times iteration, as shown in Table 4. Thus, the computational efficiency by means of the initial guess (5.1) is approximately 13 times higher than that by (3.11). This illustrates that our iteration HAM approach with the optimal convergence-control parameter can indeed greatly accelerate the convergence.
Similarly, the convergent results in arbitrary water depth are successfully obtained by means of the above-mentioned strategy, as shown in Table 5. Note that Liao (2010) suggested a general approach to gain an optimal convergence-control parameter in the frame of the HAM. According to our computation, the interval of , which guarantees the convergence of iteration, becomes larger with the increase of r 0 . It is found that, in the case of 0.05k < r 0 0.05(k + 1), where 0 k 19 is a natural number, the corresponding optimal convergence-control parameter can be expressed by the following empirical formula  Table 5. Results for a variety of water depths, given by the first-order HAM-based iteration approach.
if we use the known convergent Fourier coefficients a j in the case of r 0 = 0.05k as the initial guess. Note that a convergence-control parameter closer to 0 represents a slower convergence of solutions, i.e., a lower efficiency of computation, as shown in Figure 2. Thus, the convergence-control parameter provides us a convenient way not only to guarantee the convergence of series solutions but also to improve the computational efficiency. It is found that, for all cases considered in Table 5, a few hundred times of iteration are enough to gain convergent results of all Fourier coefficients a j . In case of extremely shallow water, a huge number of Fourier coefficients are needed to present the limiting wave with the sharp crest. Table 6 presents the convergent results given by different values of r in the case of r 0 = 0.99. It is found that r = 50000 can give the fixed results H/λ = 1.3281 × 10 −3 and (gλ)/(2πc 2 ) = 60.175 in the case of r 0 = 0.99. This indicates that our iteration HAM approach can indeed give convergent results of the limiting Stokes' waves even in extremely shallow water. Note that, in case of r = 50000, we must solve a set of 50000 coupled, highly nonlinear algebraic equations! Fortunately, this is possible nowadays by means of a supercomputer such as TH-2 at National Supercomputer Centre in Guangzhou, China. Finally, it should be emphasized that all of these convergent results are obtained directly, say, without using any extrapolation and Padé approximant techniques.
Stokes (1880) gave a famous conjecture that the limiting wave (with extreme height) should have a sharp crest with an included angle 120 • . About one hundred years later, this conjecture was independently proved in mathematics by Amick et al. (1982) and Plotnikov (2002) for Stokes' waves in arbitrary depth of water. However, to the best of our knowledge, the detailed wave profiles for limiting Stokes' wave in extremely shallow water have not been reported.  Table 7. Included crest angles in a variety of depths, given by the first-order HAM-based iteration approach.
Stokes' wave in a variety of water depths, given by our iteration HAM approach. All of the included crest angles in different depth given by the HAM are very close to the theoretic value 120 • . The wave profiles in a variety of water depths given by the HAM are shown in Figure 4. Note that the high-order Fourier coefficients a j play an important role in correctly describing the wave profile, especially the wave crest. For instance, although Cokelet's perturbation method (Cokelet 1977) can give H/d with acceptable accuracy for r 0 < 0.9, however, it fails to give accurate wave profile even for r 0 > 0.5. Fortunately, the HAM can always yield convergent results of all Fourier coefficients by means of choosing a proper convergence-control parameter . This once again illustrates the superiority of the HAM over other methods (Schwartz 1974;Cokelet 1977). According to the convergent results given by the iteration HAM approach, we have the fitted formulas of H/d versus c 2 /(gd) and λ/d: which agree quite well with our HAM results, as shown in Figure 5. Let us make some comparisons of the limiting Stokes' waves given by different analytic/numerical methods. Table 8 presents the comparison of limiting wave steepness in a variety of depths. The results of Schwartz (1974) are accurate only for r 0 < 0.7; the results of Cokelet (1977) are accurate only for r 0 0.8; the results of Williams (1981) are of high accuracy for r 0 0.9, although they are slightly smaller. However, all of these methods fail to give convergent result for r 0 > 0.9, i.e., extremely shallow water. Fortunately, the HAM can give convergent results for the limiting waves almost in arbitrary depth.  9.75 × 10 −2 9.739 × 10 −2 9.7374 × 10 −2 9.7388 × 10 −2 0.5 7.91 × 10 −2 7.910 × 10 −2 7.9072 × 10 −2 7.9084 × 10 −2 0.6 6.14 × 10 −2 6.090 × 10 −2 6.0984 × 10 −2 6.0995 × 10 −2 0.7 4.5 × 10 −2 4.374 × 10 −2 4.3975 × 10 −2 4.3983 × 10 −2 0.8 --2.79 × 10 −2 2.8258 × 10 −2 2.8263 × 10 −2 0.9 --1.5 × 10 −2 1.3667 × 10 −2 1.3670 × 10 −2 0.95 ------6.7292 × 10 −3 0.97 ------4.0128 × 10 −3 0.99 ------1.3281 × 10 −3 Table 8. Limiting wave steepness, H/λ, in a variety of depths. Figure 6 shows the comparison of the limiting wave steepness H/λ, given by Schwartz (1974), Williams (1981) and the HAM approach mentioned in this paper, respectively. It is found that the perturbation method (Schwartz 1974) is only valid for r 0 ∈ [0, 0.7] even with the aid of extrapolation and Padé approximant techniques; Williams' numerical method (Williams 1981) is only valid for r 0 ∈ [0, 0.9]. However, the HAM can give accurate convergent results even for r 0 ∈ [0, 0.99]. Figure 7 shows the comparison of H/d versus the squared Froude number, c 2 /(gd). It is found that Cokelet's perturbation method (Cokelet 1977)  accord with the results of the highest solitary wave: (5.5) This suggests that the solitary wave theory could be unified into the Stokes' wave theory. According to Hedges (1995), waves with the Ursell number Hλ 2 /d 3 > 4000 are  (Schwartz 1974); , perturbation method with the aid of both Padé approximants and Shanks's iterated e1 transformation (Schwartz 1974); : Williams' numerical method (Williams 1981). regarded as solitary waves. It is found that, in the case of r 0 = 0.99, corresponding to λ/d ≈ 600, the Hλ 2 /d 3 of the limiting Stokes' wave given by the HAM reaches 3 × 10 5 . Thus, the Stokes' wave theory is actually valid almost in arbitrary depth, as shown in Figure 8. So, in the frame of the HAM, the Stokes wave theory can describe not only the periodic waves in deep and intermediate water but also cnoidal wave in shallow water and solitary wave in extremely shallow water. In addition, the ratio of wave height to depth, H/d, of the highest solitary wave was widely studied by many researchers: H/d = 0.827 was given by Yamada (1957), Lenau (1966), Yamada & Shiotani (1968), Longuet-Higgins & Fenton (1974; but H/d = 0.8332 was given by Witting (1981), Witting & Bergin (1981), Williams (1981), Hunter & Vanden-Broeck (1983). Note that, H/d = 0.8303 > 0.827 is given by the HAM in the case of r 0 = 0.99. Hence the value H/d = 0.827 for the highest solitary wave is denied by the HAM. Figure 9 shows the wave profiles of the limiting wave in the case of r 0 = 0.99 given by the HAM from the exact wave equations, the KdV solution (Korteweg & de Vries 1895), Laitone's second-order solution (Laitone 1960) and Fenton's ninth-order solution (Fenton 1972). It is found that only the HAM gives a wave profile with a sharply pointed crest, enclosing an angle 119.2 • . So, the KdV solution (Korteweg & de Vries 1895), Laitone's solution (Laitone 1960) and Fenton's solution (Fenton 1972) are all no longer valid in the limiting case. However, compared to the famous solitary solution of KdV equation (Korteweg & de Vries 1895) and Laitone's solution (Laitone 1960), Fenton's ninth-order solution (Fenton 1972) is of higher accuracy.
In summary, using the iteration HAM approach with a proper convergence-control parameter, we gain limiting Stokes' waves almost in arbitrary water depth, without using any extrapolation techniques. Therefore, in the frame of the HAM, the Stokes' wave theory is a unified theory for all kinds of progressive waves in arbitrary depth, even including solitary waves in extremely shallow water.

Concluding remarks
Obviously, limiting Stokes' wave in shallow water is a strong nonlinear problem. Previous methods, especially the perturbation methods, usually suffer divergence either when the wave height approaches the peak value or when the water depth is extremely small. For the limiting Stokes' wave, due to the existence of singularity locating exactly at the crest, perturbation methods usually can gain convergent results only for a small part of Fourier coefficients so that the extrapolation methods such as Padé approximant techniques and Shanks' transformation had to be used.
In this paper, we employ the homotopy analysis method (HAM) to solve the limiting Stokes' wave in arbitrary depth of water. It is found that the convergence of all Fourier coefficients of the solutions can be guaranteed by choosing a proper convergence-control parameter in the frame of the HAM, as shown in Figure 2. In addition, since the Fourier series is used to represent the free surface with a sharp pointed crest, using a large number of Fourier coefficients is inevitable. For other analytic/numerical methods, this might lead to rather slow convergence of the Fourier coefficients of the solutions. Fortunately, the HAM also provides us great freedom to choose initial guesses of solutions. Based on this kind of freedom of the HAM, we proposed an iteration HAM approach to greatly accelerate the convergence of all Fourier coefficients. Note that, since we consider a large enough number of Fourier coefficients, and more importantly, all of these coefficients are convergent without using any extrapolation methods, hence we can obtain the accurate wave profile even in rather shallow water.
It should be emphasized that accurate representation of the wave profile in very shallow water is impossible for other methods, especially without using any kind of extrapolation techniques. For instance, although Cokelet's perturbation method (Cokelet 1977) can give results of H/λ with acceptable accuracy for r 0 < 0.9, however, it can only give a good wave profile for r 0 0.5. Fortunately, by means of the HAM, we gain accurate limiting wave profiles in almost arbitrary depth of water, i.e., from r 0 = 0 to r 0 = 0.99, without using any extrapolation methods such as Padé approximant techniques and Shanks' transformation. To the best of our knowledge, accurate wave profile in the case of r 0 = 0.99 has been never reported. This once again illustrates the superiority of the HAM over perturbation and traditional numerical methods for this famous problem.
According to Hedges (1995), waves with the Ursell number Hλ 2 /d 3 > 4000 are regarded as solitary waves. It is found that, in the case of r 0 = 0.99, corresponding to λ/d ≈ 600, the Hλ 2 /d 3 of the Stokes' wave given by our HAM approach reaches 3 × 10 5 . Thus, the Stokes' wave theory is actually valid almost in arbitrary depth, as shown in Figure 8. So in the frame of the HAM, the Stokes wave theory can describe not only the periodic waves in deep and intermediate water but also cnoidal wave in shallow water and solitary wave in extremely shallow water. Therefore, in the frame of the HAM, the Stokes' wave is a unified theory for all kind of progressive waves, even including the limiting (extreme) solitary waves with a sharp crest of 120 • included angle in extremely shallow water! Note that the cubic relations between a j in equations (2.12)-(2.13) were considered in this paper, although the quadratic relations between the Fourier coefficients a j were reported by Longuet-Higgins (1978). Certainly, the computational efficiency could be improved by means of using the quadratic relations (Longuet-Higgins 1985;Balk 1996), but one should obtain the same results as mentioned above in this paper, from a physical viewpoint.
From viewpoint of applied mathematics, this paper provides us an additional example to illustrate that the HAM can be indeed applied to find something completely new, such as the discovery of the steady-state exactly/nearly resonant gravity waves with time-independent wave spectrum (Liao 2011;Xu et al. 2012;Liu & Liao 2014;Liu et al. 2015;Liao et al. 2016;Liu et al. 2018a), or to attack some challenging problems with high nonlinearity. j n h n + N 1 cos θ + N 2 cos(2θ) + · · · , (A 7) in which N 1 , N 2 , · · · , N k are defined by (2.13).