2.1. Governing equations

Figure 1.

Schematic of the problem in a cylindrical coordinate system.

In this investigation, we will explore an axisymmetric column system consisting of two superimposed, incompressible and inviscid fluids. The axisymmetric cylindrical coordinate system $(x,r)$ is introduced, such that the $x$ -axis aligns with the central axis of the rod, and $r$ represents the radial coordinate (refer to the sketch of the system in figure 1). The inner layer contains a column of ferrofluid characterised by a constant density $\rho ^i$ and magnetic susceptibility $\chi ^i$ , which coats a copper rod of radius $d$ . This ferrofluid is surrounded by an outer layer of a non-magnetisable fluid (with $\chi ^\circ=0$ ) that has a density  $\rho ^\circ$ less than that of the inner fluid. In this context, the superscripts ‘ $i$ ’ and ‘ $o$ ’ refer to the properties of the inner and outer fluid layers, respectively. The two fluids are assumed to be immiscible and are separated by a sharp interface denoted as $r=S(x,t)$ , which fluctuates around a mean radius $R$ . The entire system is contained within a fixed cylindrical container with a radius $D$ . An azimuthal magnetic field, $\boldsymbol{B}=({\mu _0\boldsymbol{I}}/{2\pi r})\boldsymbol{e}_{\theta }$ , is induced by a metallic rod that carries a current $\boldsymbol{I}$ , where $\mu _0$ represents the magnetic permeability of vacuum and $\boldsymbol{e}_{\theta }$ is the unit vector pointing in the clockwise azimuthal direction. The flow of each fluid is considered irrotational; thus, we introduce velocity potentials $\phi ^{i,o}$ , whose gradients represent the velocity fields. All physical variables can be non-dimensionalised as

(2.1) \begin{equation} \tilde {r}=\frac {r}{R},\quad \tilde {x}=\dfrac {x}{R},\quad \tilde {S}=\dfrac {S}{R},\quad \tilde {t}=t\sqrt {\frac {\sigma }{\rho ^iR^3}},\quad \tilde {\phi }^{i,o}=\phi ^{i,o}\sqrt {\frac {\rho ^i}{\sigma R}}, \end{equation}

where $\sigma$ is the surface tension coefficient. Using the provided scales and omitting tildes for ease of notation, we can express the interface as a perturbation around $r=1$ , defined as $S=1+\eta$ , where $\eta$ is the displacement from equilibrium. Consequently, we derive a set of dimensionless governing equations. Within the respective fluid domains, $\phi ^{i,o}$ satisfy the axisymmetric Laplace equation:

(2.2) \begin{equation} \begin{aligned} \phi ^\circ_{rr}+\frac {1}{r}\phi ^\circ_r+\phi ^\circ_{xx}=0, \quad &\text{for} \,\,1+\eta (x,t)\lt r\lt D,\\ \phi ^i_{rr}+\frac {1}{r}\phi ^i_r+\phi ^i_{xx}=0, \quad &\text{for} \,\,d\lt r\lt 1+\eta (x,t). \end{aligned} \end{equation}

At the interface $r=1+\eta (x,t)$ , the nonlinear kinematic boundary condition reads

(2.3) \begin{equation} \frac {\partial \eta }{\partial t}=\frac {\partial \phi ^i}{\partial r}-\frac {\partial \eta }{\partial x}\frac {\partial \phi ^i}{\partial x}=\frac {\partial \phi ^\circ}{\partial r}-\frac {\partial \eta }{\partial x}\frac {\partial \phi ^\circ}{\partial x}. \end{equation}

When neglecting the effect of gravity, the magnetic field and surface tension jointly provide restoring forces for fluctuations at the interface. Under the assumption of a linear magnetisation law (see Rosensweig (Reference Rosensweig2013, § 5.1)), the dynamic boundary condition can be expressed as

(2.4) \begin{equation} \phi _t^i-\rho \phi _t^\circ +\frac 12\big [ |\boldsymbol{\nabla }\phi ^i |^2-\rho |\boldsymbol{\nabla }\phi ^\circ |^2\big ]+\kappa -1-\frac {B}{2}\left [\frac {1}{(1+\eta )^2}-1\right ]=0, \end{equation}

where $\boldsymbol{\nabla }= (\partial _x,\partial _r )$ represents the gradient operator, $\rho =\rho ^\circ/\rho ^i$ is the density ratio, $\boldsymbol{B}= {\mu _0\chi {\boldsymbol{I}}^2}/{4\pi ^2\sigma R}$ stands for the magnetic Bond number and $\kappa$ denotes the mean curvature in the axisymmetric cylindrical coordinates, namely

(2.5) \begin{equation} \kappa =\dfrac {1}{(1+\eta )(1+\eta _x^2)^{1/2}}-\dfrac {\eta _{xx}}{(1+\eta _x^2)^{3/2}}. \end{equation}

At the surface of the rod and the outer container wall, the impermeability boundary conditions must be satisfied, specifically,

(2.6) \begin{equation} \begin{aligned} \phi _r^i=0,\qquad &\text{at}\,\,r=d, \nonumber \\ \phi _r^\circ=0,\qquad &\text{at}\,\,r=D. \end{aligned} \end{equation}

Finally, (2.2)–(2.6) form a closed mathematical system for the problem.

2.2. Hamiltonian structure

In a one-layer configuration, where the outer fluid is considered hydrodynamically passive, our previous work (Xu & Wang Reference Xu and Wang2025) demonstrates that the axisymmetric ferrofluid jet problem is a Hamiltonian system. In this study, we assert that the same holds for the two-layer compound ferrofluid jet problem. The total energy of the two-layer system can be expressed as

(2.7) \begin{align} \mathcal{H}\big [\phi ^{i,o},\eta \big ]=&\underbrace {\frac 12\int \int _d^{1+\eta }r|\boldsymbol{\nabla }\phi ^i|^2\,\textrm{d}{\kern0.5pt}r\,\textrm{d}{\kern0.5pt}x}_{E_k^i}+\underbrace {\frac {\rho }{2}\int \int _{1+\eta }^Dr|\boldsymbol{\nabla }\phi ^\circ|^2\,\textrm{d}{\kern0.5pt}r\,\textrm{d}{\kern0.5pt}x}_{E_k^\circ} \nonumber \\ &\underbrace {+\frac {B}{2}\int \left [\frac 12(2\eta +\eta ^2)-\ln (1+\eta )\right ] \textrm{d}{\kern0.5pt}x+\int \left [(1+\eta )\sqrt {1+\eta _x^2}-\frac 12(2+2\eta +\eta ^2)\right ] \textrm{d}{\kern0.5pt}x}_{E_{\!P}}, \end{align}

where $E_{\!K}$ ( $=E_k^i+E_k^\circ$ ) denotes the kinetic energy and $E_{\!P}$ represents the potential energy derived from the magnetic force and surface tension. By introducing the interfacial velocity potentials $\xi ^{i,o}=\phi ^{i,o} (x,1+\eta ,t )$ , we can demonstrate that (2.7) is a Hamiltonian, with the interfacial quantities $\xi =\xi ^i-\rho \xi ^\circ$ and $\zeta ={(1+\eta )^2}/{2}$ forming a pair of canonical variables. To achieve this, it is sufficient to take the variational derivatives of $E_{\!K}$ with respect to $\xi$ and $\zeta$ . First, we perturb $E_{\!K}$ with $\delta\xi ^i$ and $\delta\xi ^\circ$ while keeping $\eta$ fixed. Upon noticing that $\delta\xi ^{i,o}=\delta\phi ^{i,o}\big |_{r=1+\eta }$ , we calculate the variation of $E_{\!K}$ with respect to $\xi$ as

(2.8) \begin{align} \delta E_{\!K}&=\frac {\delta}{2}\left (\int\! \int _d^{1+\eta }\boldsymbol{\nabla }\phi ^i\boldsymbol{\cdot }r\boldsymbol{\nabla }\phi ^i\,\textrm{d}{\kern0.5pt}r\,\textrm{d}{\kern0.5pt}x+ \rho \int\! \int _{1+\eta }^D\boldsymbol{\nabla }\phi ^\circ\boldsymbol{\cdot }r\boldsymbol{\nabla }\phi ^\circ\,\textrm{d}{\kern0.5pt}r\,\textrm{d}{\kern0.5pt}x\right ) \nonumber \\ &=\int r\!\left (\delta\phi ^i\frac {\partial \phi ^i}{\partial \boldsymbol{n}}-\rho \delta\phi ^\circ\frac {\partial \phi ^\circ}{\partial \boldsymbol{n}}\right )\sqrt {1+\eta _x^2}\,\bigg |_{r=1+\eta }\textrm{d}{\kern0.5pt}x\,=\int \delta\xi (1+\eta )\frac {\partial \phi ^i}{\partial \boldsymbol{n}}\sqrt {1+\eta _x^2}\,\textrm{d}{\kern0.5pt}x, \end{align}

where $\boldsymbol{n}=(-\eta _x,1)^\top /\sqrt {1+\eta _x^2}$ represents a unit normal vector at the interface pointing outward, and the last equality arises from the divergence theorem in cylindrical coordinates. It follows immediately that

(2.9) \begin{equation} \frac {\delta \mathcal{H}}{\delta\xi }=\frac {\delta E_{\!K}}{\delta\xi }=(1+\eta )\left (\frac {\partial \phi ^i}{\partial r}-\frac {\partial \eta }{\partial x}\frac {\partial \phi ^i}{\partial x}\right )\overset {\text{or}}{=} (1+\eta )\left (\frac {\partial \phi ^\circ}{\partial r}-\frac {\partial \eta }{\partial x}\frac {\partial \phi ^\circ}{\partial x}\right )=\zeta _t, \end{equation}

which is equivalent to the kinematic boundary condition. Next, we fix $\xi$ and vary $\zeta$ ; thus, it is readily observed that

(2.10) \begin{align} 0&=\delta\xi = (\delta\phi ^i-\rho \delta\phi ^\circ )|_{r=1+\eta }+\left (\frac {\partial \phi ^i}{\partial r}-\rho \frac {\partial \phi ^\circ}{\partial r}\right )\delta\eta \nonumber\\ &\Longrightarrow (\delta\phi ^i-\rho \delta\phi ^\circ)|_{r=1+\eta }=-\left (\frac {\partial \phi ^i}{\partial r}-\rho \frac {\partial \phi ^\circ}{\partial r}\right )\delta\eta. \end{align}

Upon noting that $\delta\zeta =(1+\eta )\delta\eta$ , we have

(2.11) \begin{align} \delta E_{\!K}&=\frac 12\int (r|\boldsymbol{\nabla }\phi ^i|^2)|_{r=1+\eta }\delta\eta \,\textrm{d}{\kern0.5pt}x+\int\! \left (\delta\phi ^ir\frac {\partial \phi ^i}{\partial \boldsymbol{n}}\right )\bigg |_{r=1+\eta }\sqrt {1+\eta _x^2}\,\textrm{d}{\kern0.5pt}x \nonumber\\ &\quad-\frac \rho 2\int (r|\boldsymbol{\nabla }\phi ^\circ|^2)|_{r=1+\eta }\delta\eta \,\textrm{d}{\kern0.5pt}x - \rho \int\! \left (\delta\phi ^\circ \frac {\partial \phi ^\circ}{\partial \boldsymbol{n}}\right )\bigg |_{r=1+\eta }\sqrt {1+\eta _x^2}\,\textrm{d}{\kern0.5pt}x \nonumber\\ &=\int\! \left [\frac{1}{2} (|\boldsymbol{\nabla }\phi ^i|^2-\rho |\boldsymbol{\nabla }\phi ^\circ|^2 )-\left (\frac {\partial \phi ^i}{\partial r}-\rho \frac {\partial \phi ^\circ}{\partial r}\right )\eta _t\right ]\delta\zeta \,\textrm{d}{\kern0.5pt}x, \end{align}

indicating

(2.12) \begin{equation} \dfrac {\delta E_{\!K}}{\delta\zeta }=\frac 12 (|\boldsymbol{\nabla }\phi ^i|^2-\rho |\boldsymbol{\nabla }\phi ^\circ|^2 )-\left (\frac {\partial \phi ^i}{\partial r}-\rho \frac {\partial \phi ^\circ}{\partial r}\right )\eta _t. \end{equation}

By combining (2.12) with the variational derivative of $E_{\!P}$ with respect to $\zeta$ , a straightforward calculation leads to

(2.13) \begin{equation} \xi _t=-\dfrac {\delta\mathcal{H}}{\delta \zeta }, \end{equation}

which is equivalent to the dynamic boundary condition at the interface. (2.9) and (2.13) affirm that $\zeta$ and $\xi$ form a pair of canonical variables that comply with Hamilton’s equations.

To further simplify the Hamiltonian, we introduce the axisymmetric DNOs, denoted as $G^i$ and $G^\circ$ for the inner and outer layers, respectively. These operators are defined in a self-adjoint manner (see Xu & Wang Reference Xu and Wang2025, § 2.4):

(2.14) \begin{align} G^i(\eta ,d)\xi ^i=(1+\eta )\big (\phi ^i_{r}-\phi ^i_{x}\eta _{x}\big ), \end{align}

(2.15) \begin{align} G^\circ(\eta ,D)\xi ^\circ=(1+\eta )\big (\phi ^\circ_{r}-\phi ^\circ_{x}\eta _{x}\big ). \end{align}

The dependence of the operators on $\eta$ and $d$ ( $D$ ) is suppressed hereafter for ease of notation. These pseudodifferential operators can be expanded as recursively convergent series, $G^{i,o}=\sum {G_{\!j}}^{\!i,o}$ ; for instance, the first three terms for $G^i$ are as follows:

(2.16) \begin{gather} G^i_{0}=\frac {|P|\left [K_1(|P|d)I_1(|P|)-I_1(|P|d)K_1(|P|)\right ]}{K_1(|P|d)I_0(|P|)+I_1(|P|d)K_0(|P|)}, \\[-30pt] \nonumber \end{gather}

(2.17) \begin{gather} G^i_{1}=-\partial _x\eta \partial _x-G^i_0\eta G^i_0, \\[-30pt] \nonumber \end{gather}

