4.1. The first limit point

The local behaviour in the vicinity of the first limit point $c_1$ is consistent with that observed in inviscid hydroelastic and capillary–gravity wave systems (Akylas Reference Akylas1984; Calvo & Akylas Reference Calvo and Akylas2002; Părău & Dias Reference Plotnikov and Toland2002). Specifically, as the load speed approaches the critical threshold where the group and phase velocities are matched, a wavepacket solitary-wave structure emerges, a regime effectively captured by weakly nonlinear asymptotic theory. To elucidate the underlying mechanism, we first examine the idealised inviscid limit ( $\mu =0$ ), focusing on the fundamental balance between nonlinearity and dispersion that precipitates the first fold. We adopt the standard asymptotic scaling

(4.1) \begin{equation} c=c_m-\epsilon ^2\alpha,\quad \zeta \sim \phi \sim O(\epsilon ),\quad P_0\sim O(\epsilon ^2), \end{equation}

where $c_m=\min c_p(k)$ denotes the minimum of the linear phase speed, $\alpha$ represents the speed detuning, and $\epsilon \ll 1$ is a small parameter. Following the perturbation procedure for the nonlinear Schrödinger (NLS) equation, we introduce the slow variables $X=\epsilon x$ and $\tau =\epsilon ^2t$ . Invoking the solvability condition at $O(\epsilon ^3)$ yields an fNLS equation for the complex wave envelope $A(X,\tau )$ . In the steady frame ( $\partial _\tau =0$ ), the envelope is governed by the stationary fNLS equation

(4.2) \begin{equation} -k_m\alpha A+\frac {\omega _{kk}}{2}\frac {\partial ^2 A}{\partial X^2}+\gamma\, |A|^2A=\frac {\tanh (k_mH)}{2c_m}P_0\,\delta (X), \end{equation}

where $\delta (X)$ is the Dirac delta function, and the interface displacement can be reconstructed as $\zeta \sim 2\,\textrm{Re} (A\,{\rm e}^{{i}k_mx} )$ , with $k_m$ the wavenumber at the minimum phase speed $c_m$ . The coefficients $\omega _{kk}$ and $\gamma$ depend on $k_m$ and the water depth $H$ (see Appendix B). Within the prescribed non-dimensionalisation (2.7), the condition $H\lt 233$ ensures that $\gamma$ and $\omega _{kk}$ possess the same sign ( $\gamma \omega _{kk}\gt 0$ ), identifying (4.2) as a focusing regime that admits bright envelope soliton solutions (Guyenne & Părău Reference Guyenne and Părău2014). For such localised states, the fNLS theory yields a leading-order relationship between the first limit-point speed and the load magnitude, derived by applying the analytical framework of Calvo & Akylas (Reference Calvo and Akylas2002):

(4.3) \begin{equation} c_1^{(u)}\approx c_m-\sqrt {\frac {\gamma }{\omega _{kk}}}\,\frac {\tanh (k_mH)}{2k_mc_m}\,P_0, \end{equation}

while the peak amplitude along the bifurcation branch is approximated by

(4.4) \begin{equation} \max {|\zeta |}\approx 2\sqrt {\frac {k_m}{\gamma }\left [(c_m-c)\pm \sqrt {(c_m-c)^2-(c_m-c_1^{(u)})^2}\right ]}\!. \end{equation}

Here, the superscript in $c_1^{(u)}$ distinguishes the prediction derived from the inviscid theory.

To elucidate the relative insensitivity of the first limit point $c_1$ to the damping coefficient $\mu$ , we examine the numerical results obtained from the fully nonlinear Euler equations. As illustrated in figure 8, the data points corresponding to a broad range of damping strengths ( $\mu \in [0.01,0.2]$ ) collapse remarkably onto the dashed line representing the idealised inviscid case ( $\mu =0$ ). This robust collapse suggests that the first transition at $c_1$ is governed by a nonlinear–dispersive balance that remains largely decoupled from dissipative effects. Consequently, the threshold for the first bifurcation is effectively determined by the inviscid limit.

Furthermore, the linear $c_1^{(u)}(P_0)$ relationship derived from the fNLS framework (4.3) is superimposed on the numerical data in figure 8. The fNLS prediction serves as an asymptotic tangent to the numerical curve as the load magnitude vanishes ( $P_0\to 0$ ). While the weakly nonlinear theory recovers the initial slope, the numerical data exhibit a gradual upward departure at larger $P_0$ , indicating that the cubic NLS regime is restricted to the immediate vicinity of the bifurcation inception. Collectively, these observations reinforce the conclusion that the inception of the first limit point is a phenomenon primarily governed by the leading-order balance between nonlinearity and dispersion under localised forcing. In this regime, dissipative effects exert only a subsidiary influence on the displacement of the bifurcation threshold.

Figure 8.

Variation of the first limit-point speed $c_1$ with load magnitude $P_0$ at water depth $H=1$ . The red solid line represents the leading-order asymptotic prediction $c_1^{(u)}$ derived from the fNLS theory (4.3), while the blue solid line corresponds to the refined asymptotic estimate (5.23). Results from the fully nonlinear Euler equations are superimposed for comparison: the dashed line denotes the inviscid case ( $\mu =0$ ), while discrete symbols represent various damping strengths (circles for $\mu =0.01$ ; squares for $\mu =0.05$ ; triangles for $\mu =0.10$ ; diamonds for $\mu =0.20$ ). The damped numerical data collapse onto a nearly identical curve, with the fNLS prediction (4.3) serving as the asymptotic tangent in the vanishing-load limit ( $P_0\to 0$ ).

