2.1. Boussinesq approximation and upside-down symmetry

The Boussinesq approximation is commonly adopted in the study of stratified fluids when the density variation is small. For the physical system under consideration, under the Boussinesq approximation, both $\varDelta _1 = (\rho _2-\rho _1)/\rho _2$ and $\varDelta _2 = (\rho _3-\rho _2)/\rho _2$ are assumed small ( $\ll 1$ ) and of the same order of magnitude, and can be neglected except where they appear in terms multiplied by $g$ . Under this approximation, the linear long-wave speeds $c_0$ are roots of

(2.18) \begin{equation} \big [ (H_1+H_2 )c^2 - g_1^\prime H_1 H_2 \big ] \big [ ( H_2 + H_3) c^2 - g_2^\prime H_2 H_3 \big ] - H_1 H_3 \,c^4 =0, \end{equation}

and $\gamma$ is defined accordingly as

(2.19) \begin{equation} \gamma ^\pm = 1+ \frac {H_2}{H_1} - \frac {g_1^\prime H_2}{{c_0^\pm }^2} = \left [ 1+\frac {H_2}{H_3} - \frac {g_2^\prime H_2}{{c_0^\pm }^2} \right ]^{-1}\!. \end{equation}

Here, $g_1^\prime = g \varDelta _1$ and $g_2^\prime = g \varDelta _2$ denote reduced gravities (the notation $g^\prime$ will be used unambiguously in the case when $\varDelta _1=\varDelta _2$ ).

We note that, under the Boussinesq approximation, both the Euler and MMCC3 systems are invariant under the so-called ‘upside-down’ symmetry (see Guan et al. Reference Guan, Doak, Milewski and Vanden-Broeck2024 and Barros et al. Reference Barros, Choi and Milewski2020, respectively). For every steady solution pair $(\zeta _1,\zeta _2)$ of a given configuration with densities $\rho _0 (1-\varDelta _1)$ , $\rho _0$ , $\rho _0 (1+\varDelta _2)$ and undisturbed thickness of layers $H_1$ , $H_2$ , $H_3$ (from top to bottom), $(-\zeta _2, -\zeta _1)$ is a solution of the modified physical configuration with densities $\rho _0 (1-\varDelta _2)$ , $\rho _0$ , $\rho _0 (1+\varDelta _1)$ , where the undisturbed thickness of layers are $H_3$ , $H_2$ , $H_1$ (from top to bottom).

In the remaining part of this section, we discuss the characterisation of solitary waves provided by weakly nonlinear theory, with particular emphasis on mode-1, also known as the first baroclinic mode.

2.2. Korteweg-de Vries theory

If a weakly nonlinear theory is adopted, by considering unidirectional waves in the limit of $\alpha = a/H_i=O(H_i^2/\lambda ^2)=O(\epsilon ^2)$ ( $\epsilon \ll 1$ ), with $a$ and $\lambda$ typical values of amplitude and wavelength, a KdV equation can be derived for the upper interface $z=\zeta _1 (x,t)$ , i.e.

(2.20) \begin{equation} \zeta _{1,t} + c_0 \,\zeta _{1,x} + c_1 \,\zeta _1 \,\zeta _{1,x} + c_2 \,\zeta _{1,xxx}=0, \end{equation}

with coefficients $c_1$ and $c_2$ given by

(2.21) \begin{align} c_1 & = \frac {3}{2} c_0 \frac {\dfrac {\rho _3}{H_3^2} \gamma ^3 + \dfrac {\rho _2}{H_2^2} (1-\gamma )^3 -\dfrac {\rho _1}{H_1^2}}{\dfrac {\rho _3}{H_3} \gamma ^2 + \dfrac {\rho _2}{H_2} (1-\gamma )^2+\dfrac {\rho _1}{H_1}}, \\[-12pt]\nonumber \end{align}

(2.22) \begin{align} c_2 & = \frac {1}{6} c_0 \frac {\rho _3 H_3 \gamma ^2 + \rho _2 H_2 (1+\gamma +\gamma ^2)+\rho _1 H_1}{\dfrac {\rho _3}{H_3} \gamma ^2 + \dfrac {\rho _2}{H_2} (1-\gamma )^2+\dfrac {\rho _1}{H_1}}, \end{align}

where $\gamma$ is defined in (2.14). Within the KdV theory, the propagation modes are fully decoupled, and solutions for mode-1 (mode-2) are obtained by replacing $c_0$ by $c_0^+$ ( $c_0^-$ ) in (2.20) and (2.14). For each mode, ISWs exist for (2.20) as long as their wave speed $c$ is greater than $c_0$ , and their polarity depends strictly on the sign of the quadratic nonlinearity coefficient $c_1$ ( $c_2\gt 0$ regardless of the physical parameters considered): the ISW is of elevation (depression) for $\zeta _1$ if $c_1\gt 0$ ( ${\lt } 0$ ). On the other hand, according to the KdV theory, no ISWs exist in the case when the coefficient $c_1$ vanishes, i.e. when

(2.23) \begin{equation} \frac {\rho _3}{H_3^2} \gamma ^3 + \frac {\rho _2}{H_2^2} (1-\gamma )^3 - \frac {\rho _1}{H_1^2} =0, \end{equation}

often referred to as the criticality condition. Typically, and regardless of the propagation mode considered, for given densities $\rho _i$ ( $i=1,2,3$ ) and a fixed ratio $H_1/H_3$ , criticality is attained at one single point in the parameter space on the $(H_1/H, H_2/H)$ -plane. An exception arises when the stratification satisfies $\rho _2/\rho _1= \rho _3/\rho _2$ . In this unique case, if mode-1 is considered, criticality is maintained along the entire line $H_2/H = - ( 1+(\rho _3/\rho _1)^{1/2} ) H_1/H + 1$ (see Appendix B in Barros et al. Reference Barros, Choi and Milewski2020).

We remark that, alternatively, a KdV equation could have been written for the lower interface displacement $\zeta _2(x,t)$ :

(2.24) \begin{equation} \zeta _{2,t} + c_0 \,\zeta _{2,x} + \tilde {c}_1 \,\zeta _2 \,\zeta _{2,x} + c_2 \,\zeta _{2,xxx}=0. \end{equation}

Here $\tilde {c}_1 = c_1/\gamma$ and the equation can be easily recovered from (2.20) by setting $\zeta _1 = \zeta _2/\gamma$ .

Figure 2. (a) Shows the locus in the parameter space where the criticality condition (2.23) is satisfied, under the Boussinesq approximation with $\varDelta _1=\varDelta _2$ . The criticality curve (in blue), along which no ISWs exist according to the KdV theory, splits up into two components on the $(H_1/H, H_2/H)$ -plane, one of which is the line $H_2/H=-2 H_1/H + 1$ , corresponding to a fixed ratio $H_1/H_3=1$ . Also marked on this line is the point with coordinates $(9/26,4/13)$ at which a transition from non-existence to the existence of solitary waves to the mKdV equation (2.25) occurs. This is made clearer in (b) where the mKdV theory is adopted by relaxing the constraint $H_1=H_3$ , for which the locus (in thick black) for the coexistence of ISWs of opposite polarity is revealed.

To obtain the expression of the KdV coefficients under the Boussinesq approximation, as well as the corresponding criticality condition, it suffices to set $\rho _1 = \rho _2 = \rho _3 = \rho _0$ , along with $c_0$ defined implicitly by (2.18), and consider $\gamma$ defined by (2.19) instead of (2.14). By doing so, we recover the KdV equation presented by Yang et al. (Reference Yang, Fang, Tang and Ramp2010), based on the work by Benney (Reference Benney1966) (see also Benjamin Reference Benjamin1966; Grimshaw Reference Grimshaw1981).

