3.1 Solitary-wave dynamics

As shown in the aforementioned references, two basic families of symmetric gravity–capillary solitary waves, depression and elevation waves, are known to exist in deep water, and we briefly review some relevant results here. The bifurcation diagrams for depression and elevation branches together with typical wave profiles are presented in figures 1 and 2, from which we conclude that both branches bifurcate from infinitesimal periodic waves at the minimum of the phase speed $c^{\ast }=\sqrt{2}$ , and exist at subcritical speeds. Figure 1 shows that both the amplitude (where ‘amplitude’ is defined as the centre elevation) and the energy (which is given by (2.26)) of depression waves are monotonic functions of the translating speed $c$ . As $c$ decreases the wave steepness increases, while the wave crest becomes more rounded and eventually the wave profile develops an overhanging structure at the trough (e.g. the profile $1b$ ). Further down the curve, the profile develops a point of contact with a ‘trapped bubble’ (which occurs at $c\approx 0.9275$ ) and even a self-intersecting structure (e.g. the profile $1c$ ). The solutions with self-intersections have no physical meaning but are mathematically admitted by the full Euler equations. In contrast to the branch of depression, the bifurcation of the elevation branch exhibits the characteristic ‘snake-like behaviour’: the speed–energy curve (figure 2 b) has stationary points where $\text{d}E/\text{d}c=0$ (e.g. $2C$ in figure 2 b), and turning points where $\text{d}E/\text{d}c=\infty$ or equivalently $\text{d}c/\text{d}E=0$ (e.g. $2B$ in figure 2 b). Together with some wave profiles, the speed–amplitude curve of the elevation branch is shown in figure 2(a), which also features the snake-like behaviour. We stopped the computations after passing through the fourth turning point when two almost completely separated depression waves appear (see $2j$ in figure 2 a). At this point, more accurate computations became prohibitive, presumably since the waves are far apart and their separation is sensitive to perturbations. We finally remark that there are infinite many branches of symmetric gravity–capillary solitary waves existing at finite amplitudes (the read is referred to the numerical computations by Dias et al. (Reference Dias, Menasce and Vanden-Broeck1996), Wang et al. (Reference Wang, Vanden-Broeck and Milewski2014), and the asymptotic analysis of Yang & Akylas (Reference Yang and Akylas1997)).

The stability characteristics of basic branches were examined by Calvo & Akylas (Reference Calvo and Akylas2002) via a linear spectral analysis and by Milewski et al. (Reference Milewski, Vanden-Broeck and Wang2010) via direct numerical simulations. They found that all single-valued depression solitary waves are stable, while the elevation waves are unstable between the bifurcation point and the first turning point ( $2B$ in figure 2 b) but regain spectral stability beyond this point where the wave profiles feature two well-separated troughs (see $2c$ in figure 2 a). The existence of stable solitary waves brings up the question of the excitation of these steady solutions. In laboratories, depression gravity–capillary solitary waves were generated under controlled conditions by Longuet-Higgins & Zhang (Reference Longuet-Higgins and Zhang1997), Diorio et al. (Reference Diorio, Cho, Duncan and Akylas2009, Reference Diorio, Cho, Duncan and Akylas2011), among others. The experimentalists excited the depression waves by blowing air towards the surface of the fluid and moving the air source with a speed close to the phase speed minimum. We present here an inviscid numerical analogue of these experiments. The equation (2.20) was simulated with a fully localised pressure distribution

The time evolution of the transient response from a quiescent state, together with wave profiles, is illustrated in figure 3. Moving the frame of reference with the speed of the forcing, we observe a periodic shedding of depression solitary waves generated downstream of the pressure distribution. A similar phenomenon occurs for pure gravity waves in shallow water, where the KdV-type solitons can be periodically generated by a moving surface pressure, and the reader is referred to the comprehensive review paper by Keller (Reference Keller1985). The difference between these two cases is that the depression gravity–capillary solitary waves trail behind the pressure distribution, while the solitons for pure gravity waves advance upstream of the disturbance in processions. Despite being different in behaviour, two phenomena share an important core feature. In the transcritical range (namely, the speed slightly below the phase speed minimum for gravity–capillary waves in deep water, or the speed slightly above the phase speed maximum for pure gravity waves in shallow water) the dispersive effect is weak. When a fully localised pressure distribution moves with the speed in this range, there is neither steady forced solution nor linear mechanism to radiate energy away. Therefore the energy accumulates within the forced region, as a consequence, the amplitude of the local wave grows. When the amplitude reaches a certain level, the nonlinear effects start to become efficient. The subtle balance between weak dispersion and weak nonlinearity results in the formation of solitary waves moving at a different velocity from that of the forcing.