4.2. The second limit point

Numerical bifurcation analysis reveals that the second limit-point speed $c_2$ is parametrically sensitive to both the damping coefficient $\mu$ and the load magnitude $P_0$ . Unlike the first fold $c_1$ , which is governed by the balance between nonlinearity and dispersion encapsulated in the fNLS description (4.3), the dynamics near $c_2$ is dictated by the balance between external power input and systemic dissipation. In this large-amplitude regime, standard weakly nonlinear approximations fail, and a closed-form expression analogous to (4.3) remains elusive. Nevertheless, a global energy balance provides a robust physical constraint, enabling a systematic estimation of $c_2$ based on power-balance requirements.

For a travelling wave moving with load speed $c$ , the total mechanical energy must remain stationary in the co-moving frame. The dimensionless total energy $E$ integrates the kinetic, gravitational potential and elastic strain contributions over the fluid domain and the ice cover:

(4.5) \begin{equation} E=\frac 12\int _{\mathbb{R}}\int _{-H}^{\zeta }|\boldsymbol{\nabla }\phi |^2\,\mathrm{d}y\,\mathrm{d}x+\frac 12\int _{\mathbb{R}}\left (\zeta ^2+\kappa ^2\sqrt {1+\zeta _x^2}\right )\mathrm{d}x. \end{equation}

By evaluating the time derivative of the total energy $E$ and invoking the divergence theorem, the energy budget is reduced to a fundamental balance between the rate of work performed by the moving load and the rate of energy dissipation due to damping:

(4.6) \begin{align} \frac {\mathrm{d}E}{\mathrm{d}t}&=\int _{\mathbb{R}}\int _{-H}^{\zeta }\boldsymbol{\nabla }\phi \boldsymbol{\cdot }\boldsymbol{\nabla }\phi _t\,\mathrm{d}y\,\mathrm{d}x+\int _{\mathbb{R}}\zeta _t\left [\frac 12|\boldsymbol{\nabla }\phi |^2+\zeta +\kappa _{ss}+\frac 12\kappa ^3\right ]_{y=\zeta }\mathrm{d}x \notag \\[5pt] &=\int _{\mathbb{R}}\sqrt {1+\zeta _x^2}\left [\phi _t\frac {\partial \phi }{\partial \boldsymbol{n}}\right ]_{y=\zeta }\mathrm{d}x-\int _{\mathbb{R}}\zeta _t\left [\phi _t+\mu \zeta _t+p(x-ct)\right ]_{y=\zeta } \mathrm{d}x \\[5pt] &=-\int _{\mathbb{R}}\big [\mu \zeta _t^2+p(x-ct)\zeta _t\big ]\mathrm{d}x. \notag \end{align}

For a steady-state profile of the form $\zeta =\zeta (x-ct)$ , the energy remains stationary, and this global balance yields the integral identity

(4.7) \begin{equation} \mu c\int _{\mathbb{R}}\zeta _x^2\,\mathrm{d}x=\int _{\mathbb{R}}p(x)\,\zeta _x\,\mathrm{d}x. \end{equation}

Physically, the integral identity (4.7) dictates that the viscous dissipation rate within the ice sheet (left-hand side) must be precisely balanced by the power performed by the moving load (right-hand side). In the limit of a concentrated point-load $p(x)=P_0\,\delta (x)$ , this global balance reduces to a pointwise constraint on the surface slope at the load’s application point:

(4.8) \begin{equation} \zeta _x(0)=c\,\frac {\mu }{P_0}\int _{\mathbb{R}}\zeta _x^2\,\mathrm{d}x\gt 0. \end{equation}

This condition yields immediate physical insights into the wave shape. Since the dissipative integral is strictly positive, the surface slope at the load position must satisfy $\zeta _x(0)\gt 0$ . In the co-moving frame, a positive slope at $x=0$ implies that the load is positioned on the rising (rear) face of the wave, such that the maximum deflection (the trough) lags behind the load. This spatial phase lag is a fundamental dissipative effect, consistent with the experimental observations of Takizawa (Reference Takizawa1985).

Following the decomposition approach of Cho & Akylas (Reference Cho and Akylas2009), we partition the state II steady-state profile $\zeta (x)$ into a dominant free solitary-wave component $\bar {\zeta }$ and a small forced correction $\widetilde {\zeta }$ :

(4.9) \begin{equation} \zeta (x)=\bar {\zeta }(x+\theta )+\widetilde {\zeta }(x), \end{equation}

where $\theta$ represents the unknown spatial lag of the solitary-wave extremum relative to the load. To establish a solvable leading-order balance, we adopt the distinguished scaling

(4.10) \begin{equation} \bar {\zeta }\sim O(1),\quad \widetilde {\zeta }\sim O(\epsilon ),\quad p(x)\sim O(\epsilon ),\quad \mu \sim O(\epsilon ),\quad \theta =\theta ^{(0)}+\epsilon \theta ^{(1)}, \end{equation}

where $\epsilon \ll 1$ characterises the relative smallness of both the dissipative and forcing effects compared to the wave amplitude. At leading order in $\epsilon$ , the energy-balance identity (4.7) is dominated by the free-wave component, reducing to

(4.11) \begin{equation} \mu c\int _{\mathbb{R}}\bar {\zeta }_x(x)^2\,\mathrm{d}x=\int _{\mathbb{R}}p(x)\,\bar {\zeta }_x(x+\theta ^{(0)})\,\mathrm{d}x. \end{equation}