It is well known that when the Boussinesq approximation is adopted, with equal density increments ( $\varDelta _1=\varDelta _2$ ) and the same undisturbed depths for the top and bottom layers ( $H_1=H_3$ ), then mode-1 solitary waves cannot exist according to KdV theory. The reason being that, for these parameters, we have $\gamma ^+=1$ , which trivially satisfies the criticality condition. On the other hand, if the modified KdV (mKdV) equation

(2.25) \begin{equation} \zeta _{1,t} + c_0^+ \,\zeta _{1,x} + \frac {1}{12} c_0^+ H_1(2H_1 + 3H_2) \,\zeta _{1,xxx} - \frac {1}{8 H_1^3} c_0^+ (8H_1 -9H_2) (\zeta _1^3)_x =0 \end{equation}

is used for the first baroclinic mode (with $c_0^+= \sqrt {g^\prime H_1}$ ), different predictions can be made (see Talipova et al. Reference Talipova, Pelinovsky, Lamb, Grimshaw and Holloway1999). When $H_1\gt (9/8) \,H_2$ (or equivalently $H_2/H\lt 4/13$ along the line $H_2/H = -2 H_1/H+1$ , as in figure 2 $a$ ), no mode-1 ISWs exist according to the mKdV theory, in agreement with the KdV theory. In contrast, when $H_1 \lt (9/8) \,H_2$ , two branches of solutions with opposite polarity can be predicted by the mKdV theory (see Grimshaw, Pelinovsky & Talipova Reference Grimshaw, Pelinovsky and Talipova1997). We note that for the same density stratification, i.e. $\delta =1$ , a mKdV equation can be written by relaxing the constraint $H_1=H_3$ . This equation is defined solely along the criticality curve. Hence, the locus for the coexistence of ISWs of opposite polarity, according to the mKdV equation, is a subset of the criticality curve. In figure 2 $(b)$ this subset is depicted in black. In the remaining subset (in blue) of the criticality curve, no ISWs can exist according to the theory.

We remark that, under the Boussinesq approximation, the KdV theory carries over the upside-down symmetry of the Euler equations. This results from (2.20) and (2.24) and the conditions

(2.26) \begin{align} c_0(H_3,H_2,H_1,1/\delta ) &= c_0(H_1,H_2,H_3,\delta ), \\[-12pt]\nonumber \end{align}

(2.27) \begin{align} c_1(H_3,H_2,H_1,1/\delta ) &= -\tilde {c}_1(H_1,H_2,H_3,\delta ), \\[-12pt]\nonumber \end{align}

(2.28) \begin{align} c_2(H_3,H_2,H_1,1/\delta ) &= c_2(H_1,H_2,H_3,\delta ), \end{align}

which are all satisfied throughout the entire parameter space.

2.3. Gardner theory

Following Choi & Camassa (Reference Choi and Camassa1999), we can extend the KdV equation (2.20) to various orders in an attempt to better describe larger amplitude waves. Namely, the Gardner equation:

(2.29) \begin{align} \zeta _{1,t} +c_0 \,\zeta _{1,x} + c_1 \zeta _1 \zeta _{1,x} + c_2 \,\zeta _{1,xxx} + c_3 \big(\zeta _1^3\big)_x = 0. \end{align}

This model is a truncation that mixes terms of different orders and, contrary to KdV, the terms may not balance as $\epsilon \rightarrow 0$ .

As expected, the coefficients $c_1$ , $c_2$ are the same as for the KdV equation, given by (2.21) and (2.22), respectively. The cubic nonlinearity coefficient $c_3$ is found (see Appendix A) as

(2.30) \begin{equation} c_3= \frac {T_1+T_2+T_3}{\dfrac {\rho _3}{H_3}\gamma ^2+\dfrac {\rho _2}{H_2}(1-\gamma )^2+\dfrac {\rho _1}{H_1}}, \end{equation}

with $T_1$ , $T_2$ , $T_3$ given by

(2.31) \begin{align} T_1 &=-\frac {1}{4} c_0 \Bigg [ 9 \frac {H_2}{\rho _2}\Bigg ( \frac {\rho _3}{H_3^2}\gamma ^2-\frac {\rho _2}{H_2^2}(1-\gamma )^2\Bigg )\Bigg (\frac {\rho _2}{H_2^2}(1-\gamma )^2-\frac {\rho _1}{H_1^2}\Bigg ) \nonumber \\ & \quad +\, 4\Bigg ( \frac {\rho _3}{H_3^3}\gamma ^4+\frac {\rho _2}{H_2^3}(1-\gamma )^4+\frac {\rho _1}{H_1^3}\Bigg )\Bigg ], \end{align}

(2.32) \begin{align} T_2 &=-\frac {1}{6} \frac {c_1}{\rho _2} \Bigg [ 10 \frac {\rho _2^2}{H_2^2}(1-\gamma )^3+\frac {3\rho _1}{H_1} \frac {\rho _3}{H_3} \left ( \frac {2 H_2}{H_1}-3\frac {H_2}{H_3}\gamma \right ) \,\gamma + 3\frac {\rho _2}{H_2}(1-\gamma )^2 \nonumber \\ & \quad \times\, \left ( 3 \frac {\rho _1}{H_1}-2\frac {\rho _3}{H_3}\gamma \right ) +\frac {\rho _1}{H_1}\frac {\rho _2}{H_1} (-1+6\gamma )+ \frac {\rho _2}{H_3} \frac {\rho _3}{H_3} \gamma ^2 \,(-9+4\gamma )\Bigg ], \end{align}

(2.33) \begin{align} T_3 &=-\frac {1}{6}\frac {c_1^2}{c_0} \Bigg [ -3\left ( \frac {\rho _3}{H_3}\gamma ^2+\frac {\rho _2}{H_2}(1-\gamma )^2+\frac {\rho _1}{H_1}\right ) \nonumber \\ & \quad +\, \frac {4 \gamma }{\rho _2 H_1 H_3} (\rho _1 \rho _2 H_3+\rho _1 \rho _3 H_2 +\rho _2 \rho _3 H_1)\Bigg ]. \end{align}

As before, to recover the expression of the Gardner coefficients under the Boussinesq approximation (see Kurkina et al. Reference Kurkina, Kurkin, Rouvinskaya and Soomere2015), it suffices to set $\rho _1 = \rho _2 = \rho _3 = \rho _0$ , along with $c_0$ defined implicitly by (2.18), and consider $\gamma$ defined by (2.19). Here, the cubic nonlinearity coefficient $c_3$ can change sign depending on the physical parameters considered. Given that the coefficient $c_2$ of the linear dispersive term $\zeta _{1,xxx}$ is always positive, it is well known that when $c_3\gt 0$ waves of either polarity exist with no maximum amplitude. Minimum wave amplitudes exist for waves of elevation (depression) for $c_1\lt 0$ ( ${\gt } 0$ ). Specifically, these have a minimum amplitude of $-2c_1/(3 c_3)$ . On the other hand, when $c_3\lt 0$ , a single branch of solutions exists with the same polarity as that predicted by the KdV equation, with the important difference that such solutions now broaden with increasing wave speeds (or amplitude), until a front is reached (see, e.g. Grimshaw, Pelinovsky & Talipova Reference Grimshaw, Pelinovsky and Talipova1999). A schematic illustration of all solution behaviours is given in figure 3 $(b)$ following Grimshaw et al. (Reference Grimshaw, Pelinovsky and Talipova1999).