(2.18) \begin{gather} G^i_{2}=\frac 12G^i_{0}\eta ^2\partial _{xx}+\frac 12\partial _{xx}\eta ^2G^i_0+\frac 12G^i_0\eta ^2G^i_0+G^i_0\eta G^i_0\eta G^i_0-\frac 12\partial _x\eta ^2\partial _x, \end{gather}

where $I_j$ and $K_j$ ( $j=0,1$ ) are the modified Bessel functions of the first and second kinds, respectively, and $P=(-\partial _{xx})^{1/2}$ . According to the properties of the modified Bessel functions (see Abramowitz & Stegun Reference Abramowitz and Stegun1965), it is not difficult to show that (2.16) is well-defined for arbitrary wavenumber $k$ . For the outer layer, the Taylor series of ${G_{\!j}}^{\!\circ}$ ( $j\geqslant 0$ ) can be directly obtained by substituting $d$ with $D$ in ${G_{\!j}}^{\!i}$ . This is feasible because both the inner and outer fluid layers share the same fundamental solution to the axisymmetric Laplace equation, along with the impermeability condition at the rigid wall. This situation differs from the two-layer interfacial wave problem in a two-dimensional geometry, where the DNOs in the upper and lower layers are formally distinct. Thus, based on the axisymmetric DNOs, it follows from the kinematic boundary condition at $r=1+\eta$ that

(2.19) \begin{equation} (1+\eta )\eta _t=G^i\xi ^i=G^\circ\xi ^\circ. \end{equation}

Recalling that $\xi =\xi ^i-\rho \xi ^\circ$ , it is straightforward to obtain

(2.20) \begin{align} \xi ^i=(G^\circ-\rho G^i)^{-1}G^\circ\xi,\quad \xi ^\circ=(G^\circ-\rho G^i)^{-1}G^i\xi. \end{align}

By applying the divergence theorem, along with (2.20), we can rewrite the Hamiltonian (2.7) in terms of the canonical variables as

(2.21) \begin{align} \mathcal{H}[\zeta ,\xi ]=&\,\frac 12\int \xi G^i\big (G^\circ-\rho G^i\big )^{-1}G^\circ\xi \,\textrm{d}{\kern0.5pt}x+\frac {B}{2}\int\! \left [\frac 12(2\eta +\eta ^2)-\ln (1+\eta )\right ]\textrm{d}{\kern0.5pt}x \nonumber\\ &\,+\int\! \left [(1+\eta )\sqrt {1+\eta _x^2}-\frac 12(2+2\eta +\eta ^2)\right ]\textrm{d}{\kern0.5pt}x. \end{align}

3.1. Shallow–shallow configuration

In the shallow–shallow case, the following scaling assumptions are made:

(3.2) \begin{equation} \partial _x=O(\epsilon ),\quad \partial _t=O(\epsilon ),\quad \eta =O(1),\quad \xi =O(1/\epsilon ),\quad B=O(1), \end{equation}

where $\epsilon =R/l$ (with $l$ representing the typical wavelength) is a small parameter reflecting the longwave approximation. Under these scaling assumptions, one expands the DNOs in powers of $\epsilon$ as

(3.3) \begin{equation} G^i_0=-\frac {1-d^2}{2}\partial _x^2+A(d)\partial _x^4+O(\epsilon ^6),\quad G^\circ_0=-\frac {1-D^2}{2}\partial _x^2 + A(D)\partial _x^4+O(\epsilon ^6), \end{equation}

where

(3.4) \begin{equation} A(d)=-\frac {1-4d^2+3d^4-4d^4\ln d}{16} \end{equation}

and $A(D)$ can be directly obtained by replacing $d$ with $D$ in (3.4). Therefore, the pseudodifferential operator ${G}^i$ in the inner layer can be expanded as

(3.5) \begin{align} G^i=&-\left (\frac {1-d^2}{2}\partial _x^2+\partial _x\eta \partial _x+\frac 12\partial _x\eta ^2\partial _x\right ) \nonumber \\ &+\left \{A(d)\partial _x^4-\frac {(1-d^2)^2}{4}\partial _x^2\eta \partial _x^2+\left [\frac {(1-d^2)^2}{8}-\frac {1-d^2}{2}\right ]\partial _x^2\eta ^2\partial _x^2\right \}+O(\epsilon ^6). \end{align}

Likewise, one can derive the expansion series of $G^\circ$ by replacing $d$ with $D$ in (3.5). By substituting the expansions of $G^{i,o}$ into (3.1) and retaining terms valid up to $O(\epsilon ^2)$ , one obtains the approximate Lagrangian,

(3.6) \begin{equation} \widetilde {\boldsymbol {L}}=\int _{\mathbb{R}}\left \{\xi \zeta _t-\frac {1}{2} \left ({T}_0+T_{2}\right )-\left [\frac {B}{2}\ln (1+\eta )-(1+\eta )\sqrt {1+\eta _x^2}\right ]\right \}\,\textrm{d}{\kern0.5pt}x, \end{equation}

where $T_0$ and $T_2$ are accurate to orders $O(\epsilon ^0)$ and $O(\epsilon ^2)$ , respectively, with the specific expressions being

(3.7) \begin{align} T_0=&\frac {1}{D^2-d^2}\Bigg [\frac {(1-d^2)(1-D^2)}{2}\xi \xi _{xx}+2\xi (\eta ^3\xi _x )_x+\frac {1}{2}\xi (\eta ^4\xi _x )_x \nonumber \\ &+ (2-D^2-d^2)\xi (\eta \xi _x)_x+\frac {6-d^2-D^2}{2}\xi (\eta ^2\xi _x )_x\Bigg ] \end{align}

and

(3.8) \begin{align} T_2=&\,\dfrac {(1-d^2)(D^2-1)(D^2+d^2-2)}{8(D^2-d^2)}\xi \big [ (2\eta -\eta ^2 )\xi _{xx}\big ]_{xx}+\dfrac {d^2A(D)+D^2A(d)}{D^2-d^2}\xi \partial _x^4\xi \nonumber\\ &\,-\dfrac {(1-d^2)(D^2-1)}{D^2-d^2}\xi (\eta ^2\xi _{xx})_{xx}-\dfrac {A(D)+A(d)}{D^2-d^2}\xi \big [(1+\eta )^2\partial _x^3\xi \big ]_x \nonumber\\ &\,-\dfrac {D^4+d^4+2(D^2+d^2)-6}{8(D^2-d^2)}\xi \big [ (\eta ^2+2\eta ) (\eta ^2\xi _{xx})_x\big ]_x \nonumber \\ &\,+\dfrac {(1-D^2)^2+(1-d^2)^2}{4(D^2-d^2)}\xi \big [ (\eta ^2+2\eta ) (\eta \xi _{xx} )_x\big ]_x. \end{align}

Minimising the approximate action $\int\! {\widetilde {\boldsymbol{L}}}\,\textrm{d}{\kern0.5pt}t$ by taking the variational derivatives with respect to $\eta$ and $\xi$ leads to governing equations, which create a closed evolution system characterised by strong nonlinearities. The specific form of the model is somewhat lengthy, and since we do not use such a high-order nonlinear model in our subsequent analyses, we have included it in Appendix A. It is important to note that the restoring forces exerted by the magnetic field and surface tension are retained precisely, without any truncation. For more details on the operation procedure of the Lagrangian perturbation method, readers can refer to § 3.1 in Xu & Wang (Reference Xu and Wang2025).

When the higher-order nonlinearity $T_2$ is omitted, a simpler model can be obtained, expressed as

(3.9) \begin{align} \eta _t+\frac {1}{1+\eta }\left \{\frac {\left [(1+\eta )^2-d^2\right ]\left [D^2-(1+\eta )^2\right ]}{2(D^2-d^2)}u\right \}_x&=0, \end{align}

(3.10) \begin{align} u_t+\left [\frac {D^2+d^2-2(1+\eta )^2}{2(D^2-d^2)}u^2\right ]_x+\frac {B\eta _x}{(1+\eta )^3}+\kappa _x&=0, \end{align}

where the governing system has been rewritten in a more revealing form by letting $u=\xi _x$ . The system (3.9)–(3.10) can be further simplified to a one-layer version,

(3.11) \begin{align} \eta _t+\frac {1}{1+\eta }\left [\dfrac {(1+\eta )^2-d^2}{2}u\right ]_x&=0, \\[-30pt]\nonumber \end{align}

(3.12) \begin{align} u_t+uu_x+\frac {B\eta _x}{(1+\eta )^3}+\kappa _x&=0, \end{align}

which was originally proposed by Cornish (Reference Cornish2018) ((4.71) and (4.72) in his thesis).

3.2. Shallow–deep configuration

In the shallow–deep case, which is a physically relevant scenario as shown in the experiment of Bourdin et al. (Reference Bourdin, Bacri and Falcon2010), one chooses the following scales:

(3.13) \begin{equation} \partial _x=O(\epsilon ),\quad \partial _t=O(\epsilon ),\quad \eta =O(1),\quad \xi =O(1/\epsilon ),\quad B=O(1),\quad D=O(1/\epsilon ). \end{equation}

Based on the scaling assumption (3.13), the DNOs can be expanded as

(3.14) \begin{align} G^{i}=&\,{-}\frac 12\big (1-d^2 \big )\partial _{xx}-\partial _{x}\eta \partial _x-\frac 12\partial _x\eta ^2\partial _x+O(\epsilon ^4), \end{align}

(3.15) \begin{align} G^{o}=&\,G_0^{o}+O\big ({1}/{(\ln (1/\epsilon ))^2}\big ), \end{align}

where

(3.16) \begin{equation} G^\circ_0=\dfrac {|P|\left [K_1(|P|D)I_1(|P|)-I_1(|P|D)K_1(|P|)\right ] }{K_1(|P|D)I_0(|P|)+I_1(|P|D)K_0(|P|)}=O\left (\frac {1}{\ln (1/\epsilon )}\right )\!. \end{equation}

By substituting the asymptotic expansions of the DNOs into Lagrangian (3.1), one arrives at the approximate Lagrangian expressed as

(3.17) \begin{align} \widetilde {\boldsymbol {L}}=&\,\int _{\mathbb{R}}\Bigg \{\xi \zeta _t+\frac 14\xi \partial _x [(1+\eta )^2-d^2 ]\partial _x\xi -\frac {B}{2}\ln (1+\eta )+(1+\eta )\sqrt {1+\eta _x^2} \nonumber \\ &\,-\frac 18\xi \partial _x\big [\big(1-d^2\big)\mathcal{T}\,\big(1+\eta \big)^2+\big(1-d^2\big)\big(1+\eta \big)^2\mathcal{T}-\big(1-d^4\big)\mathcal{T}+4\eta \mathcal{T}\,\eta \nonumber \\ &\,+2\eta \mathcal{T}\,\eta ^2+2\eta ^2\mathcal{T}\,\eta +\eta ^2\mathcal{T}\,\eta ^2\big ]\partial _x\xi \Bigg \}\,\textrm{d}{\kern0.5pt}x, \end{align}

where the pseudodifferential operator, $\mathcal{T}$ , is defined as

(3.18) \begin{equation} \mathcal{T}=-\dfrac {|P|\left [K_1(|P|D)I_0(|P|)+I_1(|P|D)K_0(|P|) \right ] }{K_1(|P|D)I_1(|P|)-I_1(|P|D)K_1(|P|)}. \end{equation}

Finally, upon inserting these expansions into (3.1), we obtain the time-evolution equations for $\eta$ and $\xi$ as follows:

(3.19) \begin{align} &(1+\eta )\eta _t+\left [\frac {(1+\eta )^2-d^2}{2}u\right ]_x=\frac {1-d^2}{4} [\mathcal{T}\,(1+\eta )^2u+(1+\eta )^2\mathcal{T}\,u \nonumber \\ &-(1+d^2)\mathcal{T}\,u ]_x+\frac 14 [4\eta \mathcal{T}\,(\eta u)+2\eta \mathcal{T}\,(\eta ^2u)+2\eta ^2\mathcal{T}\,(\eta u)+\eta ^2\mathcal{T}\,(\eta ^2u)]_x, \end{align}

(3.20) \begin{align} u_t+uu_x+\frac {B\eta _x}{(1+\eta )^3}+\kappa _x=\frac 12 [u\mathcal{T}\,(1+\eta )^2u-d^2u\mathcal{T}\,u]_x. \end{align}

When considering the outer wall boundary condition as infinitely distant from the interface, specifically with $D\rightarrow \infty$ , (3.18) reduces to

(3.21) \begin{equation} \mathcal{T}=\dfrac {|P|K_0(|P|)}{K_1(|P|)}. \end{equation}

Note that the system (3.19)–(3.20) has been rewritten in a revealing form via differentiating the dynamic boundary condition with respect to $x$ and letting $u=\xi _x$ .

A modified nonlinear model can be derived by using $\xi ^i$ instead of $\xi$ as the unknown variable. By recalling the kinematic boundary condition (2.19) and the associated operator expansions, one can express $\xi ^\circ$ in terms of $\eta$ up to the order of $O (\epsilon \ln (1/\epsilon ) )$ :

(3.22) \begin{equation} \xi ^\circ\approx \left (G^\circ_0 \right )^{-1}(1+\eta )\eta _t. \end{equation}

Substituting (3.22) into Lagrangian (3.1) and retaining terms valid up to $O(\epsilon ^2\ln (1/\epsilon ))$ , one obtains

(3.23) \begin{align} \widetilde {\boldsymbol{L}}=\int\! \Bigg \{\frac 12(G_0^\circ)^{-1}\zeta _t^2-\frac 14\xi ^i\partial _x [(1+\eta )^2-d^2]\partial _x\xi ^i-\frac {B}{2}\ln (1+\eta )+(1+\eta )\sqrt {1+\eta _x^2}\Bigg \}\,\textrm{d}{\kern0.5pt}x. \end{align}

Minimising the action of the approximate Lagrangian yields