For a concentrated point-load $p(x)=P_0\delta (x)$ , the right-hand side evaluates to $P_0\,\bar {\zeta }_x(\theta ^{(0)})$ , yielding the explicit relationship

(4.12) \begin{equation} \bar {\zeta }_x(\theta ^{(0)}) = c\,\frac {\mu }{P_0}\int _{\mathbb{R}}\bar {\zeta }_x(x)^2 \,\mathrm{d}x. \end{equation}

Figure 9.

Spatial lag characteristics of the large-amplitude (state II) branch. ( $a$ ) Graphical determination of (4.12) for $\mu /P_0=2$ and $c=0.75$ . The intersections of the free-wave slope $\bar {\zeta }_x$ (solid line) with the dissipative threshold $c(\mu /P_0)\int \bar {\zeta }_x^2\,\mathrm{d}x$ (dashed line) identify two admissible equilibrium positions, $\theta _1$ and $\theta _2$ . (b) Bifurcation of the lag distance $\theta$ as a function of the load speed $c$ . The solid line represents the asymptotic prediction obtained from the roots of (4.12), while discrete circles denote results from the full Euler equations: orange for $\mu =0.2$ , $P_0=0.1$ , and blue for $\mu =0.02$ , $P_0=0.01$ .

Equation (4.12) implicitly constrains the leading-order lag $\theta ^{(0)}$ as a function of the damping-to-load ratio $\mu /P_0$ and the propagation speed $c$ . A graphical interpretation of this relationship is provided in figure 9( $a$ ). Given that the solitary-wave slope $\bar {\zeta }_x$ is an antisymmetric, bipolar function, as characterised by the alternating positive and negative extrema in figure 9( $a$ ), the dissipative threshold $c(\mu /P_0)\int \bar {\zeta }_x^2 \,\mathrm{d}x$ (represented by the horizontal dashed line) may intersect the $\bar {\zeta }_x$ curve at multiple points. For the representative parameters $\mu /P_0=2$ and $c=0.75$ , two distinct roots, $\theta _1$ and $\theta _2$ , exist, signifying the bifurcation into a multi-valued regime in the phase-lag dynamics. This correspondence between the lag and amplitude bifurcations provides a mechanistic basis for the dual-state behaviour observed in state II. It confirms that the spatial lag is not merely a passive translation, but a fundamental dynamical signature of the forced–dissipative interaction between the moving load and the hydroelastic wave.

The multiplicity of roots in the transcendental lag (4.12) facilitates the asymptotic prediction of the second limit point. At the critical speed $c=c_2$ , the two distinct lag solutions coalesce into a degenerate double root, marking a saddle-node bifurcation in the phase-lag dynamics. Geometrically, this bifurcation occurs when the dissipative threshold on the right-hand side of (4.12) becomes tangent to the extremum of the bipolar slope function $\bar {\zeta }_x$ . We denote this leading-order asymptotic prediction as $c_2^{(0)}$ . As illustrated in figure 9( $b$ ) for a fixed ratio $\mu /P_0=2$ , the second limit speed $c_2$ extracted from the fully nonlinear Euler equations converges towards the theoretical value $c_2^{(0)}$ in the vanishing-load limit ( $P_0\to 0$ ). This agreement validates the distinguished scaling adopted in (4.10), and confirms that the energy-balance framework captures the fundamental physics governing the switching between steady states.

Because the leading-order prediction $c_2^{(0)}$ is linked to the ratio $\mu /P_0$ via a double-root (tangency) condition, imposing this criterion on (4.12) yields the implicit relationship

(4.13) \begin{equation} \frac {\mu }{P_0}=\frac {\max \left \{\bar {\zeta }_x(x;c_2^{(0)})\right \}}{c_2^{(0)}\displaystyle \int _{\mathbb{R}}\bar {\zeta }_x(x;c_2^{(0)})^2\,\mathrm{d}x}. \end{equation}

This equation provides a robust estimate for the bifurcation threshold. For the case of unit dimensionless depth ( $H=1$ ), figure 10( $a$ ) displays the functional dependence of $c_2^{(0)}$ on the damping-to-load ratio $\mu /P_0$ . The asymptotic analysis identifies a critical threshold $(\mu /P_0)_{{cr}}\approx 1.856$ below which the second fold vanishes and the S-shaped bifurcation structure is lost. This prediction is in excellent agreement with the fully nonlinear results (discrete symbols) and the limit-point data previously reported in figure 7. As $H$ increases, the threshold $(\mu /P_0)_{{cr}}$ decreases monotonically, converging towards a deep-water asymptotic plateau of approximately $0.3180$ . This trend indicates that the sensitivity of the hysteretic response to water depth is most pronounced in the shallow-water regime, where small variations in $H$ result in significant shifts in the critical ratio.

Figure 10.

Characteristics of the second limit point and the critical damping-to-load threshold. ( $a$ ) Second limit-point speed $c_2^{(0)}$ as a function of the ratio $\mu /P_0$ at unit depth ( $H=1$ ). The solid curve represents the asymptotic prediction from (4.13), while discrete symbols denote fully nonlinear results for $P_0=0.1$ (blue squares) and $P_0=0.01$ (orange circles). The vertical dashed line indicates the critical ratio $(\mu /P_0)_{{cr}}\approx 1.856$ , below which the second fold is annihilated. ( $b$ ) Variation of the critical ratio $(\mu /P_0)_{{cr}}$ with dimensionless water depth $H$ . The horizontal dashed line marks the deep-water asymptotic plateau $(\mu /P_0)_{{cr}}\approx 0.3180$ .