Figure 3. Behaviour of the quadratic and cubic nonlinearity coefficients ( $c_1$ , $c_3$ , respectively) of the Gardner equation (2.29) for the upper interface displacement $\zeta _1(x,t)$ , when restricted to mode-1 waves. The stratification is given by (2.16) with $\varDelta _1=\varDelta _2=0.1$ . It is shown how the signs for the nonlinearity coefficients may vary across the whole parameter space. The colour scheme shown in panel (b) and used in panel (a) is the same as that set by Kurkina et al. (Reference Kurkina, Kurkin, Rouvinskaya and Soomere2015), and allows us to easily distinguish the various solution behaviours. Shaded regions in (b) represent a minimum amplitude along a solution branch.

Kurkina et al. (Reference Kurkina, Kurkin, Rouvinskaya and Soomere2015) have examined in great detail how the quadratic and cubic nonlinearity coefficients can vary with the different physical parameters considered in a three-layer fluid, under the Boussinesq approximation. In particular, it was highlighted through extensive numerical tests that the coefficient $c_3$ could never change sign for mode-2 waves, and remained negative. As a consequence, Gardner theory predicts that the coexistence of ISWs of opposite polarity can only occur for mode-1 waves. As in Kurkina et al. (Reference Kurkina, Kurkin, Rouvinskaya and Soomere2015), when graphically presenting the parameter space for different solutions, we may fix the density increments $\Delta _i$ and show the result of varying the undisturbed depths $H_i$ . The domain of physically permissible parameters is then a triangle with vertices $(0,0)$ , $(0,1)$ and $(1,0)$ in the $(H_1/H,H_2/H)$ -space. We show in figure 3 how the sign for the quadratic and cubic nonlinearity coefficients changes with $\varDelta _1=\varDelta _2=0.1$ . If we consider $\varDelta _1=\delta \varDelta _2$ , with $\delta \neq 1$ and $\varDelta _2=0.1$ , and draw similar pictures, we observe that, despite the somewhat large density contrasts considered, the results obtained are in good agreement with those obtained by Kurkina et al. (Reference Kurkina, Kurkin, Rouvinskaya and Soomere2015), reproduced here in figure 4. Note that when the Boussinesq approximation is considered, sweeping the parameter space becomes easier, since it is defined only by the three parameters $\delta$ , $H_1/H$ and $H_2/H$ . We can infer from the figures that given fixed densities, the coexistence of solutions of opposite polarity can be predicted within a limited range of values of $H_2/H$ , or none at all, depending on the ratio between the undisturbed thicknesses of the top and bottom layer (see figure 6).

Figure 4. Behaviour of the quadratic and cubic nonlinearity coefficients ( $c_1$ , $c_3$ , respectively) of the Gardner equation (2.29) for the upper interface displacement $\zeta _1(x,t)$ , when restricted to mode-1 waves under the Boussinesq approximation. Different values of $\delta$ are considered, and the same colour scheme as for figure 3 is adopted.

Figure 5. Behaviour of the quadratic and cubic nonlinearity coefficients ( $\tilde {c}_1$ , $\tilde {c}_3$ , respectively) of the Gardner equation (2.34) for the lower interface displacement $\zeta _2(x,t)$ , when restricted to mode-1 waves under the Boussinesq approximation. Different values of $\delta$ are considered, and the colour scheme remains consistent with that of figure 3; however, the coordinate axes have been redefined.

We remark that, alternatively, a Gardner equation could have been written for the lower interface displacement $\zeta _2(x,t)$ :

(2.34) \begin{equation} \zeta _{2,t} + c_0 \,\zeta _{2,x} + \tilde {c}_1 \,\zeta _2 \,\zeta _{2,x} + c_2 \,\zeta _{2,xxx} + \tilde {c}_3 \big(\zeta _2^3 \big)_x=0. \end{equation}

As shown in Kurkina et al. (Reference Kurkina, Kurkin, Rouvinskaya and Soomere2015), under the Boussinesq approximation, the following relationship holds:

(2.35) \begin{equation} c_3(H_3,H_2,H_1,1/\delta ) = \tilde {c}_3(H_1,H_2,H_3,\delta ). \end{equation}

This shows that the Gardner theory does preserve the upside-down symmetry of the Euler equations. Consequently, figure 4 $(c)$ can be transformed into figure 5 $(a)$ using the mapping $(x,y)\mapsto (1-x-y,y)$ , with the caveat of a reversed colour scheme (green replacing blue and orange replacing red), as we would be describing the behaviour of $-\zeta _2(x,t)$ . What is intriguing is that figures 5 $(a)$ and 4 $(a)$ exhibit clear discrepancies, suggesting that the predicted solution behaviour may depend on which form of the Gardner equation, (2.29) or (2.34), is used; cf. Kurkina et al. (Reference Kurkina, Kurkin, Rouvinskaya and Soomere2015). Similar considerations apply to other panels in figure 4.

Figure 6. Vanishing of the quadratic (in blue) and cubic (in purple) nonlinearity coefficients ( $c_1$ , $c_3$ , respectively) of the Gardner equation (2.29) for $\zeta _1(x,t)$ , when restricted to mode-1 waves under Boussinesq approximation. The parameter $\delta$ is set to 1, as in figure 4( $b$ ). The grey lines from the top vertex to the base of the triangle correspond to different ratios $H_1/H_3$ . Note that a fixed ratio $H_1/H_3=m$ determines the line with the equation $H_2/H = -(1+m^{-1} ) H_1/H +1$ in the parameter space. From left to right we have: $H_1/H_3=0.8, \, \approx 2.5155, \,20$ . For $H_1/H_3\gt 2.5155$ , according to the Gardner theory for $\zeta _1$ no ISWs of opposite polarity can coexist. This will be shown to be in disagreement with the predictions by the MMCC3 and Euler theories (e.g. see figure 9 $a$ ).

Before moving to the next section, we make additional remarks regarding the non-uniqueness of the Gardner equations: (i) other physical variables could be used in the expansion (such as the mean layer horizontal velocity) leading to yet a different solution behaviour in parts of parameter space; (ii) despite this non-uniqueness, physical considerations may, in some cases, be invoked to guide the selection of an appropriate Gardner model for specific parameter regimes; (iii) the curves corresponding to criticality coincide in the corresponding panels of figures 4 and 5, reflecting the fact that predictions for the solution behaviour of the KdV equation are uniquely defined and do not depend on the choice of dependent variable; (iv) at criticality, the subset of the parameter space in which the cubic nonlinearity coefficient is non-negative is likewise uniquely determined, irrespective of the particular Gardner equation adopted; (v) finally, and most importantly, we note that even when different Gardner formulations agree, this agreement does not, as will be demonstrated in the following sections, guarantee an accurate description of the strongly nonlinear or Euler results.