(3.24) \begin{gather} \left (1+\eta \right )\eta _t+\left [\frac {(1+\eta )^2-d^2}{2}u\right ]_x=0, \nonumber \\[-30pt] \end{gather}

(3.25) \begin{gather} u_t+uu_x+\frac {B\eta _x}{(1+\eta )^3}+\kappa _x+\mathcal{T}\left [\frac {(1+\eta )^2-d^2}{2}u\right ]_t=0, \end{gather}

where we have defined $u=\xi _x^{i}$ as previously mentioned. Besides the simplicity in form, the system (3.24)–(3.25) has an inherent advantage over the system (3.19)–(3.20) when it comes to computing solitary waves. This is because (3.24) explicitly relates $\eta$ and $u$ for travelling waves, making the problem more concise. In addition, as demonstrated by subsequent numerical evidence, (3.24) and (3.25) constitute a suitable fully nonlinear system that will be analysed for pinch formation in the shallow–deep configuration.

3.3. Bifurcation of solitary waves: model validation

As documented in Xu & Wang (Reference Xu and Wang2025), we proposed a cubic full-dispersion model (refer to (3.39) and (3.40) in that paper) and demonstrated through numerical simulations that this model aligns well with the full Euler equations across various physical scenarios. However, the cubic full-dispersion model has two significant limitations concerning monotonic solitary waves: it is impractical to compute the branch bifurcating from a considerably large amplitude and to capture the non-monotonic feature of the speed–energy curves for depression branches. We numerically demonstrate that the strongly nonlinear model we propose in this paper can effectively address these limitations. The numerical method for solving the steady Euler equations is known as the hodograph transformation. This technique involves exchanging the independent and dependent variables and applying finite difference discretisation. In contrast, the numerical scheme for reduced models employs the standard pseudospectral method. These discretisations, along with the Newton–Raphson algorithm, have been successfully utilised to calculate axisymmetric solitary-wave solutions. The initial guess for the algorithm is obtained using the pressure vanishing method, originally proposed by Vanden-Broeck & Dias (Reference Vanden-Broeck and Dias1992). For more details on these numerical schemes, interested readers are directed to the works of Blyth & Părău (Reference Blyth and Părău2014), Doak & Vanden-Broeck (Reference Doak and Vanden-Broeck2019) and Xu & Wang (Reference Xu and Wang2025).

Figure 2.

Speed–amplitude bifurcation curves for axisymmetric depression solitary waves in a one-layer system with $d=5/11$ and two magnetic Bond numbers: $B=1.25$ (a) and $B=1.75$ (b) corresponding to the bifurcation points $c_0=0.3149$ and $c_0=0.5455$ , respectively. The solutions are obtained with different models: the full Euler equations (black dashed lines); the cubic full-dispersion model (red dotted lines); the one-layer strongly nonlinear model (blue solid lines). The inset shows the finite-amplitude bifurcation curve of the elevation branch for $B=1.75$ , computed with the strongly nonlinear model (black solid line) and the full Euler equations (black dashed line).

We first show that the strongly nonlinear model (3.11)–(3.12) for the one-layer ferrofluid problem outperforms the cubic full-dispersion model in accurately capturing large-amplitude wave profiles. It should be noted that we only consider symmetric solitary waves, with half of the solution being ‘monotonic’. It is well known that monotonic elevation and depression solitary waves have been identified in two regimes: $1\lt B\lt B_1(d)$ for elevation waves, and $B_1\lt B\lt B_2(d)$ for depression waves. Both types of waves bifurcate from the uniform stream at the global minimum of the phase speed (denoted as $c_0$ ), which can be understood from the KdV theory since these minima occur at $k=0$ . Additionally, solitary waves with opposite polarity that bifurcate from a non-zero amplitude have also been found in these regimes. For illustrative purposes, we set $d=5/11$ , leading to values of $B_1\approx 1.28$ and $B_2\approx 15.55$ . We choose two representative magnetic Bond numbers: $B=1.25$ and $B=1.75$ . The parameters selected here are consistent with those in Xu & Wang (Reference Xu and Wang2025). We first focus on the depression branch since, within the same parameter sets, the solitary wave of depression type is much narrower than the elevation wave and exhibits more pronounced nonlinearity. A comparison of the speed–amplitude bifurcation diagrams for three models – the strongly nonlinear model, the full Euler equations and the cubic full-dispersion model – is plotted in figure 2. For the one-layer problem, there is a striking quantitative agreement between the strongly nonlinear model and the Euler equations, particularly as the solitary pulse approaches the rigid bottom, a region characterised by intense nonlinearity. This scenario goes beyond the validity range of the cubic full-dispersion model. Furthermore, the strongly nonlinear model successfully captures the elevation branch that bifurcates from a considerably large finite amplitude of $S(0)\approx 2.98$ at $c=c_0$ , as depicted in the inset of figure 2(b). This capability addresses the inherent limitation of the cubic full-dispersion model, which does not allow for large-amplitude solutions. The typical profiles shown in figure 3 illustrate that the one-layer strongly nonlinear model is highly accurate at large amplitudes. Specifically, the relative difference in large-amplitude wave profiles between the Euler equations and the strongly nonlinear model is approximately $O(10^{-3})$ on this scale.

Figure 3.

(a) Typical wave profiles in the one-layer problem for $B=1.25$ and $S(0)=0.5$ . (b) Typical wave profiles for $B=1.75$ with (i) $S(0)=0.5$ and (ii) $S(0)=4.0$ .

Figure 4.

$(a)$ Speed–amplitude bifurcation curves of axisymmetric solitary waves in a two-layer system, obtained with parameters $d=5/11$ , $D=2$ and two magnetic Bond numbers: (i) $B=1.40$ and (ii) $B=3.0$ . These bifurcations emerge from $c_0=0.3542$ and $c_0=0.7921$ , respectively. The solutions are derived from the full Euler equations (black dashed lines) and the two-layer strongly nonlinear equations (blue solid lines). $(b)$ From (i) to (ii), typical wave profiles corresponding to the upward and downward triangles in $(a)$ .

Based on the numerical evidence presented, we can reasonably expect that the two-layer model, consisting of (3.9) and (3.10), will perform effectively. For the shallow–shallow configuration, we set $d=5/11$ and $D=2$ , following the choices made by Doak & Vanden-Broeck (Reference Doak and Vanden-Broeck2019). Due to the similarities, solitary waves predicted by the associated KdV equation, along with waves of opposite polarity bifurcating from finite amplitude, can also be observed in the two-layer system. We consider two typical magnetic Bond numbers: $B=1.40$ and $B=3.0$ , in the regions on either side of $B_1(d, D)\approx 1.64$ . As shown in figure 4 $(a)$ , the two-layer strongly nonlinear model successfully identifies all monotonic solitary-wave branches and aligns closely with the results from the full Euler equations. Notably, at $B=1.40$ , the elevation branch bifurcates from zero amplitude, while the depression branch bifurcates from a non-zero amplitude; this relationship is reversed for $B=3.0$ . Due to the presence of the top and bottom boundaries, the elevation and depression branches can now terminate on static profiles, which touch the upper and lower boundaries, respectively. Typical profiles are shown in figure 4 $(b)$ , where the depression and elevation solitary waves accurately approximate solutions of the Euler equations. Analogously, the newly proposed model (3.24)–(3.25) is also anticipated to effectively describe nonlinear interfacial wave phenomena in the shallow–deep configuration. We fix the value of $d=5/11$ and choose $D=12.12$ , as referenced in Bourdin et al. (Reference Bourdin, Bacri and Falcon2010), which results in $B_1(d,D)\approx 1.2840$ and $B_2(d,D)\approx 2.3004$ . We select two magnetic Bond numbers, $B=1.20$ and $B=1.50$ , to explore solution branches with varying polarities. Figure 5(a) illustrates the speed–amplitude bifurcation curves for the parameter set we examined, particularly noting the depression branches. As we push the solution towards the limiting configuration, the strongly nonlinear model demonstrates excellent agreement, even when the wave profiles nearly touch the bottom (refer to the waveform comparisons in figure 5 b with $S(0)=0.46$ , recalling that the non-dimensional thickness of the inner layer is $d\approx 0.4545$ ).

In this section, the numerical evidence regarding monotonic solitary waves shows that, in the longwave regime, the strongly nonlinear model offers a more accurate approximation of the primitive Euler equations compared with the cubic full-dispersion model. This is true regardless of whether the system consists of one layer or two. However, we must point out that strongly nonlinear models still have certain limitations. Specifically, they struggle to effectively handle wave packet solitary waves because of their weakly dispersive nature. This is where the cubic full-dispersion model comes into play. Since conducting unsteady simulations of the full Euler equations remains challenging, we will use the strongly nonlinear model as the time-dependent model for a series of numerical experiments. These experiments will focus on the stability of monotonic solitary waves and the breakup dynamics of jets.

Figure 5.

(a) Speed–amplitude bifurcation curves in the two-layer problem, with parameters set to $d=5/11$ , $D=12.12$ and two magnetic Bond numbers: $B=1.20$ and $B=1.50$ , presented from left to right. The solutions are obtained from a strongly nonlinear model for the shallow–deep case (blue solid lines) and the Euler equations (black dashed lines). Comparison of typical wave profiles, marked with a cross and a downward triangle in (a), are displayed in (b) for $B=1.20$ and $S(0)=0.46$ , and in (c) for $B=1.20$ and $S(0)=0.60$ , respectively.

4.1. Stability of solitary waves – theory

In this subsection, we analytically establish a criterion for the stability exchange of interfacial solitary waves with axisymmetric geometry, based on the two-layer full Euler equations, when subjected to longitudinal perturbations. To achieve this, we employ the asymptotic expansion method near the neutral mode, following a similar approach outlined by Kataoka (Reference Kataoka2006). To proceed, we introduce a locally confined steady solution to the dimensionless governing equations within a frame of reference moving with a constant speed $c$ along the longitudinal direction $x$ , expressed in the form

(4.1) \begin{equation} \phi ^i=\bar {\phi }^i(x-ct,r),\quad \phi ^\circ=\bar {\phi }^\circ(x-ct,r),\quad \eta =\bar {\eta }(x-ct), \end{equation}

which satisfies the magnetised Euler equations (2.2)–(2.6) exactly. The existence of solitary-wave solutions has been confirmed through numerical computations for various sets of parameters $d$ , $D$ and $B$ , as documented in Doak & Vanden-Broeck (Reference Doak and Vanden-Broeck2019). Consider an infinitesimal perturbation added to the solitary wave, expressed as

(4.2) \begin{equation} \phi ^{i,o}=\bar {\phi }^{i,o}+\hat {\phi }^{i,o}\exp (\lambda t),\quad \eta =\bar {\eta }+\hat {\eta }\exp (\lambda t), \end{equation}

where $\lambda$ is a complex constant derived from solving the linear eigenvalue problem associated with the perturbation equations. By substituting (4.2) into (2.2)–(2.6), and linearising around the fundamental state $ (\bar {\phi }^i,\bar {\phi }^\circ,\bar {\eta } )$ while imposing the decay conditions as $x\rightarrow \pm \infty$ , we derive a set of linear equations for $ (\hat {\phi }^i,\hat {\phi }^\circ,\hat {\eta } )$ ,

(4.3) \begin{align} \hat {\phi }^\circ_{rr}+\frac {1}{r}\hat {\phi }^\circ_r+\hat {\phi }^\circ_{xx}=0,\quad &\text{for}\,\,1+\bar {\eta }\lt r\lt D, \end{align}

(4.4) \begin{align} \hat {\phi }^i_{rr}+\frac {1}{r}\hat {\phi }^i_r+\hat {\phi }^i_{xx}=0,\quad &\text{for}\,\,d\lt r\lt 1+\bar {\eta }, \end{align}

subject to the boundary conditions

(4.5) \begin{align} \mathcal{L}_i [\hat {\phi }^i,\hat {\eta }]=-\lambda \hat {\eta },\quad &\text{at}\,\,r=1+\bar {\eta }, \end{align}

(4.6) \begin{align} \mathcal{L}_o [\hat {\phi }^\circ,\hat {\eta }]=-\lambda \hat {\eta },\quad &\text{at}\,\,r=1+\bar {\eta }, \end{align}

(4.7) \begin{align} \mathcal{L}_B [\hat {\phi }^i,\hat {\phi }^\circ,\hat {\eta }]=\lambda (\hat {\phi }^\circ-\hat {\phi }^i),\quad &\text{at}\,\,r=1+\bar {\eta }, \end{align}

(4.8) \begin{align} \hat {\phi }^i_r=0,\quad &\text{at}\,\,r=d, \end{align}

(4.9) \begin{align} \hat {\phi }^\circ_r=0,\quad &\text{at}\,\,r=D, \end{align}

(4.10) \begin{align} \hat {\phi }^{i,o}_x\rightarrow 0,\quad \hat {\phi }^{i,o}_r\rightarrow 0,\quad \hat {\eta }\rightarrow 0,\quad &\text{as}\,\,|x|\rightarrow \infty , \end{align}

where $\mathcal{L}_i$ , $\mathcal{L}_o$ and $\mathcal{L}_B$ are the linear operators defined as

(4.11) \begin{equation} \mathcal{L}_i=\left (\frac {\partial }{\partial r}+\bar {\eta }_x\frac {\partial }{\partial x}\right )\hat {\phi }^i+\left [\big (\!-\!c+\bar {\phi }^i_x\big )\frac {{\rm d}}{\textrm{d}{\kern0.5pt}x}-\left (\bar {\phi }^i_{rr}-\bar {\phi }^i_{xr}\frac {{\rm d}\bar {\eta }}{\textrm{d}{\kern0.5pt}x} \right )\right ]\hat {\eta }, \end{equation}

(4.12) \begin{equation} \mathcal{L}_o=\left (\frac {\partial }{\partial r}+\bar {\eta }_x\frac {\partial }{\partial x}\right )\hat {\phi }^\circ+\left [\big (\!-\!c+\bar {\phi }^\circ_x\big )\frac {{\rm d}}{\textrm{d}{\kern0.5pt}x}-\left (\bar {\phi }^\circ_{rr}-\bar {\phi }^\circ_{xr}\frac {{\rm d}\bar {\eta }}{\textrm{d}{\kern0.5pt}x}\right )\right ]\hat {\eta }\, \end{equation}