5.1. The forced–dissipative Whitham equation

Following Dinvay et al. (Reference Dinvay, Kalisch and Părău2019), the unidirectional propagation of weakly nonlinear hydroelastic waves is governed by a Whitham-type equation. Within the dimensionless framework established in § 2, this evolution of the surface displacement $\zeta (x,t)$ is expressed as

(5.1) \begin{equation} \zeta _t+\mathcal{K}[\zeta _x]+\sigma \zeta \zeta _x=0, \end{equation}

where the pseudo-differential operator $\mathcal{K}$ incorporates the full linear dispersion relation via the Fourier symbol

(5.2) \begin{equation} \mathcal{K}\big [{\rm e}^{{i}kx}\big ]=\sqrt {\left (\frac 1k+k^3\right )\tanh (kH)}\,{\rm e}^{{i}kx}. \end{equation}

The nonlinear coefficient is defined as $\sigma =3/(2\sqrt {H})$ . Whitham-type equations are well known to support solitary waves whose bifurcation branches exhibit stability exchange (see e.g. Akers & Milewski Reference Akers and Milewski2008; Akylas & Cho Reference Akylas and Cho2008). To account for external forcing and energy dissipation, we extend this formulation by incorporating a localised pressure distribution and a Burgers-type dissipative term, following the methodology of Cho & Akylas (Reference Cho and Akylas2009). In the co-moving frame $x\to x-ct$ , the forced–dissipative Whitham equation reads

(5.3) \begin{equation} \zeta _t-c\zeta _x+\mathcal{K}[\zeta _x]+\sigma \zeta \zeta _x-\nu \zeta _{xx}=-\frac 12p_x, \end{equation}

where $\nu$ denotes a phenomenological dissipation coefficient. This parameter is not an independent variable; its relationship to the damping coefficient $\mu$ in the Euler system can be established by enforcing energetic consistency between the two frameworks. Specifically, by multiplying (5.3) by $\zeta$ and integrating over the entire domain, while invoking the self-adjointness of $\mathcal{K}$ and the asymptotic decay of the solitary-wave-like profile, we arrive at the global integral balance

(5.4) \begin{equation} 2\nu \int _{\mathbb{R}}\zeta _x^2\,\mathrm{d}x=\int _{\mathbb{R}}p(x)\,\zeta _x\,\mathrm{d}x. \end{equation}

Figure 11.

Comparative bifurcation and stability analysis across the fully nonlinear Euler equations, the forced–dissipative Whitham model, and the reduced dynamical system (RDS). ( $a$ ) Bifurcation diagrams depicting the maximum deflection $\max |\zeta |$ versus the load speed $c$ for the Whitham (blue, $y=0$ plane) and Euler (red, $y=1$ plane) frameworks. Solid and dashed lines denote stable and unstable equilibrium branches, respectively, while dotted lines represent the corresponding free solitary-wave loci. Discrete black markers, where circles and crosses denote stable and unstable fixed points of the reduced dynamical system, respectively, demonstrate an excellent overlap with the Whitham bifurcation curve. ( $b$ ) Linear stability of the Whitham model. The curve tracks the maximum real part of the eigenvalues, $\textrm{Re}(\lambda )$ , along the bifurcation branch. The transition from a stable (blue) to an unstable (orange) state occurs precisely at the limit points, where the dominant eigenvalue crosses the imaginary axis, $\textrm{Re}(\lambda )=0$ . Parameters are $\mu =0.2$ , $H=1$ and $P_0=0.1$ .

This relation represents the fundamental equilibrium between dissipative loss and the power from the moving load. By comparing the Whitham-level energy budget (5.4) with the global balance derived for the full Euler system (4.7), we identify the velocity-dependent mapping $\nu =c\mu /2$ . This scaling is physically essential to restore the reflectional symmetry of the response in the static limit. In the full Euler framework, a static load ( $c=0$ ) yields a symmetric profile, as the dissipation term $\mu \zeta _t$ vanishes identically. Conversely, a Whitham-type model with a constant $\nu$ would erroneously predict an asymmetric profile and a spatial lag even as $c\to 0$ . By adopting $\nu =c\mu /2$ , the dissipative effect correctly vanishes in the static limit, consistent with the underlying physics. This reduced model qualitatively recovers the S-shaped bifurcation structure observed in the fully nonlinear system (blue lines in figure 11 $a$ ). Crucially, it also reproduces the exchange of stability at the limit points (figure 11 $b$ ), a defining dynamical feature that governs the jump phenomena in the full Euler framework. It should be emphasised that the present formulation is not intended to account for the spatial attenuation of freely propagating waves. Instead, in the unforced, inviscid limit ( $\nu =p=0$ ), the model reverts to the monotonic speed–amplitude relationship characteristic of free solitary waves (dotted lines in figure 11 $a$ ), thereby confirming its efficacy in capturing the essential nonlinear–dispersive dynamics of the hydroelastic wave problem. It is noteworthy that the application of a non-constant dissipation coefficient has precedent in the study of other water-wave problems; for instance, Longuet-Higgins (Reference Longuet-Higgins1992) and Tian et al. (Reference Tian, Perlin and Choi2010) employed an eddy viscosity within a quasi-potential framework to model the concentrated energy dissipation associated with wave breaking in surface gravity waves.