Figure 3.

(a) Shedding of depression solitary waves resulting from moving a fully localised forcing with the constant speed $c=1.41$ . (a) A forced solution in the $\unicode[STIX]{x1D709}-t$ plane clearly indicates that solitary waves are generated periodically to the left of the forcing. (b) Snapshots in the physical domain at times $t=500,1000,1500,2000$ from top to bottom respectively. The frame is moving with the forcing and only a portion of the computational domain is shown.

We now consider the generation of stable elevation solitary waves featuring two big troughs connected by a small dimple. Motivated by the structure of these waves, the pressure applied to the rest free surface consists of two identical, properly separated Gaussian functions

Snapshots in figure 4(a) show that this special pressure impinging on the fluid surface can produce a disturbance in the form of an elevation wave. The pressure was then released after a period of time ( $t=80$ ). The bottom snapshot in 4(a) and the top snapshot in 4(b) are shortly before ( $t=75$ ) and after ( $t=100$ ) the pressure was switched off respectively. The bottom snapshot in 4(b) indicates that the elevation solitary wave maintains its structure for long time in spite of the interaction with the background radiation field induced by the initial generation process.

Figure 4.

Excitation of a stable elevation solitary wave due to two fully localised, well-separated pressures moving at $c=1.405$ . The pressures are switched on at $t=0$ and switched off at $t=80$ . The solution (solid line) together with the pressures (dotted line) are shown in (a) at $t=15,45,75$ from top to bottom respectively. After the pressures are turned off, the solution at times $t=100,500,1000$ is shown in (b).

The numerical computations for generating figures 3 and 4 were performed on a periodic domain $\unicode[STIX]{x1D709}\in [-250,250)$ , with the time step $\unicode[STIX]{x0394}t=2\times 10^{-4}$ , $N=20\,000$ Fourier modes, dealiasing with a buffer of $N$ modes and no filtering. The character of the stable solution is not sensitive to the radiated waves wrapping around the domain.

The second dynamical computational experiment presented is the stability test of the depression solitary waves with an overhanging structure. As noted earlier, the numerical spectral analysis by Calvo & Akylas (Reference Calvo and Akylas2002) shows that depression waves are longitudinally stable. Unfortunately, their formulation is applicable only to the single-valued surface profiles, therefore the stability of the overhanging waves remains unknown. However this problem can be easily solved if one notices that the speed–energy curve is monotonic for the branch of depression and all perturbations are superharmonic in a sense for solitary waves. That’s because Saffman’s theorem indicates that for superharmonic perturbations stability exchanges can only occur at stationary points or turning points of the total energy. To provide more evidence for this assertion, we collided two depression waves: the larger slower wave was characterised by an overhanging structure, while the smaller faster wave served as a perturbation. Snapshots of the head-on collision in figure 5 show that both depression waves survive the collision, although they lose part of the energy, which eventually turns into a small wavepacket intercepted on the right-hand side of the larger solitary wave. The striking phenomenon is the persistence of the overhanging structure, although its shape slightly opens due to the energy loss during the collision process. In the case of overtaking collisions, the behaviours are different from those of head-on collisions, since the interacting time is much longer so that strong nonlinear effects are more likely to take place. An example is shown in figure 6, where only the larger solitary wave remains after the collision. As the smaller faster depression wave approaches the overhanging wave, it breaks up by transferring some energy to the overhanging structure and radiating the remaining part. We should emphasise that although the overhanging structure survives the overtaking collision, the gain of energy during the interaction results in a self-intersecting structure, therefore the solution satisfies the full Euler equations but has no physical significance. For the experiments of collisions, the equation (2.20) was integrated on a periodic domain $\unicode[STIX]{x1D709}\in [-100,100)$ with the time step $\unicode[STIX]{x0394}t=10^{-5}$ and $N=4000$ Fourier modes.

Figure 5.

Head-on collision between a solitary wave with an overhanging structure ( $1b$ in figure 1 a) and a small-amplitude depression wave ( $\unicode[STIX]{x1D702}(0)=-0.208$ and $c=1.405$ ): (a)  $t=0$ , (b) $t=40$ , (c) $t=50$ and (d) $t=70$ . The collision is inelastic, but both depression solitary waves survive the collision.

Figure 6.

Overtaking collision between an overhanging solitary wave ( $\unicode[STIX]{x1D702}(0)=-1.4676$ and $c=0.94$ ) and a small-amplitude depression wave ( $\unicode[STIX]{x1D702}(0)=-0.208$ and $c=1.405$ ) is shown in a frame of reference moving at the speed of the larger wave. From (a) to (d): $t=0,200,250,300$ . Only the larger wave survives the collision. Note that an unphysical wave profile shows up after the collision.