and

(4.13) \begin{align} \mathcal{L}_B=&\,\left [\left (\!-c+\frac {\partial \bar {\phi }^i}{\partial x}\right )\frac {\partial }{\partial x}+\frac {\partial \bar {\phi }^i}{\partial r}\frac {\partial }{\partial r}\right ]\hat {\phi }^i+\left [\left (\!-c+\frac {\partial \bar {\phi }^i}{\partial x}\right )\frac {\partial ^2\bar {\phi }^i}{\partial x\partial r}+\bar {\phi }^i_r\bar {\phi }^i_{rr}\right ]\hat {\eta } \nonumber \\ &\,-\!\left [\left (\!-c+\frac {\partial \bar {\phi }^\circ}{\partial x}\right )\frac {\partial }{\partial x}+\frac {\partial \bar {\phi }^\circ}{\partial r}\frac {\partial }{\partial r} \right ]\hat {\phi }^\circ-\left [\left (\!-c+\frac {\partial \bar {\phi }^\circ}{\partial x}\right )\frac {\partial ^2\bar {\phi }^\circ}{\partial x\partial r}+ \bar {\phi }^\circ_r\bar {\phi }^\circ_{rr}\right ]\hat {\eta } \nonumber \\ &\,+\frac {B\hat {\eta }}{\left (1+\bar {\eta }\right )^3}-\frac {\hat {\eta }+(1+\bar {\eta })\bar {\eta }_x\hat {\eta }_x}{(1+\bar {\eta })^2(1+\bar {\eta }_x^2)^{3/2}}-\left [ \frac {\hat {\eta }_x}{\big (1+\bar {\eta }_x^2\big )^{3/2}}\right ]_x\!. \end{align}

Equations (4.3)–(4.10) constitute an eigenvalue problem for $ (\hat {\phi }^i,\hat {\phi }^\circ,\hat {\eta } )$ , with $\lambda$ representing the corresponding eigenvalue. If this problem has a solution with $\mathrm{Re}[\lambda ]\gt 0$ , the solitary wave is considered linearly unstable under longitudinal perturbations.

4.1.1. Core solution

We begin by seeking the solution $ (\hat {\phi }^i,\hat {\phi }^\circ,\hat {\eta } )$ in the power series of $\lambda$ as follows:

(4.14) \begin{align} \hat {\phi }^{i,o}=\hat {\phi }^{i,o(0)}_{*}+\lambda \hat {\phi }^{i,o(1)}_{*}+\lambda ^2\hat {\phi }^{i,o(2)}_{*}+{\cdots},\quad \hat {\eta }=\hat {\eta }_{*}^{(0)}+\lambda \hat {\eta }^{(0)}_{*}+\lambda ^2\hat {\eta }^{(2)}_{*}+{\cdots}, \end{align}

using the asterisk subscript to specify the core solution. By substituting the expansions provided in (4.14) into (4.3)–(4.10), and then rearranging the terms according to powers of  $\lambda$ , we derive a series of equations about $\hat {\phi }^{i,o(n)}_{*}$ and $\hat {\eta }_{*}^{(n)}$ (for $n=0,1,2,{\cdots}$ ) as

(4.15) \begin{align} \hat {\phi }^{o(n)}_{*rr}+\frac {1}{r}\hat {\phi }^{o(n)}_{*r}+\hat {\phi }^{o(n)}_{*xx}=0,\quad &\text{for}\,\,1+\bar {\eta }\lt r\lt D, \end{align}

(4.16) \begin{align} \hat {\phi }^{i(n)}_{*rr}+\frac {1}{r}\hat {\phi }^{i(n)}_{*r}+\hat {\phi }^{i(n)}_{*xx}=0,\quad &\text{for}\,\,d\lt r\lt 1+\bar {\eta }, \end{align}

subject to the boundary conditions

(4.17) \begin{align} \mathcal{L}_i\big [\hat {\phi }^{i(n)}_{*},\hat {\eta }^{(n)}_{*}\big ]=G^{(n)},\quad &\text{at}\,\,r=1+\bar {\eta }, \end{align}

(4.18) \begin{align} \mathcal{L}_o\big [\hat {\phi }^{o(n)}_{*},\hat {\eta }^{(n)}_{*}\big ]=G^{(n)},\quad &\text{at}\,\,r=1+\bar {\eta }, \end{align}

(4.19) \begin{align} \mathcal{L}_B\big [\hat {\phi }^{i(n)}_{*},\hat {\phi }^{o(n)}_{*},\hat {\eta }^{(n)}_{*}\big ]=H^{(n)},\quad &\text{at}\,\,r=1+\bar {\eta }, \end{align}

(4.20) \begin{align} \hat {\phi }^{i(n)}_{*r}=0,\quad &\text{at}\,\,r=d, \end{align}

(4.21) \begin{align} \hat {\phi }^{o(n)}_{*r}=0,\quad &\text{at}\,\,r=D, \end{align}

where $G^{(n)}$ and $H^{(n)}$ are the inhomogeneous terms given by