To evaluate the efficacy of the Whitham-type model in the co-moving frame ( $\zeta _t=0$ ), we apply an asymptotic decomposition to the state II response, $\zeta =\bar {\zeta }(x+\theta )+\widetilde {\zeta }(x)$ . This expansion is consistent with the distinguished scalings established in (4.10), where the effective dissipation is scaled as $\nu \sim O(\epsilon )$ . At leading order, the governing equation (5.3) reduces to the nonlinear–dispersive balance for a free, undamped Whitham solitary wave:

(5.5) \begin{equation} -c\bar {\zeta }_x+\mathcal{K}[\bar {\zeta }_x]+\sigma \bar {\zeta }\bar {\zeta }_x=0. \end{equation}

At $O(\epsilon )$ , the first-order correction $\widetilde {\zeta }(x)$ is governed by the linearised inhomogeneous equation

(5.6) \begin{equation} \mathcal{L}[\widetilde {\zeta }]\equiv \left \{-c\,\partial _x+\mathcal{K}[\partial _x]+\sigma\, \partial _x(\bar {\zeta }{\cdot })\right \}[\widetilde {\zeta }]=\nu\, \bar {\zeta }_{xx}(x+\theta ^{(0)})-\frac 12p_x, \end{equation}

where the linear operator $\mathcal{L}$ is evaluated at the shifted leading-order solitary wave $\bar {\zeta }(x+\theta ^{(0)})$ . The adjoint operator $\mathcal{L}^\dagger$ is given by

(5.7) \begin{equation} \mathcal{L}^\dagger =c\,\partial _x-\mathcal{K}[\partial _x]-\sigma\, \bar {\zeta }(x+\theta ^{(0)})\,\partial _x. \end{equation}

It follows from (5.5) that the free solitary-wave profile $\bar {\zeta }(x+\theta ^{(0)})$ spans the null space of the adjoint operator, i.e. $\mathcal{L}^\dagger [\bar {\zeta }]=0$ , given the asymptotic decay of the solution as $x\to \pm \infty$ . Consequently, the Fredholm alternative dictates that a bounded solution for $\widetilde {\zeta }$ exists if and only if the inhomogeneous forcing is orthogonal to the kernel of $\mathcal{L}^\dagger$ . This solvability condition requires the following inner product to vanish:

(5.8) \begin{equation} \int _{\mathbb{R}}\left [\nu\, \bar {\zeta }_{xx}(x+\theta ^{(0)})-\frac 12p_x(x)\right ]\bar {\zeta }(x+\theta ^{(0)})\,\mathrm{d}x=0. \end{equation}

Upon integration by parts and simplification, this orthogonality condition is transformed into the global energy-balance relation

(5.9) \begin{equation} 2\nu \int _{\mathbb{R}}\bar {\zeta }_x^2\,\mathrm{d}x=\int _{\mathbb{R}} p(x)\,\bar {\zeta }_x(x+\theta ^{(0)})\,\mathrm{d}x. \end{equation}

This result demonstrates that evaluating the global power balance (5.4) with a shifted solitary-wave ansatz is equivalent to satisfying the solvability condition of the linearised problem. The adjoint-based derivation, however, highlights the universality of the lag-estimation relation (5.9) for a broad class of forced–dissipative dispersive systems. Analogous relations have been reported for the fifth-order Korteweg–de Vries (KdV) equation (Cho & Akylas Reference Cho and Akylas2009) and the Akers–Milewski model (Akers & Milewski Reference Akers and Milewski2009), confirming that this integral balance is a generic attribute of localised nonlinear waves interacting with moving disturbances. Furthermore, this formal consistency provides a secondary justification for the velocity-dependent scaling $\nu =c\mu /2$ . The capacity of the Whitham-type model to recover the same fundamental energy budget as the fully nonlinear Euler theory validates its role as a robust reduced framework for exploring the dynamical stability and the underlying phase-space structure of the system.

5.2. A low-dimensional dynamical system

To elucidate the dynamical bridge between the hydroelastic response and the theory of nonlinear resonators, we reduce the forced–dissipative Whitham equation (5.3) to a two-dimensional system via a Galerkin-style projection. We adopt a modulated solitary-wave ansatz of the form

(5.10) \begin{equation} \zeta =a(t)\,\bar {\zeta }(x+\theta (t);c), \end{equation}

where the time-dependent variables $a(t)$ and $\theta (t)$ represent the instantaneous amplitude and spatial phase lag, respectively. Substituting this ansatz into (5.3) and invoking the underlying free-wave relationship (5.5) to eliminate the leading-order dispersive terms, we identify the residual dynamics governing the evolution of the wave modulation:

(5.11) \begin{equation} \bar {\zeta }\,\frac {\mathrm{d}a}{\mathrm{d}t}+a\bar {\zeta }_x\,\frac {\mathrm{d}\theta }{\mathrm{d}t}+\sigma a(a-1) \bar {\zeta }\bar {\zeta }_x-a\nu \bar {\zeta }_{xx}=-\frac 12p_x. \end{equation}

Following the methodology of Cousins & Sapsis (Reference Cousins and Sapsis2015), we employ a dual projection to decouple the evolution of the amplitude $a(t)$ and the spatial lag $\theta (t)$ . First, projecting the residual onto the wave profile $\bar {\zeta }(x+\theta )$ yields the amplitude evolution (representing a global energy-like balance):

(5.12) \begin{equation} \left (\int _{\mathbb{R}}\bar {\zeta }^2\,\mathrm{d}x\right )\frac {\mathrm{d}a}{\mathrm{d}t}=-\nu a \int _{\mathbb{R}}\bar {\zeta }_x^2\, \mathrm{d}x +\int _ {\mathbb{R}}\frac 12p(x)\,\bar {\zeta }_x(x+\theta )\,\mathrm{d}x. \end{equation}