3.2 Linear spectral analysis

In this section, using numerical spectral analysis and direct numerical simulations, we revisit the branch of symmetric elevation solitary waves and investigate the stability of a typical branch of asymmetric solitary waves. Following Tanaka (Reference Tanaka1985) or Calvo & Akylas (Reference Calvo and Akylas2002), we first outline the formulation of the linear eigenvalue problem for gravity–capillary waves. We introduce a frame of reference moving to the right with speed $c$ , namely, we define $\widehat{x}=x-ct$ and $\widehat{\unicode[STIX]{x1D719}}=\unicode[STIX]{x1D719}-c\widehat{x}$ . Dropping hats, the dynamic boundary condition (2.6) can be rewritten as

while the kinematic boundary condition remains the same. As a consequence,

The steady solutions to the system (2.1), (2.2), (3.3) and (3.4) can be obtained by using hodograph transformation which exchanges the dependent variables $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D713}$ and independent variables $x$ and $y$ , and the Cauchy integral formula which solves Laplace’s equation in an explicit manner. Following Vanden-Broeck & Dias (Reference Vanden-Broeck and Dias1992), the horizontal and the vertical velocities can be expressed in terms of the complex potential $f=\unicode[STIX]{x1D719}+i\unicode[STIX]{x1D713}$ and $z=x+iy$ by the following formula

whereas $x+iy$ is an analytic function of $f$ in the upper half-plane ( $\unicode[STIX]{x1D713}>0$ ). Applying the Cauchy integral formula to the function $x_{\unicode[STIX]{x1D719}}+iy_{\unicode[STIX]{x1D719}}+(1/c)$ yields the relation between $x_{\unicode[STIX]{x1D719}}$ and $y_{\unicode[STIX]{x1D719}}$ at $\unicode[STIX]{x1D713}=0$

We rescale the Cauchy integral formula by defining $\unicode[STIX]{x1D719}\triangleq c\,\unicode[STIX]{x1D709}$ , then the expression (3.6) becomes

It is not difficult to verify that the kinematic boundary condition is satisfied automatically by substituting (3.5) into (2.2), while the dynamic boundary condition (3.3) can be rewritten as

Evidently, (3.7) and (3.8) are coincident with (2.10) and (2.25), and this fact allows us to use the obtained travelling-wave solutions. When a fully localised steady solution is given, we consider small perturbations of this solution. The solitary waves are supposed to be single valued, as well as the perturbations. Therefore we can express the free-surface displacement $\unicode[STIX]{x1D702}$ and the velocity potential $\unicode[STIX]{x1D719}$ as

where $Y(x)$ and $\unicode[STIX]{x1D6F7}(x,y)$ denote the profile of the free surface and the velocity potential of an exact solitary-wave solution, while $\widetilde{\unicode[STIX]{x1D702}}\text{e}^{\unicode[STIX]{x1D706}t}$ and $\widetilde{\unicode[STIX]{x1D719}}\text{e}^{\unicode[STIX]{x1D706}t}$ are small time-dependent perturbations superposed on these quantities respectively. Linearising the free boundary conditions (2.2) and (3.3), we obtain the eigenvalue problem for $\unicode[STIX]{x1D706}$ , $\widetilde{\unicode[STIX]{x1D702}}$ and $\widetilde{\unicode[STIX]{x1D719}}$ as follows:

Following Tanaka (Reference Tanaka1985) or Calvo & Akylas (Reference Calvo and Akylas2002), we introduce the magnitude of the velocity on the free surface, $q=\sqrt{\unicode[STIX]{x1D6F7}_{x}^{2}+\unicode[STIX]{x1D6F7}_{y}^{2}}$ , and the angle between the velocity and the $x$ -axis, $\unicode[STIX]{x1D703}=\arctan \;(\text{d}Y/\text{d}x)$ . The free surface can be parameterised by the arc length $s$ with $s=0$ being consistent with $\unicode[STIX]{x1D709}=0$ in (3.8). After some algebra, we obtain

where $\widetilde{\unicode[STIX]{x1D713}}$ is the harmonic conjugate of $\widetilde{\unicode[STIX]{x1D719}}$ . Noticing $\text{d}/\text{d}s=-q\text{d}/\text{d}\unicode[STIX]{x1D6F7}$ and $\text{d}/\text{d}\unicode[STIX]{x1D709}=c\text{d}/\text{d}\unicode[STIX]{x1D6F7}$ , one finally rewrites the linear eigenvalue problem as