Figure 7. Possible solution spaces for solitary waves with speeds $c\gt c_0^+$ . Here, all solution branches are limited in amplitude. For any set of parameters, the bifurcation of mode-1 waves for the Euler and MMCC3 systems is described by one of the panels A–F. Furthermore, for $H_2$ sufficiently large, one finds mode-2 solitary waves, with a bifurcation structure described in Doak et al. (Reference Doak, Barros and Milewski2022), represented here by panel $M2$ . When applicable, shaded regions represent a minimum amplitude along a solution branch. Panels A and B show two branches of solitary waves of opposite polarity, both evolving into tabletop solitary waves. Panel B has a finite amplitude bifurcation for the depression branch, while the elevation branch bifurcates from zero amplitude. The opposite is true for panel A. Panels C and D show a single solitary wave branch of depression and elevation, respectively, developing into a tabletop solitary wave. Panels E and F are new bifurcation structures. There exist two branches, both of depression (panel E, $c_1\lt 0$ ) or elevation (panel F, $c_1\gt 0$ ). One branch (shown in black) bifurcates from zero amplitude and limits to a tabletop solitary wave. A second branch (shown in blue) exists with speeds strictly greater than those on the first branch. The branch starts at the same limiting tabletop solitary wave as the first branch, but the corresponding profiles include an additional bump on the broadened section. As one proceeds along the branch via continuation, the broadened section gradually diminishes until the solution takes the form of a conventional solitary wave. As the amplitude continues to increase, the wave broadens again and approaches a different tabletop solitary wave with a larger limiting amplitude. These new branches are described in greater detail in § 6.

3.1. Solution space and conjugate states

The MMCC3 model has recently been examined by Barros et al. (Reference Barros, Choi and Milewski2020) and Doak et al. (Reference Doak, Barros and Milewski2022). Notably, Barros et al. (Reference Barros, Choi and Milewski2020) have shown that ISWs are governed by a Hamiltonian system with two degrees of freedom, for which the potential $V=V(\zeta _1,\zeta _2)$ is a rational function (see (4.3)). Furthermore, Doak et al. (Reference Doak, Barros and Milewski2022) established that the solutions of

(3.1) \begin{equation} V=0, \quad \frac {\partial V}{\partial \zeta _1}=0, \quad \frac {\partial V}{\partial \zeta _2}=0 \end{equation}

correspond to uniform flows to the Euler equations, which are conjugate. This problem was first addressed by Lamb (Reference Lamb2000), who arrived to a set of three algebraic equations – (12a), (12b), (14) therein – equivalent to our system (3.1). To solve the system, Lamb (Reference Lamb2000) employed a geometric approach. However, a more insightful procedure involves plotting branches of conjugate state solutions of (3.1) across different parameter regimes (see, e.g. figure 5 in Doak et al. Reference Doak, Barros and Milewski2022), following the work of Dias & Il’ichev (Reference Dias and Il’ichev2001) and Barros (Reference Barros2016). This approach is particularly useful under the Boussinesq approximation, as it allows us to analyse how the number and speed of conjugate state solutions change with increasing values of $H_2/H$ for fixed $\delta$ and $H_1/H_3$ , as illustrated in figure 8. The black dashed lines represent the curves $c=c_0^\pm$ , corresponding to the linear long-wave speeds, and for small values of $H_2/H$ , a single branch of solutions with $c\gt c_0^+$ is present. For thicker intermediate layers, two or even three such branches can be perceived. The branches are colour coded based on the quadrant of the $(\zeta _1,\zeta _2)$ -plane to which the corresponding conjugate solutions belong. Orange and blue will typically correspond to solutions of mode-1, and red and green to solutions of mode-2. However, it can be observed from the figure that solutions sought in odd quadrants do exist with speeds less than $c_0^+$ , as well as solutions sought in even quadrants for speeds larger than $c_0^+$ (see the red branch for large values of $H_2/H$ ).

Figure 8. Solutions to the conjugate state (3.1), under the Boussinesq approximation. In each panel we plot the branches of conjugate flow speeds that exist for a fixed ratio $H_1/H_3$ and varying $H_2/H$ . The physical parameters here correspond to those in figure 6. Panels (a)–(c) show results for increasing values of $H_1/H_3$ , respectively $H_1/H_3=0.8, \, \approx 2.5155, \,20$ . The colour scheme in panel $(d)$ allows us to identify the quadrant of the $(\zeta _1,\zeta _2)$ -plane to which the corresponding conjugate state belongs to. The black dashed curves are the linear long-wave speeds.

A method has recently been proposed by Lamb (Reference Lamb2023) to determine where in the parameter space ISWs of both polarities exist in a continuous double-pycnocline stratification. The same method can equally be applied to our piecewise-constant stratification and consists of identifying the physical parameters for which a blue and orange branch coexist with $c\gt c_0^+$ . According to the author, ‘This method is somewhat tedious as it requires the determination of where the conjugate flow speed is equal to the linear long-wave propagation speed’. Our analysis will confirm that the onset of the desired feature is indeed predicted by instances where the conjugate flow speed matches the linear long-wave speed $c=c_0^+$ . However, it will also become evident that not all such instances lead to the expected feature (see figure 10). Furthermore, it will be shown that the conjugate state analysis extends beyond this point, and previously unexplored aspects of the conjugate flow speed curves, which for brevity, we will refer to as the conjugate curves, will enable the discovery of new types of solution branches; cf. figure 7.

In figure 8 $(c)$ the tangency between the conjugate curve and the curve $c=c_0^+$ is what sets the coexistence of the two branches of opposite polarity (orange and blue) with $c\gt c_0^+$ . In contrast, in figure 8 $(a)$ , it is the instance when the orange branch intersects the curve $c=c_0^+$ that sets such occurrence. Although not particularly visible in the figure, this intersection precedes a tangency between the orange branch and the curve $c=c_0^+$ (at $H_2/H\approx 0.398$ ). By preserving the same value $\delta =1$ and varying the ratio $H_1/H_3$ , similar figures can be obtained. By doing so, we can cover the entire parameter space as in figure 9 $(a)$ and highlight the regions where conjugate curves of opposite polarities exist for $c\gt c_0^+$ , simply by observing the intersections (blue dots) and tangency points (blue solid line) between the conjugate curve and the curve $c=c_0^+$ . This naturally raises the question of how this relates to the coexistence of ISWs with opposite polarities. Our extensive numerical tests confirm that when conjugate states of opposite polarities coexist with $c\gt c_0^+$ , wave fronts can be exhibited for each conjugate state (in both the strongly and fully nonlinear theories), ultimately limiting the amplitude of solitary-wave solution branches. This further supports the validity of the criterion proposed by Lamb (Reference Lamb2023). An additional justification for its validity will be presented in the following sections within the framework of strongly nonlinear theory.

Figure 9. Region (shaded) of coexistence of ISWs of opposite polarity for different values of $\delta$ , according to the strongly and fully nonlinear theories, under the Boussinesq approximation. Instances where the conjugate flow speed matches the linear long-wave speed $c=c_0^+$ are represented in blue (dots for intersections and a solid line for tangency points). It can be observed that the set of tangency points coincides with the locus corresponding to the criticality condition (2.23) for the KdV equation. The purple curves correspond to the vanishing of the cubic nonlinearity coefficient for the Gardner equation for $\zeta _1$ , and are shown to illustrate the stark difference between Gardner and MMCC3 predictions. According to the Gardner theory for $\zeta _1$ , coexistence of ISWs exist only in the region enclosed by this curve (see, e.g. figure 4 b).

In figure 9 we overlay the locus (in purple) where the cubic nonlinearity of the Gardner equation vanishes, allowing for a comparison with the weakly nonlinear theory. According to the Gardner theory, ISWs of opposite polarity coexist within each closed region enclosed by a purple curve. For $\delta =1$ , some discrepancies can be observed in panel $(a)$ . However, these differences become even more pronounced as $\delta$ varies (see panel b).