Second, to isolate the phase dynamics, we project the residual onto the potential-like test function $\bar {\phi }(x+\theta )= \int _0^{x+\theta }\bar {\zeta }(\xi )\,\text{d}\xi$ . This choice is strategically motivated: the function $\bar {\phi }$ represents the spatial primitive of the elevation, and as demonstrated in Appendix C, this specific projection preserves the Hamiltonian structure of the underlying system in the dissipationless limit. Upon integration by parts and invoking the decaying boundary conditions, the evolution of the spatial lag is governed by

(5.13) \begin{equation} \left (-a\int _{\mathbb{R}}\bar {\zeta }^2\,\mathrm{d}x\right )\frac {\mathrm{d}\theta }{\mathrm{d} t}=\frac 12 \sigma a(a-1) \int _{\mathbb{R}} \bar {\zeta }^3\,\mathrm{d}x+\int _{\mathbb{R}}\frac 12p(x)\,\bar {\zeta }(x+\theta )\,\mathrm{d}x. \end{equation}

For a concentrated load $p(x)=P_0\,\delta (x)$ , the reduced dynamical system is expressed as

(5.14) \begin{align} \frac {\mathrm{d}a}{\mathrm{d}t}&=-\nu r_1a+\frac {P_0}{2}r_2\,\bar {\zeta }_x(\theta ), \\[-12pt]\nonumber \end{align}

(5.15) \begin{align} \frac {\mathrm{d}\theta }{\mathrm{d}t}&=-\frac 12\sigma r_3(a-1)-\frac {P_0}{2 a}r_2\, \bar {\zeta }(\theta ), \end{align}

where the geometric coefficients $r_j$ ( $j=1,2,3$ ) are defined by the integral moments of the basis function as

(5.16) \begin{equation} r_1=\frac {\int _{\mathbb{R}}\bar {\zeta }_x^2\,\mathrm{d}x}{\int _{\mathbb{R}}\bar {\zeta }^2\,\mathrm{d}x},\quad r_2=\frac {1}{\int _{\mathbb{R}}\bar {\zeta }^2\,\mathrm{d}x},\quad r_3=\frac {\int _{\mathbb{R}}\bar {\zeta }^3\,\mathrm{d}x}{\int _{\mathbb{R}}\bar {\zeta }^2\,\mathrm{d}x}. \end{equation}

It should be noted that the basis function $\bar {\zeta }$ depends parametrically on the load speed $c$ . The resulting system is structurally analogous to the averaged equations of a forced–damped Duffing oscillator (Nayfeh & Mook Reference Nayfeh and Mook1995). This structural identity reveals that the spatial lag $\theta$ is the hydrodynamic analogue of the temporal phase shift in classical nonlinear resonators. Consequently, the reduced-order system provides a robust dynamical framework, characterising the observed jump phenomena as a generic manifestation of a folded equilibrium manifold in the system’s phase space.

We now examine the fixed points $(a_0,\theta _0)$ of the reduced dynamical system, and their associated stability characteristics. By imposing the steady-state conditions on (5.14)–(5.15), and eliminating $a_0$ , we obtain the implicit equation for the equilibrium manifold $G(\theta ,c)=0$ :

(5.17) \begin{equation} G(\theta ,c)\equiv \frac 12\frac {r_3\sigma }{ r_1\nu }\,\bar {\zeta }_x(\theta )\left [\frac {r_2P_0}{2r_1\nu }\,\bar {\zeta }_x(\theta )-1 \right ]+\bar {\zeta }(\theta )=0. \end{equation}

The bifurcation curve defined by the fixed points of (5.17) is superimposed on the numerical results from the Euler and Whitham frameworks in figure 11( $a$ ). The reduced-order system successfully recovers the characteristic S-shaped bifurcation structure, and exhibits excellent quantitative agreement with the parent forced–dissipative Whitham model. This overlap confirms that the projection onto the amplitude and lag coordinates captures the system’s essential degrees of freedom, thereby validating the dimensionality reduction and the chosen ansatz for the hydroelastic response.

The linear stability of these fixed points is determined by the spectrum of the Jacobian matrix $\boldsymbol{J}$ , defined as

(5.18) \begin{equation} \boldsymbol{J}=\begin{bmatrix} -\nu r_1 & \nu r_1a_0\dfrac {\bar {\zeta }_{xx}(\theta _0)}{\bar {\zeta }_x(\theta _0)}\\[11pt] \sigma r_3 \dfrac {1-2a_0}{2a_0} & -\nu r_1 \end{bmatrix}\!. \end{equation}

Notably, the trace is strictly negative, $\text{tr}(\boldsymbol{J})=-2\nu r_1\lt 0$ , which constrains the equilibria to be either stable nodes/foci or unstable saddle points. This condition precludes the existence of Hopf bifurcations, ensuring that the exchange of stability occurs exclusively via a zero-eigenvalue crossing at the limit points. The determinant of the Jacobian is found to be

(5.19) \begin{equation} \det (\boldsymbol{J})=\nu ^2r_1^2+\frac 12\nu \sigma r_1r_3\frac {\bar {\zeta }_{xx}(\theta _0)}{\bar {\zeta }_x(\theta _0)}(2a_0-1). \end{equation}

By differentiating the equilibrium relation (5.17), we establish the fundamental identity