where $X_{\unicode[STIX]{x1D709}}=-1+{\mathcal{H}}[Y_{\unicode[STIX]{x1D709}}]$ and $J=X_{\unicode[STIX]{x1D709}}^{2}+Y_{\unicode[STIX]{x1D709}}^{2}$ . We remark that under a different scaling, a similar formulation for left-going waves was previously used to study the stability of solitary waves for pure gravity waves in shallow water by Tanaka (Reference Tanaka1985) and for gravity–capillary waves in deep water by Calvo & Akylas (Reference Calvo and Akylas2002). In terms of the left- and right-going waves, the formulations given in Tanaka (Reference Tanaka1985) and presented here are essentially the same.

For consistency, the eigenvalue problem (3.14)–(3.15) is discretised at the uniform grid points of a periodic domain with the derivative and the Hilbert transform being approximated by the spectral matrices. Specifically, the square matrix of the discrete Hilbert transform, denoted by $\mathbb{H}$ , is given by

and the first-derivative matrix $\mathbb{D}$ is given by

where $L$ is the length of the computational domain and $N$ is the number of grid points (see Trefethen (Reference Trefethen2000) for detailed derivations). Therefore, the original eigenvalue problem is approximated by a matrix eigenvalue problem, which can be solved using the standard QR algorithm in MATLAB.

To validate our numerical code, we checked the results produced by our program against the results in Calvo & Akylas (Reference Calvo and Akylas2002). For the elevation gravity–capillary solitary waves, Calvo & Akylas (Reference Calvo and Akylas2002) considered the stability exchange near the first turning point, and found that the waves are unstable before that point. For $\unicode[STIX]{x1D6FC}$ (which is defined as $1/c^{4}$ ) equals $0.26$ , $0.3$ and $0.4$ , the growth rate of the most unstable mode is respectively $0.009$ , $0.0153$ and $0.0077$ . While our numerical results read $0.009$ , $0.0154$ and $0.01096$ after being converted to their scaling. The small discrepancies are most likely due to the different way to approximate infinitely long solitary waves, and the different matrices of the Hilbert transform and derivatives arising therefrom. The profile of the most unstable mode for $\unicode[STIX]{x1D6FC}=0.3$ was also checked (see figure 7), which is in excellent agreement with Calvo & Akylas’ result (figure 10 in their paper), confirming the accuracy of our numerical method. Throughout this work, the linear eigenvalue problem (3.14)–(3.15) was solved with $L=200$ and $N=8000$ for all waves we tested.

Figure 7.

The most unstable normal mode for the elevation wave with $\unicode[STIX]{x1D6FC}=0.3$ (which corresponds to $c=1.3512$ ). The profiles are modified in scale in order to compare with figure 10 of Calvo & Akylas (Reference Calvo and Akylas2002).

Table 1.

Instability growth rate $\unicode[STIX]{x1D706}$ of elevation solitary waves for $c=1.36$ and $c=1.4$ . The value of the middle point and the energy of the wave are listed for reference.

We start to explore the stability characteristics of the elevation branch near other limiting points. In figure 2(b), away from the turning point $2B$ no instability can be detected until the stationary point $2C$ , beyond which the elevation waves become unstable again. Table 1 shows the instability growth rate of unstable elevation waves at $c=1.36$ and $c=1.4$ . The most striking result is the instability of the waves after the fourth turning point, whose profiles feature two big troughs connected by two small dimples (see $2j$ in figure 2 a). The structure of these waves is similar to the known stable elevation waves except for the local shape between two big troughs. However, the increase of the distance between the troughs weakens their coherence, therefore a ‘collision-like’ instability occurs. An example of this phenomenon is presented in figure 8, whereby the small perturbation accelerates the left trough and decelerates the right one. Figure 8(a) depicts the profiles of the perturbation, which is the most unstable normal mode with the corresponding growth rate $\unicode[STIX]{x1D706}=0.00614$ . We present in (b,c) the comparisons between the numerical simulation of the full Euler equations and the prediction of the linear theory at $t=600,750$ , which show great agreement. After $t=750$ , the nonlinear effects become significant and a quasi-elastic collision is observed where both waves survive and the primary effect of the collision is a rapid phase shift (see the wave trajectories in 8 d and typical wave profiles in 8 e). The computation shown in figure 8 was carried out with $L=300$ , $N=8000$ and $\unicode[STIX]{x0394}t=4\times 10^{-4}$ .

Figure 8.