Returning to the point made earlier on unexplored aspects of the conjugate curves allowing us to unveil new types of ISWs, we first highlight the intricate behaviour of these curves, as depicted in figure 10. In the figure, with $\delta =0.5$ , $H_1/H_3=0.1$ , in addition to the intersections and tangencies between the conjugate curves and the curve $c=c_0^+$ , we also track the self-intersections (grey squares) of the conjugate curves, and where they first emerge (black square), provided $c\gt c_0^+$ . By varying $H_1/H_3$ while keeping $\delta$ constant, the entire parameter space can be explored. All these markers can then be compiled, as in figure 11 $(a)$ , to delimit regions where distinct ISWs properties arise. This classification is further illustrated in figure 11 $(b)$ , following the colour scheme outlined in figure 7 (see also figure 17 with $\delta =1$ to see how the complete solution space compares with figure 9 $a$ ). How these markers separate the bifurcation spaces, and why, is explored in the following sections.

Figure 10. Solutions to the conjugate state (3.1), under Boussinesq approximation, for $\delta =0.5$ and $H_1/H_3=0.1$ , using the same colour scheme as figure 8. The figure also includes vertical dotted curves, which deliminate the various solution space behaviours given by panels A–F and $M2$ in figure 7. Panel $(b)$ shows a blow-up of the region marked by a box in panel $(a)$ . In panel $(b)$ we highlight the appearance of two additional conjugate state branches (a maximum and a slower saddle) with $c\gt c_0^+$ with a black square, the self-intersection of two mode-1 conjugate states (both maxima) with $c\gt c_0^+$ with a grey square and mode-1 conjugate states passing through $c=c_0^+$ with a blue circle.

Figure 11. All details of our conjugate state analysis encapsulated in one figure. In panel $(a)$ we identify in the parameter region for $\delta =0.5$ all relevant instances for the conjugate curves. The main markers are defined in figures 9 and 10 $(b)$ . In addition, the red line denotes the set of all intersections between a red conjugate curve and the linear long-wave speed $c=c_0^+$ (see figure 10 $a$ for large values of $H_2/H$ ). This information can be used to delimit regions where different properties for mode-1 ISWs can be found, as in panel $(b)$ . The colour scheme used here is outlined in figure 7. The purple lines represent the condition $c_3=0$ , included for comparison with predictions from Gardner theory for $\zeta _1$ .

4.1. Geometric locus of the critical points

Consider fixed values of $\rho _i$ and $H_i$ ( $i=1,2,3$ ). Then, for a given wave speed $c$ , the set of critical points is found by solving $ {\partial V}/{\partial \zeta _1}= {\partial V}/{\partial \zeta _2}=0$ . The admissible solutions are those within the triangular region $\mathcal{R}$ :

(4.6) \begin{equation} {\mathcal{R}} = \{ (\zeta _1,\zeta _2)\in {\mathbb{R}}^2: \, h_i\gt 0, \quad i=1,2,3\}. \end{equation}

Given that the potential $V$ is a rational function, finding its critical points is equivalent to finding the intersection of two plane algebraic curves:

(4.7) \begin{align} {\mathcal{C}}_1 & \equiv \frac {1}{2} c^2 \big [ (\rho _2-\rho _1) h_1^2 h_2^2 + \rho _1 H_1^2 h_2^2 - \rho _2 H_2^2 h_1^2 \big ] -(\rho _2-\rho _1)g \,\zeta _1 h_1^2 h_2^2 =0, \\[-12pt]\nonumber \end{align}

(4.8) \begin{align} {\mathcal{C}}_2 & \equiv \frac {1}{2} c^2 \big [ (\rho _3-\rho _2) h_2^2 h_3^2 + \rho _2 H_2^2 h_3^2 - \rho _3 H_3^2 h_2^2 \big ] -(\rho _3-\rho _2)g \,\zeta _2 h_2^2 h_3^2 =0. \end{align}

The fact that the origin remains a critical point for all values of $c$ can be interpreted geometrically by noting that both ${\mathcal{C}}_1 = 0$ and ${\mathcal{C}}_2 = 0$ are homogeneous curves. Eliminating $c$ from (4.7) and (4.8) yields the following relationship between $\zeta _1$ and $\zeta _2$ :

(4.9) \begin{equation} P(\zeta _1,\zeta _2) \equiv (\rho _2-\rho _1) \zeta _1 \,h_1^2 \,G_3(h_2,h_3) - (\rho _3-\rho _2) \zeta _2 \,h_3^2 \,G_1 (h_1,h_2) =0. \end{equation}

Here

(4.10) \begin{align} G_1(h_1,h_2) & =(\rho _2-\rho _1)h_1^2 h_2^2 + \rho _1 H_1^2 h_2^2 - \rho _2 H_2^2 h_1^2, \\[-12pt]\nonumber \end{align}

(4.11) \begin{align} G_3(h_2,h_3) &= (\rho _3-\rho _2)h_2^2 h_3^2 + \rho _2 H_2^2 h_3^2 - \rho _3 H_3^2 h_2^2. \end{align}

Within the admissible region $\mathcal{R}$ , the two sets of (4.7), (4.8) and (4.8), (4.9) are equivalent. In particular, any critical point of $V$ lies on the geometric locus defined by the curve $P=0$ in (4.9). Conversely, it can be shown (see Lemma1 of Appendix B) that for any point $(\zeta _1, \zeta _2)$ on the curve $P=0$ there exists a real value of $c$ for which (4.8), (4.9) are satisfied, and hence $(\zeta _1, \zeta _2)$ is a critical point of $V$ .

Figure 12. Plots for the curve $P=0$ in (4.9). We set $(H_1,H_2,H_3)=(3,1,6)$ and adopt the stratification (2.16) with $\varDelta _1=\delta \varDelta _2$ and $\varDelta _2=0.1$ . $(a)$ $\delta =1$ , $(b)$ $\delta \approx 1.70824$ , $(c)$ $\delta =2$ . When $c\approx 0$ , there are seven critical points of $V$ : three local maxima (M), three saddles (S), one minimum (m).

The curve $P=0$ can be readily visualised in the $(\zeta _1,\zeta _2)$ -plane once all physical parameters are specified, revealing two main distinct configurations; cf. panels $(a)$ and $(c)$ of figure 12. In either case, the curve $P=0$ consists of three distinct branches: branch I, which connects vertices $(-H_2-H_3,-H_3)$ and $(H_1, H_1+H_2)$ ; branch II, which connects the point $ ( -H_2/ (1+\delta ), H_2 \delta /(1+\delta ) )$ to one of the two points $(0,-H_3)$ or $(H_1,0)$ ; branch III, which connects $(H_1,-H_3)$ to either $(0,-H_3)$ or $(H_1,0)$ . All endpoints of these branches correspond to critical points in the limit when $c\rightarrow 0^+$ . As $c$ increases, all points (except the origin) move inward. Three maxima (M) depart from the vertices of the triangle $\mathcal{R}$ ; three saddles (S) depart from each side of the triangle; the origin remains a minimum (m) until $c$ reaches the value $c_0^-$ . These points move according to two basic rules: (i) if two points ever meet at any location other than the origin, they cannot pass each other; (ii) as $c$ increases, any given point on a branch of $P=0$ must remain on that same branch. The first rule follows from a simple observation: from (4.7), or (4.8), it is evident that each point on $P=0$ , apart from the origin, corresponds to a unique value of $c$ . The second rule is more subtle and requires a technical justification, which is provided in Appendix D.