(5.20) \begin{equation} \det (\boldsymbol{J})=\frac {P_0\nu r_1r_2}{2a_0}\frac {\partial G}{\partial \theta }. \end{equation}

The exchange of stability requires $\det (\boldsymbol{J})=0$ , which, according to the aforementioned identity, coincides with the criticality condition $G_{\theta }=0$ . By invoking the implicit function theorem, the slope of the equilibrium branch is given by $\mathrm{d}c/\mathrm{d}\theta =-G_{\theta }/G_c$ . Consequently, provided that the manifold is transversal to the speed axis ( $G_c\ne 0$ ), the points of stability exchange correspond precisely to the turning points (or limit points) of the $\theta{-}c$ bifurcation curve, where the tangent becomes vertical. While a local bijective mapping between the spatial phase lag $\theta$ and the amplitude $a$ near the turning points remains a theoretical hypothesis, the resulting framework provides a highly plausible dynamical mechanism for the observed jump phenomena: the bifurcation-induced instability at these limit points is predicted to trigger a catastrophic transition towards a coexisting stable equilibrium (Strogatz Reference Strogatz2018). Strikingly, this theoretical prediction is strongly vindicated by numerical simulations, showing excellent agreement with both the fully nonlinear Euler computations (figure 4 $a$ ) and the forced–dissipative Whitham model (figure 11 $b$ ) across a wide spectrum of load magnitudes and damping coefficients (figure 7).

While the steady-state scaling in § 4 provided preliminary estimates for the limit points, the reduced dynamical system permits a comprehensive prediction of $c_1$ that transcends the small-amplitude regime. The fixed points $(a_0,\theta _0)$ satisfy the coupled algebraic relations

(5.21) \begin{align} a_0&=\frac {P_0r_2}{2\nu r_1}\,\bar {\zeta }_x(\theta _0), \\[-12pt]\nonumber \end{align}

(5.22) \begin{align} a_0&=\frac 12\pm \sqrt {\frac 14-\frac {r_2P_0}{r_3\sigma }\,\bar {\zeta }(\theta _0)}. \end{align}

The first limit point $c_1$ is characterised by the coalescence of steady states at small lag distances ( $\theta _0\ll 1$ ). By expanding the wave profile near its extremum, where $\bar {\zeta }_x\approx \bar {\zeta }_{xx}(0)\,\theta _0$ and $\bar {\zeta }\approx \bar {\zeta }(0)$ , the condition for a fold bifurcation requires the discriminant in (5.22) to vanish. This yields the leading-order asymptotic estimate for $c_1^{(0)}$ :

(5.23) \begin{equation} \frac {\sigma }{4} \frac {1}{P_0} = \frac {r_2}{r_3}\, \bar {\zeta }(0) = \frac {\bar {\zeta }(0;c_1^{(0)})}{\int _{\mathbb{R}} \bar {\zeta }(x;c_1^{(0)})^3\,\mathrm{d}x}. \end{equation}

This result offers two significant physical insights. First, the limit-point speed $c_1$ is primarily governed by the forcing magnitude $P_0$ and exhibits a notable dissipative immunity, remaining independent of $\nu$ at leading order. This theoretical independence is in excellent agreement with the fully nonlinear observations presented in figure 7( $a$ ). As demonstrated in figure 8, the present estimate (blue solid line) shows markedly superior agreement with the Euler numerical results (discrete symbols) compared to the linear prediction from fNLS theory (4.3) (red solid line). Second, at the onset of the jump, the response amplitude is approximately half that of the corresponding free solitary wave, $a_0^{(0)}=1/2$ . This half-amplitude relation remains robust in the fully nonlinear Euler framework, with the ratio consistently falling within a narrow range between 0.5 and 0.6. Finally, the validity of this asymptotic estimate is contingent upon the parameter constraint $|\sigma r_3/\nu r_1|\gg 1$ , representing the predominance of nonlinearity over dissipation at the bifurcation threshold.

The estimation of the second limit point $c_2$ hinges on the observation that at this stage of the bifurcation, the response has fully transitioned into a large-amplitude state. Near $c_2$ , the wave profile closely mirrors that of a free solitary wave, justifying the approximation $a_0=1+\delta a$ with $\delta a\ll 1$ . In this regime, the coupled fixed-point relations reduce to

(5.24) \begin{align} 1=\frac {P_0 r_2}{2\nu r_1}\,\bar {\zeta }_x(\theta _0^{(0)}),\quad \delta a=-\frac {r_2P_0}{2r_3\sigma }\,\bar {\zeta }(\theta _0^{(0)}). \end{align}

The fold bifurcation at $c_2$ occurs when the first equation in (5.24) admits a double root for the spatial lag. Given the antisymmetric, bipolar structure of $\bar {\zeta }_x$ , this criticality condition corresponds to the point where the external forcing reaches the extremal slope of the free-wave profile. This yields the leading-order estimate

(5.25) \begin{equation} \frac {2\nu }{P_0}=\frac {r_2}{r_1}\max \{\bar {\zeta }_x\}=\frac {\max \{\bar {\zeta }_x(x;c_2^{(0)})\}}{\int _{\mathbb{R}} \bar {\zeta }_x(x;c_2^{(0)})^2\,\mathrm{d}x}. \end{equation}