Instability of an elevation wave which consists of two depression waves connected by two small dimples. The steady solution, which corresponds to the point $2A$ in figure 2(a) ( $\unicode[STIX]{x1D702}(0)=0.0425$ and $c=1.39$ ), was disturbed by its most unstable mode ( $\unicode[STIX]{x1D706}=0.00614$ ) shown in (a) for $\widetilde{\unicode[STIX]{x1D702}}$ (top) and $\widetilde{\unicode[STIX]{x1D719}}$ (bottom). (b,c) Comparison between the linear theory (dashed line) and the nonlinear evolution code (solid line) at $t=600$ and $t=750$ respectively. (d) Time-dependent solution in the $\unicode[STIX]{x1D709}-t$ plane. (e) Four free-surface profiles are shown from top to bottom at times $t=750,1000,2000,8000$ (solid line), compared with the unperturbed steady solution (dotted line). The solution is shown in a frame of reference moving to the right with the speed of the undisturbed wave.

Lastly, we consider the stability of asymmetric gravity–capillary solitary waves. Following Wang et al. (Reference Wang, Vanden-Broeck and Milewski2014), the initial guess for finding an asymmetric solitary wave was obtained by superposing a depression wave ( $1a$ in figure 1 a) and a stable elevation one ( $2d$ in figure 2 a), both of which translate at $c=1.4$ . With this initial guess, the Newton’s algorithm converged to the solution 9f in figure 9. It is natural to ask whether the obtained asymmetric solution is stable since its constituent parts are both stable.

Using 9f as the starting point, the bifurcation curve can be followed, which is shown in figure 9. The spectral stability analysis was then carried out using the method described above, and the growth rate is presented in table 2. We conclude from table 2 that all asymmetric waves on this branch are unstable, even when the wave profiles resemble two stable solitary waves placed side by side. It is noted that the investigation of the global bifurcation of asymmetric gravity–capillary solitary waves is beyond the scope of this paper, and the interested reader is referred to Wang et al. (Reference Wang, Vanden-Broeck and Milewski2014).

Table 2.

Instability growth rate $\unicode[STIX]{x1D706}$ of asymmetric solitary waves for various values of the speed $c$ . The energy of the wave is listed for reference.

Figure 9.

Speed–energy bifurcation diagram of a branch of asymmetric solitary waves, together with typical wave profiles for $c=1.4$ (9c and 9f) and $c=1.39$ (9a and 9h). The branch arises from an asymmetric wave with the initial guess being composed of a stable elevation wave and a depression one propagating at $c=1.4$ .

Figure 10.

Time evolution of the unstable mode for an asymmetric solitary wave (9c in figure 9). All profiles are shown in a frame of reference moving to the right with $c=1.4$ . (a) The profiles of the most unstable normal mode. (b,c) Comparison between the linear theory (dashed line) and the nonlinear computation (solid line) at $t=125$ and $t=160$ , respectively. (d) The evolution of the perturbation using the nonlinear evolution code. (e) Wave profiles at $t=160,200,500,1200$ from top to bottom, with dotted line corresponding to the unperturbed steady solution.

Figures 10 and 11 show two numerical experiments of asymmetric solitary waves (9c and 9f respectively) perturbed with their most unstable normal modes, which were performed with $L=300$ , $N=8000$ and $\unicode[STIX]{x0394}t=4\times 10^{-4}$ . Figure 10(b–d) show the nonlinear evolution of the initial disturbance whose profiles are given in figure 10(a). Figure 10(c) strongly suggests that the linear theory captures the development of the unstable mode until $t\approx 160$ , a value of order $1/\unicode[STIX]{x1D706}$ (here $\unicode[STIX]{x1D706}=0.02608$ ). Snapshots in figure 10(e) demonstrate the long-time behaviour of the instability, which manifests itself by a visible change of the central troughs and the subsequent breakdown of the coherent structure. Figure 11 depicts the instability of the wave 9f. The decrease of the amplitude of the elevation component increases its speed, while the increase of the amplitude of the depression component reduces its speed. Hence the small difference in speed between two components leads to changes of the connecting portion during the interaction. The comparisons between the nonlinear simulation based on (2.20) and the linear theory are made in figure 11(b,c), which show remarkable agreement. In contrast to the instability shown in figure 8, the elevation component is totally destroyed in the long-time dynamics.

Figure 11.

Time-dependent computation of an asymmetric solitary waves (9f in figure 9) subject to the unstable normal mode. All profiles are shown in a frame of reference moving to the right with $c=1.4$ . (a) Profiles of the most unstable normal mode. (b,c) Comparison between the linear theory (dashed line) and the nonlinear computation (solid line) at $t=800$ and $t=1000$ respectively. (d) From top to bottom, the snapshots show the evolution of the free surface at time $t=1000,1500,2750,5000$ (solid line), and the unperturbed profile for reference.