The origin is a somewhat special point. It is the only critical point that will not change its position as the wave speed varies. Moreover, it is not a regular point of $P=0$ . In fact, at this double point (or crunode), the tangents to the curve $P=0$ can be found as $\zeta _2=\gamma ^\pm \,\zeta _1$ as in (2.13) (see Lemma2 in Appendix B).

For certain exceptional values of the physical parameters, the curve $P=0$ may exhibit an additional singular point. Assuming the stratification given in (2.16) with $\varDelta _1=\delta \varDelta _2$ , and fixing the parameters $\varDelta _2$ , $H_i$ ( $i=1,2,3$ ), the specific value of $\delta$ that gives rise to an extra crunode can be determined by solving the system $P=0$ , $ {\partial P}/{\partial \zeta _1}=0$ and $ {\partial P}/{\partial \zeta _2}=0$ , simultaneously. In the particular case depicted in figure 12, this critical value is approximately $\delta \approx 1.70824$ , as shown in panel $(b)$ . If, instead, the Boussinesq approximation were applied, this transition between the two curve configurations would occur precisely at $\delta =2$ (corresponding to $=H_3/H_1$ ), as demonstrated in Appendix B.

4.2. Collision of critical points

Consider fixed values of $\rho _i$ and $H_i$ ( $i=1,2,3$ ). As the wave speed varies, the surface of the potential $V$ is smoothly deformed and the number and nature of critical points remains invariant, unless a collision takes place. Examining when collisions occur is thus a crucial step in the critical point analysis.

4.2.1. Trivial collisions

We begin by examining the collision of critical points at the origin, which we refer to as ‘trivial’ collisions. It is important to note that the origin is a regular point of both curves ${\mathcal{C}}_1 = 0$ and ${\mathcal{C}}_2 = 0$ . Specifically, the tangents to each curve at the origin are given, respectively, by

(4.12) \begin{equation} \begin{aligned} \big [ (\rho _1 H_2 + \rho _2 H_1)c^2 - g(\rho _2-\rho _1)H_1 H_2 \big ] \zeta _1 - \rho _2 H_1 c^2 \,\zeta _2 &=0,\\ \rho _2 H_3 c^2 \,\zeta _1 - \big [(\rho _2 H_3 + \rho _3 H_2)c^2 - g(\rho _3-\rho _2)H_2 H_3 \big ] \zeta _2 &=0. \end{aligned} \end{equation}

Moreover, it can be shown that $ {\partial {\mathcal{C}}_k}/{\partial \zeta _2} \neq 0$ for any point on the curve ${\mathcal{C}}_k$ ( $k=1,2$ ). We can therefore use the following result (Coolidge, Theorem 4, pp. 16–17): in the vicinity of an ordinary finite point where the tangent is not vertical, $\zeta _2$ may be expressed as a convergent power series in terms of $\zeta _1-\zeta _{1,0}$ , where $\zeta _{1,0}$ is the abscissa of the point in question. For our purposes, we set $\zeta _{1,0}=0$ and seek expansions:

(4.13) \begin{equation} \zeta _2 = p_k \,\zeta _1 + q_k \,\zeta _1^2 + r_k \,\zeta _1^3 + O(\zeta _1^4), \quad k=1,2. \end{equation}

The coefficients $p_k$ , $q_k$ , $r_k$ ( $k=1,2$ ) are all given in (C4) of Appendix C. A double collision at the origin holds when $p_1=p_2$ , which is precisely when (2.12) is satisfied. In other words, this type of collision arises only at the linear long-wave speeds $c=c_0^{\pm }$ , at which the tangent lines given in (4.12) become identical.

For a triple collision at the origin to occur, not only the first- but also the second-order approximation to the curves must coincide, i.e.

(4.14) \begin{equation} p_1=p_2, \quad q_1=q_2. \end{equation}

It is shown in Appendix C that (4.14) holds precisely at criticality, i.e. when (2.23) is satisfied.

Lastly, a quadruple collision at the origin occurs when

(4.15) \begin{equation} p_1=p_2, \quad q_1=q_2, \quad r_1 = r_2. \end{equation}

As shown in Appendix C, this situation arises precisely when, at criticality, $T_1=0$ , with $T_1$ defined by (2.31). In other words, a quadruple collision at the origin occurs exactly when both the quadratic and cubic coefficients, $c_1$ and $c_3$ , respectively, of Gardner’s equation vanish simultaneously. To understand this, note that the coefficient $c_3$ in (2.30) vanishes when $T_1+T_2+T_3=0$ . Since $T_2$ and $T_3$ are both proportional to $c_1$ , we deduce that, at criticality (i.e. when $c_1=0$ ), the cubic nonlinearity coefficient can only vanish if $T_1=0$ . Another important consequence is that quadruple collisions can only occur at $c=c_0^+$ .

In figure 3 we can identify in panel $(b)$ three instances in the parameter space where a quadruple collision occurs. In all other panels, only one instance of this type of collision is possible. To gain a better understanding of how many such collisions can occur in the parameter space, it is helpful to resort to the Boussinesq approximation. Under this reduction, a quadruple collision is equivalent to

(4.16) \begin{equation} \tilde {T}_1=0, \quad H_1^2 H_2^2 \,\gamma ^3 + H_1^2 H_3^2 \,(1-\gamma )^3 - H_2^2 H_3^2=0, \quad \gamma ^2+ \varLambda \gamma - \delta =0, \end{equation}

with $\varLambda =\delta (1+H_2/H_3)-1-H_2/H_1$ . The first condition in (4.16) results from viewing $T_1=0$ as the vanishing of a multivariate polynomial, similarly to (C21), but under Boussinesq approximation. The second condition is equivalent to the criticality condition. The third results from imposing that the collision occurs at the linear long-wave speeds (see (B.4) in Barros et al. Reference Barros, Choi and Milewski2020). Solving this system simultaneously through use of the Gröbner basis (see, e.g. Cox et al. Reference Cox, Little, O’shea and Sweedler1997; Cox, Little & O’Shea Reference Cox, Little and O’Shea2005), we arrive at a condition (more specifically, an algebraic curve) involving only $\delta$ and $H_1/H_3$ . This is very convenient, since with one single plot a full description can be provided, as in figure 13. For values of $\delta$ close to 1, there will be three instances where a quadruple collision occurs, whereas in general only one quadruple collision is expected. Specifically, for $\delta =1$ , quadruple collisions occur for the following values of $H_1/H_3$ : $\approx 0.397536, \, 1.0, \, \approx 2.515$ . In figure 6, quadruple collisions correspond to intersections between the blue and purple curves. Quadruple collisions help define regions within the parameter space where distinct ISW properties arise, as illustrated in figure 11 and further demonstrated in figures 17 and 18.

Figure 13. Instances in the parameter space where a quadruple collision occurs at the origin. Here, the Boussinesq approximation is adopted and, for each value of $\delta =\delta _0$ , the number of quadruple collisions may be read from the plot by counting how many times the horizontal line $\delta =\delta _0$ intersects the graph.

4.2.2. Non-trivial collisions

Whenever two critical points coalesce, they form a degenerate critical point, which can be identified by the vanishing of the determinant of the Hessian matrix in (4.4) at that specific location. Let $(\zeta _1,\zeta _2)\neq (0,0)$ be the coordinates of such a point. Then,