This result is of central theoretical importance as it unifies three distinct analytical paths: it recovers the global energy balance of the full Euler system (4.13), coincides with the solvability condition for the Whitham-type model (5.9), and identifies $c_2$ as the point of degeneracy where the reduced Jacobian becomes singular. Consistent with the fully nonlinear data in figure 7( $b$ ), (5.25) confirms that the location of the second limit point is primarily governed by the damping-to-load ratio. Finally, the validity of this estimate remains subject to the same consistency constraint as $c_1$ , namely $|\sigma r_3/\nu r_1|\gg 1$ .

The reduced dynamical system further facilitates the classification of stable equilibria, providing a dynamical basis for the observation in the fully nonlinear Euler simulations. The nature of a stable fixed point (i.e. whether it is a node or a focus) is determined by the sign of the discriminant ${\varDelta }= \text{tr}(\boldsymbol{J})^2-4\det (\boldsymbol{J})$ . By substituting the Jacobian elements and invoking the steady-state conditions, we obtain

(5.26) \begin{equation} {\varDelta }=-\frac {P_0r_2\sigma r_3}{a_0}\,\bar {\zeta }_{xx}(\theta _0)\,(2 a_0-1). \end{equation}

For hydroelastic depression waves, the third-order integral moment $\int _{\mathbb{R}}\bar {\zeta }^3\,\mathrm{d}x$ is typically negative, implying $r_3\lt 0$ . Given $r_2\gt 0$ and the positive nonlinear coefficient $\sigma$ , the sign of the discriminant simplifies to

(5.27) \begin{equation} \text{sign}({\varDelta })=\text{sign}\left \{\bar {\zeta }_{xx}(\theta _0)\,(2 a_0-1)\right \}\!. \end{equation}

This criterion, coupled with our asymptotic estimates for the limit points, reveals that both stable branches consist of stable foci. In the small-amplitude regime, where stable states lie on the branch below $c_1$ (implying $a_0\lt 1/2$ ) and $\theta _0\to 0^+$ , the sign of ${\varDelta }$ is determined by $-\bar {\zeta }_{xx}(0)$ . Since a depression wave is locally concave-up at its trough, it follows that ${\varDelta }\lt 0$ . Similarly, in the large-amplitude regime ( $a_0\gt 1/2$ ), the equilibrium lag $\theta _0$ typically resides in the region where $\bar {\zeta }_{xx}(\theta _0)\lt 0$ (as exemplified by $\theta _2$ in figure 9 a), again yielding ${\varDelta }\lt 0$ . Consequently, both the small- and large-amplitude stable solutions are classified as stable foci, implying that the system approaches the steady state via damped oscillations in the $(a,\theta )$ phase plane. This result provides an asymptotic dynamical explanation for the ‘spiralling’ settling behaviour observed in the time-dependent Euler trajectories (e.g. figure 5 $b$ ), where the wave amplitude and spatial lag undergo mutual adjustment during the transient relaxation process. This qualitative consistency between the reduced-order dynamics and the full governing equations thus completes the dynamical characterisation of the hydroelastic jump phenomenon.

The limit points $c_1$ and $c_2$ partition the dynamical landscape of the hydroelastic wave system into three distinct regimes. For subcritical load speeds ( $c\lt c_2$ ), the system admits a unique small-amplitude steady state. In the reduced phase plane, this attractor is identified as a stable focus (figure 12 $a$ ), corresponding to the quasi-linear response observed in the fully nonlinear Euler simulations (see the brown dashed line in figure 3 a). Conversely, for speeds exceeding the first limit point ( $c\gt c_1$ ), the system converges towards a unique large-amplitude steady state. This state is likewise classified as a stable focus (figure 12 $c$ ), and manifests the characteristic profile of a free solitary wave (see the green dashed line in figure 3 a).

In the intermediate speed range ( $c_2\lt c\lt c_1$ ), the system exhibits bistability (figure 12 $b$ ), mirroring the coexistence of three equilibrium branches observed in the fully nonlinear Euler framework (figure 3). Within this regime, the small-amplitude state (i) and the large-amplitude state (iii) are both characterised as stable foci (represented by red and blue markers in figure 12 $b$ , respectively), with the latter distinguished by a significantly increased spatial lag. The intermediate-amplitude solution (ii) is an unstable saddle point (black marker), which serves as the critical dynamical threshold for state selection. In the reduced phase space, this threshold is defined by the stable manifold of the saddle – the separatrix – which partitions the phase plane into two distinct basins of attraction (shaded regions in figure 12 b). Initial conditions residing below this manifold asymptotically decay towards the quasi-linear equilibrium (i), while those exceeding this threshold undergo transient growth before settling into the solitary-wave-like state (iii). This classification offers an alternative explanation for the amplitude- and phase-lag-dependent switching behaviour observed in the unsteady Euler trajectories (see figure 5).

Figure 12.

Phase space and basins of attraction for the reduced dynamical system (5.14)–(5.15) at $H=1$ . Phase portraits in the $(\theta , a)$ plane are shown for three characteristic speed regimes: ( $a$ ) $c\lt c_2$ , featuring a unique stable focus (red circle) corresponding to the small-amplitude state I; ( $b$ ) $c_2\lt c\lt c_1$ , the bistable regime where two stable foci (red and blue circles) coexist with an intermediate saddle point (black marker); ( $c$ ) $c\gt c_1$ , where the system converges to a unique large-amplitude stable focus (blue circle), representing state II. In ( $b$ ), the thick black curve denotes the separatrix – defined by the stable manifold of the saddle – which partitions the phase space into two distinct basins of attraction (shaded regions). Grey streamlines represent the global vector field, while faint coloured curves depict sample trajectories illustrating the spiralling relaxation towards the equilibria.