(4.22) \begin{equation} G^{(n)}=\left \{ \begin{aligned} &0, &&\text{for}\,\,n=0,\\ &{-}\hat {\eta }^{(n-1)}_*, &&\text{for}\,\,n=1,2,{\cdots} \end{aligned}\right . \end{equation}

and

(4.23) \begin{equation} H^{(n)}=\left \{ \begin{aligned} &0, &&\text{for}\,\,n=0,\\ &\hat {\phi }^{o(n-1)}_* - \hat {\phi }^{i(n-1)}_*, &&\text{for}\,\,n=1,2,{\cdots}. \end{aligned} \right . \end{equation}

At $n=0$ , (4.15)–(4.21) form a homogeneous system. It is straightforward to demonstrate that the leading-order system has a non-trivial solution satisfying the decay conditions (4.10),

(4.24) \begin{equation} \big (\hat {\phi }^{i(0)}_{*},\hat {\phi }^{o(0)}_{*},\hat {\eta }^{(0)}_{*}\big )=\big (\bar {\phi }^i_x,\bar {\phi }^\circ_x,\bar {\eta }_x\big ), \end{equation}

which arises from the invariance of the steady Euler equations under horizontal shifts. For $n=1,2,{\cdots}$ , the system represented by (4.15)–(4.21) becomes inhomogeneous. Therefore, in order to find a solution that does not diverge exponentially as $|x|\rightarrow \infty$ , the inhomogeneous terms $G^{(n)}$ and $H^{(n)}$ must satisfy the solvability condition:

(4.25) \begin{equation} \sum _{m=1}^{n}\lambda ^m\int _{-\infty }^{+\infty }(1+\bar {\eta })\big [\big (\bar {\phi }^i_x-\bar {\phi }^\circ_x\big )G^{(m)}-\bar {\eta }_xH^{(m)}\big ]_{r=1+\bar {\eta }}dx=o (|\lambda |^n). \end{equation}

At $n=1$ , the solvability condition (4.25) is automatically satisfied, and the solution to (4.15)–(4.21) reads

(4.26) \begin{equation} \big (\hat {\phi }^{i(1)}_{*},\hat {\phi }^{o(1)}_{*},\hat {\eta }^{(1)}_{*}\big )=\big (\!-\!\bar {\phi }^i_c,-\bar {\phi }^\circ_c,-\bar {\eta }_c\big ), \end{equation}

where the homogeneous solution is omitted in (4.26) because it can be incorporated into the leading-order solution (4.24). It should be noted that both $\hat {\phi }_*^{i(1)}$ and $\hat {\phi }_*^{o(1)}$ are non-zero as $x \rightarrow \pm \infty$ , since

(4.27) \begin{align} & [\hat {\phi }^{i(1)}_*]_{x\rightarrow -\infty }= [\hat {\phi }^{i(1)}_* ]_{x\rightarrow \infty } -\frac {2}{1-d^2}\frac {\rm d}{{\rm d}c}\left (\frac {2{E}_k^i}{c}-\textit{cM}\right )\!, \end{align}

(4.28) \begin{align} & [\hat {\phi }^{o(1)}_*]_{x\rightarrow -\infty }= [\hat {\phi }^{o(1)}_* ]_{x\rightarrow \infty } -\frac {2}{D^2-1}\frac {\rm d}{{\rm d}c}\left (\frac {2 {E}_k^\circ}{c}+\textit{cM}\right )\!, \end{align}

derived from integrating $r\hat {\phi }^{o(1)}_{*x}$ and $r\hat {\phi }^{i(1)}_{*x}$ over the inner and outer fluid domains, respectively, and $M$ represents the total mass of a solitary wave, defined by

(4.29) \begin{equation} M=\int _{-\infty }^{+\infty }\dfrac {\bar {\eta }^2+2\bar {\eta }}{2}\textrm{d}{\kern0.5pt}x. \end{equation}

At $n=2$ , the solvability condition (4.25) states that

(4.30) \begin{equation} \frac {\lambda ^2}{c}\frac {d\mathcal{H}}{dc}=O(|\lambda |^3). \end{equation}

From this condition, we conclude that $( {d\mathcal{H}}/{dc})$ must be $O(\left |\lambda \right |)$ for a non-trivial solution to exist. When condition (4.30) is satisfied, a non-trivial solution $ (\hat {\phi }^{i(2)}_{*},\hat {\phi }^{o(2)}_{*},\eta ^{(2)}_* )$ to (4.15)–(4.21) exists that does not diverge exponentially as $x\rightarrow \pm \infty$ . However, this solution does not satisfy the boundary conditions (4.10). This inconsistency can be observed from (4.19) evaluated as $x\rightarrow \pm \infty$ :

(4.31) \begin{equation} c [\hat {\phi }^{o(2)}_*-\hat {\phi }^{i(2)}_* ]+(B-1-\partial _{xx})\hat {\eta }^{(2)}_*\rightarrow \hat {\phi }^{o(1)}_*-\hat {\phi }^{i(1)}_*. \end{equation}

The terms on the right-hand side of (4.31) are non-zero as $|x|\rightarrow \infty$ , due to (4.27) and (4.28), leading to a solution that conflicts with the boundary conditions (4.10). This issue persists at higher orders. The root cause of this mismatch is the failure to balance the terms involving $\lambda$ with those involving $-c{\partial }/{\partial x}$ in (4.17)–(4.19). To achieve the balance, we need to introduce a scaled coordinate with respect to $x$ and seek a solution that varies slowly with $x$ . This solution will be referred to as the far-field solution. By matching the core solution with the far-field solution, we can develop asymptotic solutions to (4.15)–(4.21) that meet the boundary condition (4.10). The current analysis is limited to the core region, while the far-field conditions (4.10) will be addressed separately in the next part.

At $n=3$ , the solvability condition (4.25) is given by

(4.32) \begin{equation} \frac {\lambda ^2}{c}\frac {d \mathcal{H}}{dc}+\lambda ^3\big \{ [ (\hat {\phi }^{o(1)}_*-\hat {\phi }^{i(1)}_* )\hat {u}^{(2)}_*]_{x\rightarrow -\infty }-\! [(\hat {\phi }^{o(1)}_*-\hat {\phi }^{i(1)}_*)\hat {u}^{(2)}_*]_{x\rightarrow \infty }\big \}=O (|\lambda |^4), \end{equation}

where the quantities in the square brackets with the subscripts $x\rightarrow \pm \infty$ are evaluated as

(4.33) \begin{equation} [\hat {u}^{(2)}_*]_{\pm \infty }=\left [\int ^1_d r\frac {\partial \hat {\phi }^{i(2)}_*}{\partial x}\textrm{d}{\kern0.5pt}x-c\hat {\eta }^{(2)}\right ]_{\pm \infty }=\left [\!-\!\int ^D_1r\frac {\partial \hat {\phi }^{o(2)}_*}{\partial x}\textrm{d}{\kern0.5pt}x-c\hat {\eta }^{(2)}\right ]_{\pm \infty }\!. \end{equation}

To derive (4.32), it is essential to carefully analyse the terms associated with $O(|\lambda |^3)$ in (4.25) by noticing that

(4.34) \begin{equation} \begin{aligned} -\bar {\eta }_x=-\hat {\eta }_*^{(0)}=G^{(1)}=\mathcal{L}_i\big [\hat {\phi }_*^{i(1)},\hat {\eta }_*^{(1)}\big ]=\mathcal{L}_o\big [\hat {\phi }_*^{o(1)},\hat {\eta }_*^{(1)}\big ], \\ \nonumber\bar {\phi }_x^\circ-\bar {\phi }_x^i=\hat {\phi }_*^{o(1)}-\hat {\phi }_*^{i(1)}=H^{(1)}=\mathcal{L}_B\big [\hat {\phi }_*^{i(1)},\hat {\phi }_*^{o(1)},\hat {\eta }_*^{(1)}\big ]. \end{aligned} \end{equation}

Then, by integrating the $O(|\lambda |^3)$ terms by parts and applying the following two equations:

(4.35) \begin{equation} \int _{-\infty }^{+\infty }dx\int _d^{1+{\bar {\eta }}}r\hat {\phi }^{i(1)}_*\boldsymbol{\nabla} ^2\hat {\phi }^{i(2)}_*\textrm{d}{\kern0.5pt}r=0,\quad \int _{-\infty }^{+\infty }dx\int ^D_{1+{\bar {\eta }}}r\hat {\phi }^{o(1)}_*\boldsymbol{\nabla} ^2\hat {\phi }^{o(2)}_*\textrm{d}{\kern0.5pt}r=0, \end{equation}

we can derive condition (4.32). In addition, integrating (4.4), multiplied by $r$ , over the entire inner fluid domain at $n=2$ yields a crucial identity:

(4.36) \begin{align} [\hat {u}_*^{(2)}]_{x\rightarrow -\infty }= [\hat {u}_*^{(2)}]_{x\rightarrow \infty } -\frac {\textrm{d}{\kern0.5pt}M}{\textrm{d}{\kern0.5pt}c}. \end{align}

A similar process applied to the outer fluid region produces the same result. It is important to note that condition (4.32) will provide the specific value of $\lambda$ , which is necessary to obtain the core solution up to the orders required to establish the criterion for the stability exchange of solitary waves. The specific values of $ [\hat {\phi }^{i(1)}_* ]_{x\rightarrow \infty }$ , $ [\hat {\phi }^{o(1)}_* ]_{x\rightarrow \infty }$ and $ [\hat {u}_*^{(2)} ]_{x\rightarrow \infty }$ , needed to determine the value of $\lambda$ , will be provided after matching the core solution with the far-field solution outlined in § 4.1.2.

4.1.2. Far-field solution

We aim to find a solution by analysing the balance between the terms involving $\lambda$ and those involving $-c\partial /\partial x$ in (4.5)–(4.7). To facilitate this, we introduce a scaled coordinate with respect to $x$ , defined as $X=|\lambda |x$ , where $|\lambda |$ is used to adjust the scale of the coordinate; we do this instead of directly using $\lambda$ since it may be complex or negative. In this formulation,  $X$ is of order $O(1)$ . Next, we seek a solution to (4.3)–(4.10) that exhibits moderate variation in both $X$ and $r$ , expressed in the following power series in terms of $\lambda$ :

(4.37) \begin{equation} \hat {\phi }^{i,o}_F(X,r)=\lambda \hat {\phi }_F^{i,o(1)}+\lambda ^2\hat {\phi }_F^{i,o(2)}+{\cdots}, \quad \hat {\eta }_F(X)=\lambda ^2\hat {\eta }_F^{(2)}+\lambda ^3\hat {\eta }_F^{(3)}+{\cdots}, \end{equation}

where the subscript ‘ $F$ ’ indicates the far-field solution. The asymptotic series begins from $O(\lambda )$ for $\hat {\phi }_F^{i,o}$ and $O(\lambda ^2)$ for $\hat {\eta }_F$ , since the core solutions exhibit non-vanishing behaviours as $|x|\rightarrow \infty$ from these orders; refer to (4.27), (4.28) and (4.36).

Substituting (4.37) into (4.3)–(4.10), and rearranging the terms of the same order in $\lambda$ , we obtain a series of equations for $\hat {\phi }_F^{i,o(n)}$ (for $n=1,2,{\cdots}$ ),

(4.38) \begin{align} \dfrac {\partial ^2\hat {\phi }^{o(n)}_F}{\partial r^2}+\dfrac {1}{r}\frac {\partial \hat {\phi }^{o(n)}_F}{\partial r}=J^{o(n)},\quad &\text{for}\,\,1\lt r\lt D, \end{align}

(4.39) \begin{align} \dfrac {\partial ^2\hat {\phi }^{i(n)}_F}{\partial r^2}+\dfrac {1}{r}\frac {\partial \hat {\phi }^{i(n)}_F}{\partial r}=J^{i(n)},\quad &\text{for}\,\,d\lt r\lt 1, \end{align}

subject to the boundary conditions

(4.40) \begin{align} \dfrac {\partial \hat {\phi }_F^{o(n)}}{\partial r}=0,\quad &\text{at}\,\,r=D, \end{align}

(4.41) \begin{align} \dfrac {\partial \hat {\phi }_F^{i(n)}}{\partial r}=0,\quad &\text{at}\,\,r=d, \end{align}

(4.42) \begin{align} \dfrac {\partial \hat {\phi }_F^{i,o(n)}}{\partial r}=K^{(n)},\quad &\text{at}\,\,r=1, \end{align}

(4.43) \begin{align} \dfrac {\partial \hat {\phi }_F^{i,o(n)}}{\partial X},\,\,\dfrac {\partial \hat {\phi }_F^{i,o(n)}}{\partial r},\,\,\big [\hat {\phi }_F^{i(n)}-\hat {\phi }_F^{o(n)}\big ]_{r=1}\rightarrow 0,\quad &\text{as}\,\,|X|\rightarrow \infty, \end{align}

where

(4.44) \begin{equation} J^{i,o(n)}= \left \{ \begin{aligned} &0,\quad &&\text{for}\,\,n=1,2\\ &-\!\dfrac {|\lambda |^2}{\lambda ^2}\dfrac {\partial ^2\hat {\phi }_F^{i,o(n-2)}}{\partial X^2},\quad &&\text{for}\,\,n\geqslant 3 \end{aligned} \right . \end{equation}

and

(4.45) \begin{equation} K^{(n)}= \left \{ \begin{aligned} &0,\quad &&\text{for}\,\,n=1,2\\ &\left (-c\dfrac {|\lambda |}{\lambda }\dfrac {d}{dX}+1\right )\hat {\eta }_F^{(n-1)},\quad &&\text{for}\,\,n\geqslant 3 \end{aligned} \right . \end{equation}

with

(4.46) \begin{align} \hat {\eta }_F^{(n)}=\frac {1}{B-1}\left (c\frac {|\lambda |}{\lambda }\frac {d}{dX}-1\right )\left [\hat {\phi }_F^{i(n-1)}-\hat {\phi }_F^{o(n-1)}\right ]_{r=1}+O(|\lambda |^2)\,,\quad \text{for}\,\,n\geqslant 2. \end{align}

At $n=1,2$ , (4.38)–(4.42) for $ (\hat {\phi }^{i(n)}_F,\hat {\phi }^{o(n)}_F )$ are homogeneous and yield a solution that is independent of $r$ , specifically $\hat {\phi }^{i,o(n)}_F=\hat {\phi }^{i,o(n)}_F(X)$ . At $n\geqslant 3$ , the system becomes inhomogeneous; thus, the terms $J^{i(n)}$ , $J^{o(n)}$ and $K^{(n)}$ must satisfy the following solvability conditions:

(4.47) \begin{equation} -\!\int _1^DrJ^{o(n)}\textrm{d}{\kern0.5pt}r=\int _d^1rJ^{i(n)}\textrm{d}{\kern0.5pt}r=K^{(n)}. \end{equation}

In particular, (4.47) at $n=3$ can be expressed as

(4.48) \begin{gather} \dfrac {\partial ^2\hat {\phi }^{i(1)}_F}{\partial X^2}-\frac {1}{c_0^2}\left (c\frac {d}{dX}-\frac {\lambda }{|\lambda |} \right )^2\big [\hat {\phi }^{i(1)}_F-\hat {\phi }^{o(1)}_F\big ]=0, \end{gather}

(4.49) \begin{gather} \dfrac {\partial ^2 \hat {\phi }^{i(1)}_F}{\partial X^2}+\dfrac {D^2-1}{1-d^2}\dfrac {\partial ^2\hat {\phi }^{o(1)}_F}{\partial X^2}=0. \end{gather}

These two equations, along with the far-field boundary conditions (4.43), determine the dependence of $\hat {\phi }^{(1)}_F$ on $X$ as

(4.50) \begin{align} &\hat {\phi }^{i(1)}_F=a_{\pm }\exp (kX)+\bar {a}_{\pm }\exp (\bar {k}X)+b_{\pm }, \end{align}

(4.51) \begin{align} &\hat {\phi }^{o(1)}_F=\dfrac {1-d^2}{1-D^2}\big [a_{\pm }\exp (kX)+\bar {a}_{\pm }\exp (\bar {k}X)\big ]+ b_{\pm }, \end{align}

where $k$ and $\bar {k}$ are given by

(4.52) \begin{equation} k=\frac {\lambda }{|\lambda |(c+c_0)},\quad \bar {k}=\frac {\lambda }{|\lambda |(c-c_0)}, \end{equation}

and $c_0$ represents the dimensionless linear longwave speed, defined as

(4.53) \begin{equation} c_0=\sqrt {\frac {(1-d^2)(D^2-1)(B-1)}{2\left (D^2-d^2\right )}}\gt c. \end{equation}

Here, $a_{\pm }$ and $\bar {a}_{\pm }$ in (4.50)–(4.51) remain undetermined and independent of the coordinates $X$ and $r$ , and those with the subscript ‘+’ are applicable for $X\gt 0$ , while those with the subscript ‘ $-$ ’ apply when $X\lt 0$ . From condition (4.43), $a_{\pm }$ and $\bar {a}_{\pm }$ must satisfy

(4.54) \begin{align} a_+=0,\quad \bar {a}_{-}=0,\quad &\text{for}\quad \text{Re}[\lambda ]\gt 0, \end{align}

(4.55) \begin{align} {a}_-=0,\quad \bar {a}_{+}=0,\quad &\text{for}\quad \text{Re}[\lambda ]\lt 0. \end{align}

Substituting (4.50)–(4.51) into (4.46) for $ n=2$ , we find that

(4.56) \begin{eqnarray} \hat {\eta }_F^{(2)}=\frac {1-d^2}{2c_0}\left [-\frac {a_{\pm }}{c+c_0}\exp ({k}X)+\frac {\bar {a}_{\pm }}{c-c_0}\exp (\bar {k}X)\right ]\!. \end{eqnarray}

Consequently, the leading-order far-field solution is obtained as (4.50)–(4.51) and (4.56).

4.1.3. Matching the core solution with the far-field solution

In this part, we focus on matching the core solution, $ (\hat {\phi }_*,\hat {\eta }_* )$ , with the far-field solution, $ (\hat {\phi }_{F},\hat {\eta }_F )$ . Within the core region, we rearrange the ordering of the far-field solution. Specifically, we expand the far-field solution $ (\hat {\phi }_F,\hat {\eta }_F )$ in a power series of $X$ (or $|\lambda |x$ ), yielding

(4.57) \begin{align} &\hat {\phi }^{i,o}_F=\lambda \big [\hat {\phi }_F^{i,o(1)}\big ]_0+\lambda ^2\big [\hat {\phi }_F^{i,o(2)}\big ]_0+\lambda |\lambda |x\left [\frac {\partial \hat {\phi }_F^{i,o(1)}}{\partial X}\right ]_0+{\cdots}, \end{align}

(4.58) \begin{align} &\hat {\eta }_F=\lambda ^2\big [\hat {\eta }_F^{(2)}\big ]_0+\lambda ^3\big [\hat {\eta }_F^{(3)}\big ]_0+ \lambda ^2|\lambda |x\left [\frac {{\rm d}\hat {\eta }_F^{(2)} }{{\rm d}X}\right ]_0+{\cdots}, \end{align}

where the subscript ‘0’ denotes evaluation at $X=0$ . By collecting terms of the same order in $|\lambda |$ , we arrive at a reordered form $ (\hat {\phi }_{F}^{i(n)^*},\hat {\phi }_{F}^{o(n)^*},\hat {\eta }_F^{(n)^*} )$ . For instance, in this new arrangement, $\hat {\phi }_F^{o(2)*}$ is of order $O(|\lambda |^2)$ and is expressed as $\lambda ^2 [\hat {\phi }^{o(2)}_F ]_0+\lambda |\lambda |x [\partial _X\hat {\phi }^{o(1)}_F ]_0$ . Next, we compare the core and far-field solutions at each order and perform the necessary matching. The matching is accomplished when the conditions

(4.59) \begin{equation} \big [\hat {\phi }^{(n)}_*\big ]_{x\rightarrow \pm \infty }=\hat {\phi }_F^{(n)^*}\quad \text{and}\quad \big [\hat {\eta }^{(n)}_*\big ]_{x\rightarrow \pm \infty }=\hat {\eta }_F^{(n)^*} \end{equation}

are satisfied. At $n=1$ , since $\hat {\phi }^{i(1)^*}_F= [\hat {\phi }^{i(1)}_F ]_0$ and $\hat {\phi }^{o(1)^*}_F= [\hat {\phi }^{o(1)}_F ]_0$ , matching is achieved if

(4.60) \begin{align} &\big [\hat {\phi }^{o(1)}_*\big ]_{x\rightarrow \pm \infty }=\frac {1-d^2}{1-D^2}\left (a_{\pm }+\bar {a}_{\pm } \right )+b_{\pm }, \end{align}

(4.61) \begin{align} &\big [\hat {\phi }^{i(1)}_*\big ]_{x\rightarrow \pm \infty }=a_{\pm }+\bar {a}_{\pm }+b_{\pm }, \end{align}

where (4.50) and (4.51) have been used. At $n=2$ , given that

(4.62) \begin{gather} \hat {\phi }^{i(2)^*}_F=\big [\hat {\phi }^{i(2)}_F\big ]_0+\frac {|\lambda |x}{\lambda }\big [\partial _X\hat {\phi }^{i(1)}_F\big ]_0,\quad \hat {\phi }^{o(2)^*}_F=\big [\hat {\phi }^{o(2)}_F\big ]_0+\frac {|\lambda |x}{\lambda }\big [\partial _X\hat {\phi }^{o(1)}_F\big ]_0, \end{gather}

the matching conditions for the velocity potentials include two types of terms: those independent of $x$ and those proportional to $x$ . These two cases must be addressed separately. The matching for the latter type is particularly crucial at this stage, which yields

(4.63) \begin{align} &\big [\hat {\phi }_{*x}^{i(2)}\big ]_{x\rightarrow \pm \infty }=\frac {a_{\pm }}{c+c_0}+\frac {\bar {a}_{\pm }}{c-c_0}, \end{align}

(4.64) \begin{align} &\big [\hat {\phi }_{*x}^{o(2)}\big ]_{x\rightarrow \pm \infty }=-\dfrac {1-d^2}{D^2-1}\left (\frac {a_{\pm }}{c+c_0}+\frac {\bar {a}_{\pm }}{c-c_0}\right )\!, \end{align}

where (4.50) and (4.51) have been used again. Next, by expressing the remaining two conditions in terms of $ [u^{(2)}_* ]_{x\rightarrow \pm \infty }$ defined by (4.33), we have

(4.65) \begin{equation} \big [u^{(2)}_*\big ]_{x\rightarrow \pm \infty }=\left [\int _{d}^{1}r\frac {\partial \hat {\phi }^{i(2)}_*}{\partial x}\textrm{d}{\kern0.5pt}r-c\hat {\eta }_*^{(2)}\right ]_{x\rightarrow \pm \infty }=\frac {1-d^2}{2}\frac {a_{\pm }-\bar {a}_{\pm }}{c_0}, \end{equation}

where (4.50) and (4.56) have been used as before. The matching conditions for $\hat {\eta }$ are automatically satisfied if (4.60)–(4.61) and (4.65) hold. In addition to the aforementioned conditions, $a_{\pm }$ and $\bar {a}_{\pm }$ must satisfy either (4.54) or (4.55), and considering the arbitrariness with respect to $b_{\pm }$ in the original set of governing (4.3)–(4.10), we can set $b_{+}=0$ .

Finally, one obtains, for $\mathrm{Re}[\lambda ]\gt 0$ ,

(4.66) \begin{align} \big [\hat {\phi }^{i(1)}_*\big ]_{x\rightarrow \infty }=\bar {a}_{+},\quad \big [\hat {\phi }^{o(1)}_*\big ]_{x\rightarrow \infty }=-\frac {1-d^2}{D^2-1}\bar {a}_{+},\quad [\hat {u}^{(2)}_*]_{x\rightarrow \infty }=-\frac {1-d^2}{2c_0}{{\bar {a}}_{+}}, \end{align}

(4.67) \begin{align} a_{-} =-\frac {c_0}{1-d^2}\frac {\rm d}{{\rm d}c}\left [\frac {c_0}{B-1}\varOmega +M \right ]\!,\quad \bar {a}_{+}=\frac {c_0}{1-d^2}\frac {\rm d}{{\rm d}c}\left [\frac {c_0}{B-1}\varOmega -M\right ]\!, \end{align}

(4.68) \begin{align} a_{+}=\bar {a}_{-}=0,\quad b_+=0,\quad b_-=-\frac {2}{D^2-d^2}\frac {\rm d}{{\rm d}c}\left [\frac {2}{c}({E}_k^i+{E}_k^\circ)\right ]\!, \end{align}

while for $\mathrm{Re}[\lambda ]\lt 0$ ,

(4.69) \begin{align} \big [\hat {\phi }^{i(1)}_*\big ]_{x\rightarrow \infty }={a}_{+},\quad \big [\hat {\phi }^{o(1)}_*\big ]_{x\rightarrow \infty }=-\frac {1-d^2}{D^2-1}\bar {a}_{-},\quad [\hat {u}^{(2)}_*]_{x\rightarrow \infty }=\frac {1-d^2}{2c_0}{{a}_{+}}, \end{align}

(4.70) \begin{align} a_{+}=\frac {c_0}{1-d^2}\frac {\rm d}{{\rm d}c}\left [\frac {c_0}{B-1}\varOmega +M\right ]\!,\quad \bar {a}_{-}=-\frac {c_0}{1-d^2}\frac {\rm d}{{\rm d}c}\left [\dfrac {c_0}{B-1}\varOmega -M\right ]\!, \end{align}

(4.71) \begin{align} a_{-}=\bar {a}_{+}=0,\quad b_+=0,\quad b_-=-\frac {2}{D^2-d^2}\frac {d}{dc}\left [\frac {2}{c}\big({E}_k^i+{E}_k^\circ\big)\right ]\!, \end{align}

with

(4.72) \begin{equation} \varOmega =\frac {2}{c}\left [\frac {2}{1-d^2}{E}_k^i-\frac {2}{D^2-1}{E}_k^\circ\right ]-\left [\frac {2}{1-d^2}+\frac {2}{D^2-1}\right ]\textit{cM}\!. \end{equation}

By substituting these relations into the solvability condition (4.32), we find that $\lambda$ is given by

(4.73) \begin{equation} \lambda =\left \{ \begin{aligned} &0,\,\,\pm \frac {1}{c\gamma }\frac {d\mathcal{H}}{dc},\quad&&\text{for}\,\frac {d\mathcal{H}}{dc}\gt 0,\\ &0,\,\,&&\text{for}\, \frac {d\mathcal{H}}{dc}\lt 0, \end{aligned} \right . \end{equation}

where

(4.74) \begin{equation} \gamma =\frac {c_0(D^2-d^2)}{2(D^2-1)}\left [\left (\frac {c_0}{B-1}\frac {d\varOmega }{dc}\right )^2+\left (\frac {{\rm d}M}{{\rm d}c}\right )^2\right ]\!. \end{equation}

The core and far-field solutions are matched at the leading order in $|\lambda |$ . The corresponding eigenfunction for $\lambda =0$ is given by (4.14) with (4.24). For $\lambda =\pm ({1}/{\gamma c}) ({{\rm d}\mathcal{H}}/{{\rm d}c})$ , the eigenfunctions are provided by (4.14) with (4.24) and (4.26) for the core solution. While the far-field solution is represented by (4.37) and (4.66)–(4.68) for Re $[\lambda ]\gt 0$ , or (4.69)–(4.71) for Re $[\lambda ]\lt 0$ , up to the first order in $|\lambda |$ . The expression for $\lambda$ in (4.73) suggests that an exchange of stability may occur at the stationary value of $\mathcal{H}$ , which is a necessary but insufficient criterion. We should emphasise that this result for radially symmetric waves extends Saffman’s finding for pure gravity waves within a two-dimensional geometry. Finally, we note that the above argument also holds true for the case when $\rho \neq 1$ .

4.2. Stability of solitary waves – numerics

Figure 6.

Speed–energy bifurcation curves in the one-layer problem computed with a fixed $d=5/11$ and varying values of $B$ : (a) the depression branch with $B=1.25$ ; (b) and (c) the depression and elevation branches with $B=1.75$ , respectively. The solution branches are computed with different models: the full Euler equations (black dashed lines); the strongly nonlinear equations (blue solid lines); the cubic full-dispersion model (red dotted lines). The inset illustrates the stationary points obtained from both the full Euler equations and the strongly nonlinear model.

In this subsection, we anticipate that the strongly nonlinear model can effectively address the stability of large-amplitude waves. We will then numerically examine the dynamics of axisymmetric solitary waves, as computed in § 3.3, in response to longitudinal perturbations based on the result discussed in § 4.1. This will involve numerically integrating the strongly nonlinear model for either a one-layer or a two-layer configuration, using slightly perturbed solitary waves as the initial conditions. According to the necessary condition for stability exchanges under superharmonic longitudinal perturbations, which we proved in § 4.1, the stability characteristics of solitary waves on a monotonic segment of the speed–energy bifurcation diagram should be consistent. Therefore, it is sufficient to examine the stability property of one typical wave situated between two successive stationary points on the speed–energy curve. In our numerical experiments, we typically choose mesh spacing and time step as $\textrm{d}{\kern0.5pt}x=0.05$ and $\textrm{d}{\kern0.5pt}t=1\times 10^{-4}$ in most computations. For convenience, we adopt a frame of reference that moves with the speed of the unperturbed solitary wave.

In the case of the one-layer ferrofluid problem, the previously proposed cubic full-dispersion model fails to capture the stationary points of the depression branches observed in the primitive Euler equations. This limitation arises because the solutions become severely nonlinear as their amplitudes increase, approaching the rigid bottom. To address this issue, we utilise a one-layer strongly nonlinear model, described by (3.11)–(3.12), to investigate the stability of the depression branch of monotonic solitary waves. Figures 6(a) and 6(b) show comparisons of the speed–energy curves derived from the one-layer strongly nonlinear model, the full Euler equations and the cubic full-dispersion model for two magnetic Bond numbers, $B=1.25$ and $B=1.75$ , respectively. The results reveal that the strongly nonlinear model not only provides quantitatively superior agreement with the full Euler equations but also accurately identifies the stationary point. For instance, with a Bond number of $B=1.25$ , the maximum Hamiltonian of the depression branch, $\mathcal{H}^{*}=0.2032$ , occurs at $c^*=0.2295$ according to the Euler equations, while our model achieves a maximum energy of $\widetilde {\mathcal{H}}^*=0.2069$ at ${\widetilde {c}\,}^*=0.2298$ , as demonstrated in the inset of figure 6(a). While there is a slight quantitative discrepancy in identifying the energy extremum point, we will demonstrate in the following numerical simulations that the strongly nonlinear model effectively qualitatively captures the stability exchange of solitary waves. For the depression branch with $B=1.25$ , superharmonic perturbations with amplitudes ranging from 1 % to 5 % of the initial wave, within the range of $0.2295\lt c\lt c_0$ , show no signs of instability. Aside from translation in the $x$ -direction, the wave structure and amplitude remain nearly unchanged from the initial wave throughout the time-evolution process. This behaviour is qualitatively consistent with that reported in § 5.2 of Xu & Wang (Reference Xu and Wang2025); therefore, the results are not explicitly shown here for brevity. In the region where $0\lt c\lt 0.2295$ , solitary waves exhibit relatively large amplitudes. However, for a depression wave with $S(0)\lt 0.6$ , the numerical scheme struggles to effectively resolve and preserve the wave steepness during the time evolution. Specifically, the numerical solution without any initial perturbation was observed to break down after a short run time. To address this numerical issue, a linear spectral stability analysis is then conducted using the method described in Appendix B, and the resulting growth rates are presented in table 1. Beyond the stationary point at $c^*=0.2295$ , we observe that the instability growth rate increases with the wave amplitude. Based on these findings, we conclude that monotonic depression solitary waves within the range $0\lt c\lt 0.2295$ are unstable. In addition, as illustrated in figure 6(c), the one-layer strongly nonlinear model demonstrates a high level of quantitative agreement with the full Euler equations for the elevation branch that bifurcates from the significantly large amplitude $S(0)=2.98$ , which the cubic full-dispersion model fails to capture. All the growth rates tested for the waves on the elevation branch in the strongly nonlinear model are less than $10^{-11}$ ; thus, the branch characterised by a monotonic speed–energy relation remains stable. However, for the depression branch with $B=1.75$ , as shown in figure 6 $(b)$ , the examinations of stability properties through unsteady simulation and linear spectral analysis follow a straightforward and similar approach; therefore, the details are omitted for brevity. Instability occurs only when $c\lt c^*=0.3286$ , as indicated by the instability growth rates listed in table 1, suggesting that the depression branch experiences an exchange of stability at the energy maximum point.

Table 1.

Instability growth rate $\lambda$ of depression solitary waves in the one-layer problem with a fixed $d=5/11$ . This analysis is conducted for two cases, $B=1.25$ and $B=1.75$ , and the corresponding regions of interest are for $c\leqslant c^*=0.2298$ and $c\leqslant c^*=0.3286$ , respectively. The midpoint value and wave energy are listed for reference.

Figure 7.

Speed–energy bifurcation curves computed with $d=5/11$ , $D=2$ and varying values of $B$ : (a i) elevation and (a ii) depression branches for $B=1.40$ ; (b) the same layout as in (a) for $B=2$ . The solution branches are computed using the full Euler equations (black dashed lines) and the two-layer strongly nonlinear model (blue solid lines).

Figure 8.

The time evolution of perturbed solitary waves with the following parameters: $d=5/11$ , $D=2$ and $B=1.40$ . (a) Snapshots of an elevation wave with $c=0.2329\gt c^*=0.1213$ , subjected to the 1 % amplitude-decreasing perturbation. The initial perturbed profile is plotted with a dashed curve, while the terminal profile at $t=5000$ is shown as a solid line. (b) Snapshots of a depression wave with $c=0.3474\gt c^*=0.2702$ , subjected to the 3 % amplitude-decreasing perturbation. The initial perturbed profile is displayed as a dashed curve, and the terminal profile at $t=4000$ is depicted by a solid line.

Next, we examine the numerical stability for monotonic solitary waves in the two-layer problem using the corresponding strongly nonlinear model described by (3.9)–(3.10). In the shallow–shallow configuration, as illustrated in figure 7, we observe stationary points in the speed–energy bifurcation curves on the full Euler equations. Similar to the one-layer scenario, the two-layer strongly nonlinear model performs effectively for the elevation branch, providing a good quantitative match at the stationary point. However, for the depression branch, small deviations from the full Euler solutions persist. Specifically, for $B=1.40\lt B_1\approx 1.64$ , the elevation branch, which emerges from zero amplitude, reaches a maximum energy of $\mathcal{H}^*=0.6285$ at $c^*=0.1213$ , as depicted in figure 7(ai). In contrast, the depression branch shown in figure 7(aii), bifurcates from finite amplitude, with the model identifying a stationary point at $\widetilde {\mathcal{H}}^*=0.2461$ and ${\widetilde {c}\,}^*=0.2774$ . These values show slight quantitative differences compared with the corresponding values of $\mathcal{H}^*=0.2409$ and $c^*=0.2702$ obtained from the full Euler equations. The numerical experiments and linear spectral analyses reveal that both the elevation and depression branches exchange stability near their stationary points, as indicated in figure 8 and table 2. For the shallow–deep configuration, the speed–energy curves for the depression branches with $B=1.20$ and $B=1.40$ are shown in figures 9(a) and 9(b), respectively. The time-evolution dynamical behaviours of these waves closely resemble those of their counterparts in the shallow–shallow case and thus are not repeated here.

Table 2.

Instability growth rate $\lambda$ of elevation $(c\leqslant c^*=0.1213)$ and depression $(c\leqslant c^*=0.2774$ ) solitary waves in the two-layer problem with $d=5/11$ , $D=2$ and $B=1.40$ . The midpoint value and wave energy are listed for reference.

Figure 9.

Speed–energy curves computed with $d=5/11$ , $D=12.12$ and varying values of $B$ : (a) the depression branches for $B=1.40$ , bifurcating from non-zero amplitude; (b) the depression branches for $B=1.50$ , bifurcating from zero amplitude. The results are derived from the full Euler equations (black dashed lines) and the strongly nonlinear model (blue solid lines).

5.1. Numerical scheme

Due to the more complex coupling between the inner and outer layers, a more efficient numerical scheme, as proposed by Choi, Goullet & Jo (Reference Choi, Goullet and Jo2011), needs to be adopted. Specifically, the system (3.24)–(3.25) can be rewritten in a more compact form as

(5.1) \begin{align} \varXi _t+\big [\big (\varXi -d^2\big )u\big ]_x=0, \end{align}

(5.2) \begin{align} \varPsi _t=-uu_x+\left [\frac {B}{2\varXi }-\kappa \right ]_x\!, \end{align}

where $\varXi =(1+\eta )^2$ and

(5.3) \begin{equation} \varPsi =u+\frac {1}{2}\mathcal{T}\, [ (\varXi -d^2 )u ]. \end{equation}

The full curvature $\kappa$ can be expressed as

(5.4) \begin{equation} \kappa =\dfrac {2}{\sqrt {4\varXi +\varXi _x^2}}\left (1-\dfrac {2\varXi \varXi _{xx}-\varXi _x^2}{4\varXi +\varXi _x^2}\right )\!, \end{equation}

which can also be approximated using its asymptotic expansion, valid up to $O(\epsilon ^2)$ , given by

(5.5) \begin{equation} \kappa \approx \frac {1}{\sqrt {\varXi }}+\frac {1}{8}\frac {\varXi _x^2}{\sqrt {\varXi }^3}-\frac {1}{2}\frac {\varXi _{xx}}{\sqrt {\varXi }}. \end{equation}

From this point forward, (5.1)–(5.2) with the curvature term $\kappa$ specified as either the full expression given in (5.4) or its asymptotic expansion in (5.5) are referred to as the full curvature model and the asymptotic slice model for the shallow–deep scenario, respectively.

It is noted that $\varPsi$ does not directly appear on the right-hand side of (5.2); therefore, it does not directly contribute to wave evolution. Next, we briefly outline the numerical scheme used for updating $u$ . The system is solved numerically at each time step using the pseudospectral method. After integrating (5.1) and (5.2) in time for $\varXi$ and $\varPsi$ with the fourth-order Runge–Kutta scheme, one needs to invert (5.3) to update $u$ . Inspired by the work of Choi et al. (Reference Choi, Goullet and Jo2011) and Wang, Hu & Wang (Reference Wang, Hu and Wang2024), we implement an iterative scheme, which can be expressed as

(5.6) \begin{equation} \big (1+\mathcal{C}\mathcal{T}\, \big )u^{(n+1)}=\varPsi -\mathcal{T}\left [\left (\frac {\varXi -d^2}{2}-\mathcal{C}\right )u^{(n)}\right ],\quad \text{for}\,n=0,1,{\cdots} \end{equation}

where the superscript $(j)$ indicates the value of $u$ at the $j$ th iteration step, and

(5.7) \begin{equation} \mathcal{C}=\frac 12\mathrm{max}\big [\varXi -d^2\big ]. \end{equation}

Taking the Fourier transform of (5.6) yields

(5.8) \begin{equation} \big (1+\mathcal{C}\widehat {\mathcal{T}}\,\big )\widehat {u}^{(n+1)}=\widehat {\varPsi }- \widehat {\mathcal{T}}\mathcal{F}\left [\left ( \dfrac {\varXi -d^2}{2}-\mathcal{C}\right )u^{(n)}\right ]\!, \end{equation}

where both $\mathcal{F}$ and $\widehat {\boldsymbol{\cdot }}$ represent the Fourier transform. The iterative scheme is repeated until

(5.9) \begin{equation} \dfrac {\big \|u^{(n+1)}-u^{(n)}\big \|_{\infty }}{\big \|u^{(n)}\big \|_{\infty }}\lt 10^{-12}, \end{equation}

is satisfied. To enhance spatial accuracy, the numerical solution is spectrally interpolated whenever the amplitude of any of the last 20 % of the modes exceeds a chosen tolerance level. This tolerance is dependent on the number of modes $N$ and typically increases as $N$ increases. In our simulations, $N$ varies from $2^{11}$ to $2^{15}$ , with tolerance level ranging from $10^{-10}$ for $N=2^{11}$ to $10^{-5}$ for $N=2^{15}$ . The time step is adaptive and decreases as $N$ increases. To ensure numerical stability, we impose a constraint on the time step $dt$ such that $k^2\,dt\lt \nu$ , where $\nu$ is between $0.1$ and $0.8$ in our calculations. The restriction on the time step arises from the linear dispersion relation for large wavenumbers (see Barannyk et al. (Reference Barannyk, Papageorgiou, Petropoulos and Vanden-Broeck2015) for details). Furthermore, we employ smooth spectral filtering as proposed by Hou & Li (Reference Hou and Li2007) to eliminate aliasing errors and smooth the highest Fourier modes.

Figure 10.

(a) The evolution of $\eta$ in the stable regime with $d=0.4$ , $D=50$ and two different magnetic Bond numbers: (i) $B=1.5$ and (ii) $B=2.0$ . (b) The corresponding maximal values of $\eta$ over time for the interval $[0, 50]$ .

Before presenting the numerical results, it is necessary to discuss the linear stability of the system. By linearising around the trivial solution $\eta =u=0$ and introducing the normal modes into the perturbations – specifically, by writing $\eta =\hat {\eta }e^{\mathrm{i}(kx-\omega t)}+\text{c.c.}$ and $u=\hat {u}e^{\mathrm{i}(kx-\omega t)}+\text{c.c.}$ – we derive the following dispersion relation:

(5.10) \begin{equation} \omega ^2=\frac {1-d^2}{2}\frac {(B-1+k^2)}{1+\frac {1-d^2}{2}\hat {\mathcal{T}}}, \end{equation}

where the symbols $k$ and $\omega$ represent the wavenumber and frequency, respectively, and

(5.11) \begin{equation} \hat{\mathcal{T}}=-\dfrac {|k|\left [K_1(|k|D)I_0(|k|)+I_1(|k|D)K_0(|k|)\right ]}{K_1(|k|D)I_1(|k|)-I_1(|k|D)K_1(|k|)}\gt 0,\quad \text{for}\,\,|k|\neq 0, \end{equation}

is a monotonically increasing function of $|k|$ , providing shortwave regularisation to the system. As indicated in (5.10), the term associated with the pseudodifferential operator $\mathcal{T}$ and the radius of the copper rod ( $0\lt d\lt 1$ ) influences only the magnitude of $\omega$ , without altering the stability characteristics. The stability of the system is primarily determined by the magnetic Bond number $B$ . The linear regime is considered stable if $B\gt 1$ for all $k$ , and unstable otherwise.

5.2. Linearly stable regime

In this part, we conduct numerical simulations in the linearly stable regime of the magnetic Bond number $B\gt 1$ . We use a grid resolution of $N=2^{10}$ and a time step of $dt=10^{-4}$ . The initial trivial solution over the domain $[0,2\pi ]$ is perturbed by a small spatially periodic disturbance, specifically $\eta =0.1\cos (x)$ and $u=0$ . The numerical solutions, derived from the full curvature model or the asymptotic slice model, exhibit a persistent oscillatory motion that appears to be quasiperiodic over time. This oscillation behaviour is qualitatively similar to that observed by Cornish (Reference Cornish2018) in a one-layer simulation. To illustrate, we use the asymptotic slice model. Figure 10(a) shows, from figure 10(a i) to figure 10(a ii), the time evolution of the interface $\eta$ over the time interval $[0, 50]$ for $B=1.5$ and $B=2.0$ , where standing waves with quasiperiodic oscillations are clearly observed. Additionally, as indicated in figure 10(b), the frequency of oscillation in the dynamics increases as the magnetic Bond number rises. Further computations are conducted for a single period with a fixed time step, as verified by the results depicted in figure 10(a).

Subsequently, we will focus on the unstable regime where $\omega$ from (5.10) no longer admits a real solution for certain wavenumbers, causing the destabilisation of the fluid system. To gain a deeper understanding of the formation of singularities, we will examine the dynamics of a monochromatic periodic wave influenced by a weak magnetic field ( $0\lt B\lt 1$ ).

5.3. Formation of infinite-slope singularities after a finite time

In this part, we conduct numerical experiments using the leading-order system derived from (5.1) and (5.2) by neglecting higher-order terms. Consequently, the governing equations form a first-order quasilinear system, which can be expressed in matrix form as

(5.12) \begin{equation} \begin{pmatrix} \varXi \\ u \end{pmatrix}_t + \begin{pmatrix} u & \varXi -d^2\\ \frac {1}{2\varXi ^2}\big (B-\sqrt {\varXi }\big ) & u \end{pmatrix} \begin{pmatrix} \varXi \\ u \end{pmatrix}_x = \begin{pmatrix} 0 \\ 0 \end{pmatrix}\!. \end{equation}

As reported by Cornish (Reference Cornish2018), when the parameter $d$ is set to $0$ , the introduction of a magnetic field still causes numerical solutions of the leading-order system to exhibit an infinite-slope singularity after a finite time. For our analyses, we set the value of $d$ to $0.395$ , a value used in the experiments of Bourdin et al. (Reference Bourdin, Bacri and Falcon2010). We then proceed to accurately evaluate the asymptotic predictions of the critical time of infinite-slope formation, as presented in Appendix C. The leading-order system is initiated with

(5.13) \begin{equation} \varXi (x,0)=1+\delta \cos (x),\quad u(x,0)=\delta \sin (x), \end{equation}

where $\delta \ll 1$ represents an initial perturbation amplitude. This choice is not necessary for the theoretical analysis but serves to clarify the presentation. Computations are stopped when the maximum magnitude of $\varXi _x$ or $u_x$ exceeds a specified threshold. The allowable maximum slopes range from 20 to 50, depending on the spatial resolution used. The timing of this occurrence is referred to as the singular time, denoted as $T_c$ . According to Cornish (Reference Cornish2018), a critical magnetic Bond number, represented by $B_c$ , was found to lie between $0.61$ and $0.64$ for the single-layer problem. For $0\lt B\lt B_c$ , the interface develops two symmetric peaks. However, as $B_c\lt B\lt 1$ , a third central peak emerges, becoming increasingly pronounced as the value of $B$ rises, ultimately dominating the instability growth. To illustrate this, we consider two specific magnetic Bond numbers, $B=0.5$ and $B=0.8$ , with $\delta =0.1$ , anticipating the emergence of two distinct infinite-slope singular structures in the leading-order system. Figures 11 and 12 depict the time-evolution snapshots of $\varXi (x,t)$ and $u(x,t)$ for $B=0.5$ and $B=0.8$ , respectively. It is observed that, after some time, the axial velocity $u(x,t)$ develops infinite slopes at two symmetric positions, as shown in figure 11(b) for $B=0.5$ , and at the centre in figure 12(b) for $B=0.8$ . At these positions, infinite slopes of $\varXi (x,t)$ also appear and likely take the form of cusps, as illustrated in figures 11(a) and 12(a) for $B=0.5$ and $B=0.8$ , respectively. Additionally, the time variations of the quantities $|\varXi _x|_{\textit{max}}$ and $|u_x|_{\textit{max}}$ are provided in figures 11 and 12, offering further evidence of shock formation. The asymptotic prediction for the blowup time, denoted as $t_c$ , is given by

(5.14) \begin{equation} t_c\approx \sqrt {\frac {2}{(1-d^2)(1-B)}}\,\ln \frac {1}{\delta }, \end{equation}

and the detailed derivation can be found in Appendix C. Meanwhile, as exemplified by $B=0.5$ , we record the critical time when the maximum value of $\varXi _x$ or $u_x$ exceeds $25$ . This is done across a broad range of values for $\delta$ , specifically, $\delta =0.1,0.08,{\cdots} ,0.001$ . Comparing these values with the asymptotic predictions from (5.14) is straightforward. The numerical results from the leading-order equations and the theoretical predictions are plotted in figure 13, showing acceptable agreement.

Figure 11.

(a,b) Snapshots of $\varXi$ and $u$ obtained by solving the leading-order equations with $B=0.5$ . (c,d) Time evolutions of the maximal values of $|\varXi _x|$ and $|u_x|$ .

Figure 12.

(a,b) Snapshots of $\varXi$ and $u$ obtained by solving the leading-order equations with $B=0.8$ . (c,d) Time evolution of the maximal values of $|\varXi _x|$ and $|u_x|$ .

Figure 13.

A comparison of the singular time between numerical solutions with the leading-order equations (black circles) and asymptotic predictions (black solid curve) for various values of $\delta$ .

5.4. Numerics for the regularised models

To generate the pinch phenomenon, we will explore two models given by (5.1)–(5.2) alongside either (5.4) or (5.5) for the shallow–deep case. In the following numerical experiments, we use the initial conditions

(5.15) \begin{equation} \varXi (x,0)=1,\quad u(x,0)=-\delta \sin (x). \end{equation}

We consider a simplified configuration in which the radius of the metallic rod is assumed to be negligible ( $d=0$ ), and the outer boundary is placed at an infinite distance ( $D\rightarrow \infty$ ). Numerical solutions are presented to demonstrate the existence of pinching singularities and to compare the solutions of the asymptotic slice model with those of the full curvature model.

In the asymptotic slice model, the initial value problem is solved with the parameters $B=0.5$ , $\delta =1$ and $\epsilon ^2=0.1$ . The time-evolution snapshots of $\varXi (x,t)$ and $u(x,t)$ are shown in figure 14. The computations are halted when the minimum of $\varXi$ (denoted by $\varXi _{\textit{min}}$ ) drops below $0.01$ . The results provide strong evidence that the interface touches the bottom in a finite time frame. This occurs as $\varXi$ approaches zero simultaneously at two symmetric points, while the corresponding velocity $u$ near these points becomes unbounded due to the fluid being pushed into droplets. In what follows, we are interested in quantifying these nearly terminal states and establishing the dominant solutions at breakup by performing a local analysis through constructing self-similar solutions. Assuming that after a time $t_s$ , the jet radius vanishes at some position $x_s$ , we look for self-similar solutions of the form

(5.16) \begin{equation} \varXi (x,t)=\tau ^{\alpha }\!f(\xi ),\quad u(x,t)=\tau ^{\beta }\!g(\xi ),\quad \xi =\dfrac {x-x_s}{\tau ^{\gamma }}, \end{equation}

where $\tau =t_s-t\gt 0$ , the constants satisfy $\alpha \gt 0$ and $\beta \lt 0$ (in accordance with the numerical evidence mentioned earlier). The balance in (5.1) as $\tau \rightarrow 0$ gives $u\sim \tau ^{\gamma -1}$ . The value of $\alpha$ emerges from a balance of the leading-order terms in (5.2) near the singular time, resulting in $\alpha =2(1-\gamma )$ . Therefore, the self-similar solutions can be expressed as

Figure 14.

Snapshots of $\varXi$ (a) and $u$ (b) with $B=0.5$ , $\delta =1$ and $\epsilon ^2=0.1$ , obtained by solving the asymptotic slice model.

(5.17) \begin{equation} \varXi (x,t)=\tau ^{2(1-\gamma )}\!f(\xi ),\quad u(x,t)=\tau ^{\gamma -1}g(\xi ), \end{equation}

where $0\lt \gamma \lt 1$ . A reasonable scaling exponent for $\gamma$ can be determined by considering the dominant balance among the higher-order terms in (5.2). As we approach a singularity and spatial scales shorten, shorter scales in real space correspond to larger values of the wavenumber $k$ in Fourier space. Hence, the behaviour of $\widehat {\mathcal{T}\,[u]}$ at large $k$ becomes significant. From the definition of (3.21), we find that

(5.18) \begin{equation} \widehat {\mathcal{T}\, \left [u\right ]}=\dfrac {|k|K_0(|k|)}{K_1(|k|)}\hat {u}(k)\sim |k|\hat {u}(k)=\widehat {P[u]}, \end{equation}

where the operator $P=(-\partial _{xx})^{1/2}$ is a Fourier multiplier, as defined in § 2.2. By applying (5.18) and conducting an order-of-magnitude analysis, it is relatively straightforward to conclude that the orders of magnitude for the non-local term and the higher capillary terms are $O(\tau ^{-2\gamma })$ and $O(\tau ^{1-4\gamma })$ , respectively. This leads to two candidates for $\gamma$ : $\gamma =2/3$ or $\gamma =3/5$ , while letting $\tau \rightarrow 0$ . However, we speculate that the pinch phenomenon is primarily driven by higher capillary terms, suggesting that the value of $\gamma$ is approximately $3/5$ . This inference will be supported by the numerical results presented later.

To numerically validate the self-similar scalings, we note that, according to (5.17), the minima of $\varXi$ and $u^{-2}$ over the $2\pi$ -periodic domain exhibit asymptotic behaviours,

(5.19) \begin{align} \mathop{min} _{0\lt x\leqslant 2\pi }\varXi \sim (t_s-t)^{2(1-\gamma )},\quad \mathop{min} _{0\lt x\leq 2\pi }{u^{-2}}\sim (t_s-t)^{2(1-\gamma )}, \end{align}

indicating that these quantities are proportional to each other. This relationship is important because we expect the plot of $\varXi _{\textit{min}}$ versus $u^{-2}_{\textit{min}}$ , with time as a parameter, to yield a straight line as we approach the singularity. Numerical results for two sets of parameters – $B=0.5$ , $\delta =1.0$ , $\epsilon ^2=0.1$ and $B=0.55$ , $\delta =0.5$ , $\epsilon ^2=0.1$ – are presented in figure 15. The plots clearly show straight lines in both cases, confirming the self-similar scalings described in (5.17).

Figure 15.

Linear relationships between $\varXi _{\textit{min}}$ and $({1}/{|u|^2_{\textit{max}}})$ with (a) $B=0.5$ , $\delta =1.0$ , $\epsilon ^2=0.1$ , and (b) $B=0.55$ , $\delta =0.5$ , $\epsilon ^2=0.1$ .

To obtain an accurate estimate of $\gamma$ , it is crucial to have a reliable method for determining the time of singularity formation, denoted as $t_s$ . The technique employed here appears to be effective, particularly because it is not affected by the specific type of ansatz being evaluated. From (5.17), we expect that as a pinch forms, the value of $u_x$ increases dramatically as $\tau ^{-1}$ , irrespective of the exponent value $\gamma$ , i.e.

(5.20) \begin{equation} \frac {1}{|u_x|_{\textit{max}}}\sim (t_s-t),\quad \text{as}\,\,t\rightarrow t_s^-. \end{equation}

The time evolution of $({1}/{|u_x|_{\textit{max}}})$ is depicted in figure 16(a) with $\varXi _{\textit{min}}\lt 0.001$ at the last computed time in the parameter space defined by $B=0.5$ , $\delta =1$ and $\epsilon ^2=0.1$ . To achieve accurate solutions with smaller values of $\varXi _{{min}}$ near the blowup, spectral interpolation is necessary to improve spatial accuracy during numerical simulations. As we approach the singularity, we apply a least-squares fitting technique to determine that $t_s\approx 1.02954$ , which is marked on the inset of figure 16(a) with an open circle. As $\tau \rightarrow 0^-$ , the exponent for $|u_x|_{\textit{max}}$ is extracted through polynomial interpolation (commonly quadratic) as shown in figure 16(c), which demonstrates good agreement with the relation (5.20).

With the estimate of $t_s$ in hand, we now proceed to estimate $\gamma$ . In figure 17, we present plots of $\ln (\varXi _{\textit{min}})$ and $\ln (u_{\textit{max}})$ against $\ln (\tau )$ , which should yield straight lines when $t$ is close to $t_s$ . The slopes of these lines are equal to $2(1-\gamma )$ and $\alpha -1$ , respectively. The fits are obtained by starting with the asymptotic behaviour described in (5.17) and noting, for example, that for the linear relation $\ln (\varXi _{\textit{min}} )=C+2(1-\alpha )\ln (t_s-t)$ (which holds in a similar way for $\ln (u_{\textit{max}})$ ), we can use the least-squares method to estimate the constant $C$ and the slope $2(1-\gamma )$ . In this analysis, we find that $\gamma \approx 0.60352$ , which shows remarkable consistency with the corresponding results from the least-squares fit of the $u_{\textit{max}}$ data. The latter yields $u_{\textit{max}}\approx (t_s-t)^{-0.39591}$ , and hence, $\gamma \approx 0.60401$ . Both estimates of $\gamma$ agree well, differing by only $5\times 10^{-4}$ . Following the same procedure, we also validate numerically the similarity scalings (5.17) when the non-local term is not present in the same parameter space. As shown in figure 18, we independently estimate values of $\gamma$ by fitting the evolutions of $\varXi _{\textit{min}}$ and $u_{\textit{max}}$ , resulting in $\gamma \approx 0.6004$ and $\gamma \approx 0.59971$ . The difference between these obtained values of $\gamma$ is approximately $7\times 10^{-4}$ . These numerical results provide further confirmation that the non-local term has a negligible quantitative effect on the similarity scaling (5.17) as $\tau \rightarrow 0^-$ . For comparison, we extract the value of $\gamma$ obtained from $\varXi _{\textit{min}}$ and $u_{\textit{max}}$ for various $B$ with the same parameters $\epsilon ^2=0.1$ and $\delta =1.0$ , starting from the initial condition (5.15). The results are summarised in table 3.

Table 3.

A comparison of $\gamma$ extracted from $\varXi _{\textit{min}}$ and $u_{\textit{max}}$ for various magnetic Bond numbers, with consistent parameters: $\epsilon ^2=0.1$ and $\delta =1.0$ . The simulations are halted when $\varXi _{\textit{min}}=10^{-3}$ .

Figure 16.

Estimate of singular times obtained by solving the asymptotic model with $B=0.5$ , $\delta =1$ and $\epsilon ^2=0.1$ : (a) time-evolution of $|u_x|_{\textit{max}}$ ; (b) time-evolution of $1/|u_x|_{\textit{max}}$ . The inset in (b) displays the blowup curve as the singularity is approached, along with a superimposed fitted curve (red dashed line) derived from a least-squares fit. The open circle marks the estimated singular time, $t_s\approx 1.02954$ . (c) Logarithmic plot of $1/|u_x|_{\textit{max}}$ evolution with the fitted curve $(t_s-t)^{0.99995}$ . The slopes (green lines) are also plotted for comparison.

To proceed, we substitute the ansatz (5.17) into the evolution system (5.1)–(5.2). This results in the following equations for the scaling functions $f$ and $g$ :

(5.21) \begin{align} -2(1-\gamma )f+\xi \gamma f'+(fg)'=0, \end{align}

(5.22) \begin{align} (1-\gamma )g+(\xi \gamma +g)g'-\frac {B}{2}\left (\frac {1}{f}\right )'+\epsilon ^2\left (\frac {{f'}^2}{8f^{3/2}}-\frac {f''}{2f^{1/2}}\right )'=0,\end{align}

where the primes stand for differentiation with respect to $\xi$ . The behaviours of the scaling functions $f$ and $g$ for large $|\xi |$ can be found by noticing that, to leading order, (5.21)–(5.22) imply

(5.23) \begin{equation} -2(1-\gamma )f+\xi \gamma f'=(1-\gamma )g+\xi \gamma g'=0\quad \text{as}\,\,|\xi |\rightarrow \infty , \end{equation}

yielding

(5.24) \begin{equation} f\sim |\xi |^{\frac {2(1-\gamma )}{\gamma }},\quad g\sim |\xi |^{\frac {\gamma -1}{\gamma }}. \end{equation}

Note that these scalings can be anticipated by recognising that, far from the singular point, the solution behaves quasistatically, to leading order, and hence, the functions $f$ and $g$ must possess algebraic forms (5.24) to effectively cancel out the time-dependent scaling factors, denoted as $\tau ^{2(1-\gamma )}$ and $\tau ^{\gamma -1}$ , respectively. The analysis also suggests that the spatial structure of the singularity corresponds to

Figure 17.

Estimate of the singular time by solving the asymptotic slice model with parameters $B=0.5$ , $\epsilon ^2=0.1$ and $\delta =1$ . (a) Log–log plot of $\varXi _{\textit{min}}$ alongside the fitted curve $(t_s-t)^{0.79295}$ (red dashed line). (b) Log–log plot of $u_{\textit{max}}$ with the fitted curve $(t_s-t)^{-0.39591}$ (red dashed line). In both (a) and (b), the corresponding slopes (green lines) are superimposed for comparison.

Figure 18.

Estimate of the singularity time derived from solving the asymptotic slice model without including the non-local term. The chosen parameters are $B=0.5$ , $\epsilon ^2=0.1$ and $\delta =1$ , leading to the singular time is $t_s\approx 1.26706$ . (a) A log–log plot of $\varXi _{\textit{min}}$ with superimposed fitted curve $(t_s-t)^{0.79920}$ . (b) Log–log plot of $u_{\textit{max}}$ , with a fitted curve $(t_s-t)^{-0.40029}$ superimposed.

(5.25) \begin{equation} \varXi \sim |x-x_s|^{\frac {2(1-\gamma )}{\gamma }}\,,\quad u\sim |x-x_s|^{\frac {\gamma -1}{\gamma }}. \end{equation}

Based on the calculated value of $\gamma \approx 0.6$ , the spatial structure at the singularity (5.25) becomes

(5.26) \begin{equation} \varXi \sim |x-x_s|^{4/3}\,,\quad u\sim |x-x_s|^{-2/3}. \end{equation}

This implies that the interface $\varXi$ touches the wall with bounded slope $\varXi _{x}$ , while exhibiting infinite curvature (meaning $\varXi _{xx}$ becomes unbounded). It is also inferred that in terms of the asymptotic behaviour of $\varXi$ at infinity, $S\sim |x-x_s|^{2/3}$ as $x\rightarrow x_s$ . However, it is important to note that the corresponding horizontal velocity, $u$ , becomes singular at the moment of touchdown due to the fluid being squeezed away from the touchdown region. At this point, the longwave models are no longer valid and should be regarded as providing an outer view of the dynamics. As shorter length scales begin to develop, all terms in the primitive equations must be considered to establish a new asymptotic regime.

For the full curvature model, it is found that this model leads to infinite-slope singularities for $\varXi$ after a finite period. Additionally, the minimum jet radius at which this occurs can be finite if $\delta$ is considerably large, or it may tend towards zero when $\delta$ is small. This scenario contrasts with the breakup observed in the asymptotic slice model, where the pinching characteristics differ; in the full curvature model, $\varXi _x$ becomes unbounded. Below, we present numerical evidence supporting these findings. Numerical solutions are terminated whenever the absolute value of $\varXi _x$ exceeds a preset tolerance level of 50, or when $\varXi _{\textit{min}}$ falls below a preset tolerance level, usually set at 0.01. We conduct refined numerical experiments with either higher allowable maximum slopes or smaller minimum radii, depending on the specific context of the study. The primary concern in this discussion is the comparison with the asymptotic slice model, particularly in demonstrating different terminal states. We start with a case that mirrors the numerical experiment used for the asymptotic slice model, specifically with parameters $B=0.5$ and $\delta =1$ . It is observed that an infinite-slope singularity forms after a finite time, before the minimum jet radius decreases to zero. When the computation is halted, we record $\varXi _{\textit{min}}=0.15214$ for this particular run. Figure 19 shows the snapshots of $\varXi$ and $u$ at $t=1.00991$ , along with the evolution of the maxima of $\varXi _x$ and $u_x$ . These plots illustrate a situation qualitatively similar to the leading-order problem addressed in § 5.3. Next, we consider the initial amplitude $\delta$ , ranging from $1.0$ to $0.1$ , with $B=0.3$ . This produces final states with successively decreasing values of the minimum jet radius. The results are shown in figure 20, illustrating the evolution of $\varXi _{\textit{min}}$ for various initial conditions. This suggests that pinching may occur when initial amplitudes are sufficiently small. The stopping condition for all runs remains as previously described.

Figure 19.

Simulations of the full curvature model with $B=0.5$ , $\epsilon ^2=0.1$ and $\delta =1$ . (a,b) Snapshots of $\varXi$ and $u$ , respectively. (c,d) Time series of maximal values of $|\varXi |_x$ and $|u|_x$ , respectively.

Figure 20.

Time evolution of $\varXi _{\textit{min}}$ for $B=0.3$ with varying inital amplitudes ( $\delta =1.0, 0.8,0.6,0.4,0.2,0.1$ from left to right).