(4.17) \begin{equation} {\text{det}\,} {\mathcal{H}} (\zeta _1,\zeta _2)=0 \end{equation}

is a condition that also depends on the value of speed $c$ . However, given that $(\zeta _1,\zeta _2)$ is a critical point, the corresponding value of $c$ can be easily determined from (4.8) and substituted back into (4.17). Eliminating denominators and simplifying the expression leads to a family of plane algebraic curves given by $F(\zeta _1,\zeta _2)=0$ , which depends on the parameters $\rho _i$ , $H_i$ ( $i=1,2,3$ ). Once these parameters are fixed, the coordinates $(\zeta _1,\zeta _2)$ corresponding to non-trivial collisions can be found as the intersection points of the two plane algebraic curves $P=0$ and $F=0$ . Figure 14 shows some demonstrative examples of the locations of non-trivial collisions. This figure will be revisited in § 5 to explain the bifurcation of panels $A$ – $D$ in figure 7. Note that while the explicit expression for $F$ is unfortunately too cumbersome to be presented here, the results in figure 14 can be reproduced using any symbolic programme, such as Mathematica. In order to enable a more thorough investigation of collisions, we may employ the same methodology used for the conjugate states (see figure 8) to generate figure 15. This approach reveals not only the quadrants of the $(\zeta _1,\zeta _2)$ -plane in which collisions occur, but also the corresponding wave speeds at which these events take place. As illustrated, for fixed values of the density and undisturbed thickness of each layer, there exist either three or five non-trivial collisions.

Figure 14. Plots of the curve $P=0$ in (4.9) and points where non-trivial collisions occur. The stratification is given by (2.16), with $\varDelta _1= \varDelta _2=0.1$ , and the undisturbed thicknesses of each layer parameters are set as follows: $(a)$ $(H_1,H_2,H_3)=(2,6,1)$ , $(b)$ $(H_1,H_2,H_3)=(1,7,2)$ , $(c)$ $(H_1,H_2,H_3)=(3,1,6)$ , $(d)$ $(H_1,H_2,H_3)=(4,1,2)$ . The points in black (red) correspond to collisions leading to a decrease (increase) of the number of critical points. In each panel we overlay the zero contour of the potential $V$ , obtained for $(c/c_0^+)^2=1.01$ . The mode-1 ISWs for the physical parameters in panels $(a)$ – $(d)$ resemble, respectively, those illustrated in panels $A$ – $D$ of figure 7.

Figure 15. Collision curves on the $(H_2/H,c/\sqrt {gH})$ -plane for the stratification (2.16), with $\varDelta _1= \varDelta _2=0.1$ . $(a)$ $H_1/H_3=0.5$ , $(b)$ $H_1/H_3=2$ . The black dashed curves are the linear long-wave speeds, at which trivial collisions occur. Non-trivial collisions are represented by dotted (solid) lines when leading to a decrease (increase) of the number of critical points. We adopt the same colour scheme as in figure 8 to reveal the quadrant of the $(\zeta _1,\zeta _2)$ -plane where the collisions occur.

In general, a non-trivial collision leads to a decrease in the number of critical points. This is an immediate consequence of the first rule of how critical points evolve: since two points cannot move past one another, they must annihilate upon collision. This also poses a constraint on the nature of the critical points that collide. As a result of a theorem by Kronecker (see Theorem 2.4.6 in Emelyanov et al. Reference Emelyanov, Korovin, Bobylev and Bulatov2007), the total sum of topological indices of critical points is 1 and must remain invariant as the value of $c$ is increased (see Barros Reference Barros2016 for more details). Thus, no local maximum (or minimum) of the potential $V$ can collide with another point other than a saddle. In addition, after all collisions have taken place, only a single critical point remains: the origin, as a local maximum of $V$ . Let $c=c_{\!f}$ ( $\gt c_0^+$ ) denote the speed at which the final collision takes place. Beyond this speed, no solitary-wave solutions can exist. This critical speed thus serves as an upper bound for the propagation speed of mode-1 ISWs, as illustrated by the blue dotted curve in figure 15 $(a)$ , and it inherently limits their amplitude.

We now turn to the exceptional case in which, unexpectedly, two critical points emerge spontaneously. These occurrences are indicated by red markers in figure 14 and by solid lines in figure 15. Panels $(b)$ and $(c)$ of figure 14 are related to figure 15 $(a)$ , since in both cases $H_1/H_3=0.5$ . Similarly, panels $(a)$ and $(d)$ of figure 14 are related to figure 15 $(b)$ , since in both cases $H_1/H_3=2$ . From figure 15, emergent collisions are to be expected for large values of $H_2/H$ , which explains the contrast between the bottom and top panels of figure 14. While `typical’ collisions – where critical points annihilate – occur along all branches of the curve $P=0$ , our extensive numerical investigations (see figures 15 and 16) indicate that such emergent collisions occur exclusively along branch I. From this, we conclude that these events correspond to the simultaneous emergence of a saddle point and a local maximum. In the two cases considered in figure 15, such occurrences are observed only for speeds below $c_0^+$ . However, as illustrated in figure 16, emergent collisions can also arise at speeds exceeding $c_0^+$ .

Figure 16. In panel $(a)$ , similar to figure 15, collision curves are displayed on the $(H_2/H,c/\sqrt {gH})$ -plane for the parameters $\varDelta _1=0.2$ , $\varDelta _2=0.1$ and $H_1/H_3=10$ . The two other panels illustrate distinct situations in panel $(a)$ , obtained with different values of $H_2/H$ : $(b)$ $H_2/H=0.28$ , $(c)$ $H_2/H=0.6$ . In both plots we set $H_1=1$ and overlay the zero contour of the potential $V$ , obtained for $(c/c_0^+)^2=1.01$ .

5.1. Criteria for the coexistence of ISWs of opposite polarity (panels A and B)

Based on the discussion above, we conclude that an emergent collision is a necessary condition for the coexistence of mode-1 waves with opposite polarities. This is exemplified by the physical parameters corresponding to the top panels of figure 14. However, as illustrated in figures 16 $(b)$ and 16 $(c)$ , the presence of an emergent collision alone is not sufficient: solutions of opposite polarity fail to materialise in both cases. To guarantee the coexistence of solutions of opposite polarity, we require that, at $c=c_0^+$ , two local maxima of $V$ exist on either side of branch I, with both satisfying $V\gt 0$ . In light of how critical points evolve, this implies that (i) an emergent collision (yielding a saddle and a maximum) must occur at $c\lt c_0^+$ , (ii) the resulting maximum must reach the origin as $c$ approaches $c_0^+$ . Furthermore, the polarity of the solution branch bifurcating from zero amplitude at the linear long-wave speed $c=c_0^+$ is determined by the quadrant of the $(\zeta _1,\zeta _2)$ -plane where the emergent collision takes place. Specifically, in figure 14 $(a)$ the emergent collision occurs in the first quadrant, resulting in a depression branch bifurcating from zero amplitude, while the elevation branch bifurcates from finite amplitude. A similar situation, but with the roles of the elevation and depression branch reversed, can be observed in figure 14 $(b)$ , where the emergent collision occurs in the third quadrant. In the regimes considered, the solutions resemble those illustrated in panels $A$ and $B$ of figure 7, respectively.

Identifying the parameter region for the coexistence of solutions with opposite polarity is quite challenging when relying solely on the analysis of collisions. Based on figure 15, one might be inclined to think that the desired feature holds for large values of $H_2/H$ , since collision curves in orange and blue exist for $c\gt c_0^+$ (following an emergent collision at $c\lt c_0^+$ ). This is a clear indication that, at $c=c_0^+$ , $V$ has two maxima on either side of branch I, which eventually collide with saddle points as $c$ increases. However, this approach does not a priori provide a guarantee of the desired outcome, as we cannot confirm if, originally at $c=c_0^+$ , $V\gt 0$ at the two maxima. This can be overcome by observing that, for that to happen, the two loops on either side of branch I should both shrink to a point, as $c$ increases from the value $c=c_0^+$ . When that happens (for each loop), a maximum becomes an equilibrium at the same energy level as the origin, i.e. $V= {\partial V}/{\partial \zeta _1}= {\partial V}/{\partial \zeta _2}=0$ , and thus, a conjugate state is reached and (3.1) is recovered. Therefore, if conjugate states of opposite polarity coexist for $c \gt c_0^+$ , the presence of ISWs with opposite polarity is assured. This finding supports the validity of the criterion proposed by Lamb (Reference Lamb2023) (see § 3.1), at least within the scope of the strongly nonlinear theory.

5.2. Criteria to distinguish panels C–F

The above highlights the benefits of using a conjugate state analysis to identify specific behaviours of solutions. Accordingly, the criteria presented hereafter are based on the properties of the conjugate states for this physical system.

Assume, for a given set of parameters, that conjugate states of opposite polarity, both with $c\gt c_0^+$ , do not exist. In this case, extensive numerical tests reveal that, for all $c\gt c_0^+$ , all critical points of the potential $V$ along branch I lie in the first (third) quadrant if $c_1\gt 0$ ( $c_1\lt 0$ ). Hence, the KdV theory correctly predicts the polarity of the solitary waves in panels $C$ – $F$ . We note that it is also possible, for the chosen parameters, that a conjugate state exists on branch II. As a result, in addition to the waves discussed above, there are other solutions exhibiting mode-2 characteristics, which will be addressed separately in § 5.3.

In terms of conjugate states for the system, a few cases must be considered. The simplest case is that there is only one mode-1 conjugate state with $c\gt c_0^+$ . In such cases, there exists a single solution branch of mode-1 solitary waves (panel $C$ or $D$ ), starting from zero amplitude and limiting to a tabletop solitary wave.

The case when there are two or three mode-1 conjugate states with $c\gt c_0^+$ in the first or third quadrant requires more delicate consideration. Without loss of generality, consider the case where two or three mode-1 conjugate states with $c\gt c_0^+$ exist in the first quadrant. For this to occur, we must have $c_1\gt 0$ , and the emergent collision of critical points must have occurred in the first quadrant along branch I, such as observed in figure 16(b). Note that this emergent collision occurs in between the origin and the maximum departing from $(H_1,H_1+H_2)$ in the limit when $c\rightarrow 0^+$ , since each point along the curves $P=0$ corresponds to a unique value of $c$ . The criterion for distinguishing whether we have one branch of solutions (panel $D$ ) or two branches (panel $F$ ) is as follows. As long as the speed at which the emergent maximum becomes a conjugate state is less than the speed at which the other maximum does, we have panel $F$ . Otherwise, we have panel $D$ . Analogous consideration applies to the case when two or three mode-1 conjugate states with $c\gt c_0^+$ exist in the third quadrant. This is shown, for example, in figure 10. At the black square, where the emergent collision occurs at $V=0$ with $c\gt c_0^+$ , the speed of the emergent maximum conjugate state (given by the black square) is less than the other maximum (the blue curve above the black square). However, as $H_2$ increases, the two front curves cross at the grey square. Hence, past the grey square, the speed of the emergent maximum as a conjugate state is faster than that of the other maximum. We see here that the black square and grey square signal changes in behaviour, while the crossing of the front curves through $c=c_0^+$ (blue circle) in this case do not. Greater detail on the solution space for panels $E$ and $F$ is given in § 6. Next, we use the above criteria to present the Boussinesq solution space for all parameters.

Figure 17. Partition of the (Boussinesq) parameter space into regions associated with distinct solution types. The dotted curves in each panel are the instances in the parameter space at which the front curve crosses the curve $c=c_0^+$ . The blue curve corresponds to criticality, at which the front curve and the curve $c=c_0^+$ are tangent. The other markers in the panels (when applicable) are the black and grey squares discussed in figure 11 $(b)$ . The red curve represents the instance when the front curve bearing characteristics of mode-2 crosses the curve $c=c_0^+$ . In the region above this red line, not only coexistence of ISWs of opposite polarity exist, but also large amplitude solutions with mode-2 characteristics. The panels show results for different values of $\delta$ . The colour scheme used here is outlined in figure 7. The purple lines represent the condition $c_3=0$ , included for comparison with predictions from Gardner theory for $\zeta _1$ .

Figure 18. Same as in figure 17, but for $\delta \gt 1$ .

5.3. Parameter sweep of the Boussinesq solution space

The above provides a mechanism for determining the solution space from the behaviour of the conjugate states. Figures 17 and 18 demonstrate how the solution space evolves as one varies $\delta$ within the Boussinesq regime. The colour scheme corresponds to panels A–F in figure 7.

Plots such as figure 8 show conjugate states as a function of $H_2$ for given $\delta$ and varying $H_1/H_3$ , where each value of $H_1/H_3$ represents a ray in the parameter space, as shown in figure 6. To create figures 17 and 18, a representative set of rays spanning the parameter space is selected. From each ray, key markers are extracted as illustrated in figure 10, including blue circles (mode-1 conjugate state with $c=c_0^+$ ), black squares (emergent collision coinciding with being a conjugate state), grey squares (two conjugate states in the same quadrant having the same speed) and the red curve (mode-2 conjugate state with $c=c_0^+$ ). All of these points are recovered numerically. The plots also display blue and purple lines, which are obtained analytically. As can be seen from figures 8 and 10, each point of tangency between conjugate curves and $c=c_0^+$ is marked by a colour change (from blue to orange or vice versa). Within the MMCC3 theory, this means that tangencies correspond to a triple collision at the origin, which holds at criticality – specifically when $c_1=0$ in the KdV equation (see Appendix C) – and is depicted as a blue line. The purple lines represent the condition $c_3=0$ , included for comparison with predictions from Gardner theory.

For small values of $H_2$ , the solution space always corresponds to either panel $C$ or $D$ , depending on the sign of $c_1$ . As $H_2$ increases, various events – such as the emergence of critical points, the crossing of front curves, their intersection with the curve $c=c_0^+$ and criticality – can significantly affect the structure of the solution space. Figures 17 and 18 capture all the key topological transitions observed in the solution space. More precisely, between two consecutive panels (i.e. successive values of $\delta$ ) the structure of the solution space remains qualitatively similar to that corresponding to the lower value of $\delta$ . It is important to note that, although not indicated in the colour scheme, for any point in parameter space lying above the red curve, along with the solutions described by panel $A$ or $B$ we also have solutions with characteristics of mode-2 corresponding to panel $M2$ in figure 7 (see Doak et al. Reference Doak, Barros and Milewski2022 for further details).

Finally, it is worth noting that any coloured panel corresponding to $\delta \gt 1$ can be transformed into one with $1/\delta$ via the mapping $(x,y)\mapsto (1-x-y,y)$ , with the caveat that the colour scheme is reversed. This reflects the upside-down symmetry inherent in both the Euler and MMCC3 systems, and also highlights the fact that any behaviour exhibited by one interface is mirrored by the other.