Dynamics of solitary waves on a ferrofluid jet: the Hamiltonian framework

Abstract

The stability and dynamics of solitary waves propagating along the surface of an inviscid ferrofluid jet in the absence of gravity are investigated analytically and numerically. For the axisymmetric geometry, the problem is shown to be a conservative system with total energy as the Hamiltonian; however, one of the canonical variables differs from those in the classic water-wave problem in the Cartesian coordinate system. The Dirichlet–Neumann operator appearing in the kinetic energy is then expanded as a Taylor series, described in homogeneous powers of the surface displacement. Based on the further analysis of the Dirichlet–Neumann operator, a systematic procedure is proposed to derive reduced model equations of multiple scales in various asymptotic limits from the full Euler equations in the Hamiltonian/Lagrangian framework. Particularly, a fully dispersive model arising from retaining terms valid up to the quartic order in the series expansion of the kinetic energy, which results in quadratic and cubic algebraic nonlinearities in Hamilton's equations and henceforth is abbreviated as the cubic full-dispersion model, is proposed. By comparing bifurcation curves and wave profiles of various types of axisymmetric solitary waves among different model equations, the cubic full-dispersion model is found to agree well with the full Euler equations, even for waves of considerably large amplitudes. The stability properties of axisymmetric solitary waves subjected to longitudinal disturbances are verified with the newly proposed model. Our analytical results, consistent with Saffman's theory, indicate that in the axisymmetric cylindrical system, the stability exchange subjected to superharmonic perturbations also occurs at the stationary point of the speed-energy bifurcation curve. A series of numerical experiments for the stability and dynamics of solitary waves are performed via the numerical time integration of the model equation, and collision interactions between stable solitary waves show non-elastic features.

1. Introduction

Ferrofluids are colloidal liquids made of nanoscale ferromagnetic or ferrimagnetic particles suspended in a carrier fluid that exhibit strong magnetic properties. These physical properties and behaviour of ferrofluids can be controlled and manipulated by an external magnetic field. They are fascinating materials with a wide range of applications, including magnetic resonance imaging, dynamic loudspeakers, magnetooptic sensors, heat transfer or dissipation, etc. (see Kole & Khandekar (2021) and references therein). Since the analytic work of Bashtovoi & Krakov (1978) and the experiments carried out by Arkhipenko et al. (1980), it has been well known that magnetic fields can stabilize the Rayleigh–Plateau instability when considering a column of ferrofluid. Motivated by these pioneering works, solitary waves propagating along the surface of an axisymmetric ferrofluid jet under the action of an external magnetic field have been found and investigated (Rannacher & Engel 2006; Bourdin, Bacri & Falcon 2010; Blyth & Părău 2014; Guyenne & Părău 2016; Doak & Vanden-Broeck 2019). In addition to the aforementioned numerical and experimental results, rigorous mathematical studies have been established (Groves & Nilsson 2018; Wang & Yang 2019).

Solitary waves were initially discovered in surface water waves, where two types had ever been found. They are gravity solitary waves in shallow water and gravity-capillary solitary waves in deep water, whose bifurcation mechanisms are completely different. Gravity solitary waves in shallow water, which exist above the phase speed maximum, decay monotonically in the direction of wave propagation (precisely, a symmetric wave profile with half being monotonic) and can be approximated by the sech-squared solution of the celebrated Korteweg–de Vries (KdV) equation. However, gravity-capillary solitary waves in deep water, which exist below the phase speed minimum, feature oscillatory decaying tails and bifurcate from infinitesimal periodic waves (see Vanden-Broeck & Dias (1992) for more details). Akylas (1993) elucidated that the envelopes of gravity-capillary solitary waves can be well described by the cubic nonlinear Schrödinger (NLS) equation. Therefore, in the literature, these waves are usually called wavepacket solitary waves.

Both types of solitary waves mentioned above have already been discovered in axisymmetric ferrofluid jets. We consider an irrotational motion of an axisymmetric, inviscid and incompressible ferrofluid column, coating an electric-current-carrying metallic rod of radius $d$ mounted along the $x$-axis. The ferrofluid is surrounded by a non-magnetizable and immiscible fluid of a similar density to neglect the effect of gravity, as experimented by Bourdin et al. (2010). The current-carrying rod can induce an azimuthal magnetic field to resist the Rayleigh–Plateau instability. For the one-layer problem, i.e. the outer fluid is hydrodynamically passive, Rannacher & Engel (2006) confirmed via a linear stability analysis that the jet was indeed stabilized and derived the KdV equation describing axisymmetric weakly nonlinear disturbances in the case when the radius of the rod was negligible (namely $d=0$). They found that the features of solitary waves depend broadly on the magnetic Bond number (denoted by $B$), which is a dimensionless parameter. Specifically, if $1< B<3/2$, these are the elevation waves with a hump at the centre of waves, are depression waves with a dip at the centre if $3/2< B<9$ and are again elevation waves if $B>9$. More recently, Blyth & Părău (2014) expanded the problem into the fully nonlinear regime and conducted a numerical investigation of solitary-wave solutions to the full Euler equations in a more general setting where $d\neq 0$. Based on the technique of the hodograph transformation, they found numerically that, for $1< B< B_1(d)$, solitary waves bifurcating from zero amplitude are elevation waves, while for $B_1(d)< B\leqslant B_2(d)$, these solutions are depression waves. Moreover, the branches of depression and elevation solitary waves both bifurcating from non-zero amplitude were found within the region of $1< B< B_1(d)$ and $B_1(d)< B <2$, respectively. This is surprising since such bifurcations cannot be well elucidated by the KdV equation and have not yet been reported in free-surface water-wave problems. However, for $B>B_2(d)$, wavepacket-type solitary waves can be expected since the minimum phase speed is attained at a non-zero wavenumber, akin to capillary-gravity waves in deep water. For typical sets of parameters, two fundamental branches, elevation (a positive free-surface elevation at the centre) and depression (a negative free-surface elevation at the centre), were found to bifurcate from the phase speed minimum and terminate at a limiting configuration: static wave, trapped toroidal-shaped bubble or cusp. In the context of wavepacket-solitary waves in ferrofluid jets, a rigorous existence theorem has been established in the small-amplitude regime by Groves & Nilsson (2018). If the surrounding fluid is considered to be hydrodynamically active, Doak & Vanden-Broeck (2018) modified the numerical method to accommodate two flow domains and calculated interfacial wave solutions for the fully nonlinear problem. They discovered various types of progressive waves, including periodic, solitary and generalized solitary waves.

In contrast to computing solitary waves, which can be viewed as steady solutions in a co-moving frame of reference, there are relatively fewer studies on the stability and dynamics of solitary waves in axisymmetric ferrofluid jets. In the shallow water limit and for the magnetic Bond number less than $9$, there were theoretical descriptions of soliton collisions based on the KdV equation by Rannacher & Engel (2006). For more general situations, unsteady simulations for ferrofluid jets were considered by Guyenne & Părău (2016), who derived the convergent Taylor series expansion of the Dirichlet–Neumann operator (DNO) in the axisymmetric geometry for $d=0$. Unlike the hodograph transformation, their numerical approach reduces the problem to a lower-dimensional computation involving surface variables only and is easily combined with the fast Fourier transform. Those authors used the proposed time-dependent algorithm to simulate head-on and overtaking collisions between two KdV-type solitons for $B\leq 9$. Nevertheless, to the best of our knowledge, the stability properties of solitary waves on ferrofluid jets, including both monotonic and oscillatory types, have not yet been known, as well as the dynamics of wavepacket solitary waves (namely when $B>B_2(d)$), which are one of the focuses of the present work.

To understand the stability properties of solitary waves, the concept of superharmonic perturbation needs to be invoked. This is because, for solitary waves, all the longitudinal perturbations have the same horizontal scale as the fundamental waves and, hence, are superharmonic. In free-surface water waves, it was Longuet-Higgins (1978) who first considered the superharmonic instability of periodic gravity waves. Subsequently, Tanaka (1983) found numerically that the superharmonic instability occurs at the first stationary value of the total wave energy, kinetic plus potential, as a function of wave speed. He also confirmed that the same principle holds for solitary waves (Tanaka 1986). However, Saffman (1985) made another significant analytical work based on Zakharov's Hamiltonian formulation of water waves, which demonstrated that the stability exchange of periodic waves in deep water subjected to superharmonic perturbations could only manifest at stationary points of total energy, a necessary but insufficient condition (a similar argument can be found from Zufiria & Saffman (1986) for the finite-depth case). As pointed out in his original paper, the argument can be easily generalized to include the surface tension term without essential modifications. However, the results discussed above pertain to solitary and periodic waves in the classical two-dimensional water-wave problem. In this paper, we aim to extend Saffman's theory to include solitary waves on the surface of an axisymmetric ferrofluid jet. To achieve this, we must first determine the Hamiltonian of the jet system and identify the corresponding canonical variables.

After obtaining the Hamiltonian structure of the axisymmetric system, as a first step towards a comprehensive understanding of axisymmetric solitary waves, we derive several reduced models based on the Hamiltonian/Lagrangian perturbation theory. Among all the models, the cubic full-dispersion model (see (3.39)–(3.40)) is the most accurate and efficient. Using the Taylor series expansion of the DNO, the model can be achieved by truncating the kinetic energy at the fourth-order nonlinearity. It should be emphasized that to ensure self-adjointness, the DNO is not defined directly by using the kinematic boundary condition and, hence, differs from that of Guyenne & Părău (2016). This is ascribed to the particularity of the cylindrical coordinate system (see § 2.4 for details). Since the cubic full-dispersion model retains the complete dispersion relation, it can provide a quantitatively good approximation to both monotonic and oscillatory solitary waves in the full Euler equations, including wave profiles of small- and moderate-amplitude solutions and bifurcation diagrams near bifurcation points. Subsequently, using the cubic full-dispersion model, the analytical verification of the superharmonic instability for axisymmetric solitary waves is conducted via a perturbation method for small eigenvalues, achieving the same result as Saffman's theory. Guided by these theoretical findings, we then employ a time integration technique to numerically study the stability and collision of these waves.

The rest of the paper is structured as follows. In § 2, we present the mathematical formulation of the problem, discuss its Hamiltonian structure and define the axisymmetric DNO along with the Taylor expansion in recursion form. In § 3, we propose a systematic procedure based on the Hamiltonian/Lagrangian framework to derive a series of weakly and strongly nonlinear dispersive models. The numerical results for various types of solitary waves are presented in § 4, which includes the comparison of bifurcations between the Euler equations and the models. Building on the cubic full-dispersion model, § 5 employs an asymptotic method to derive a necessary condition for the stability exchange of axisymmetric travelling waves subjected to superharmonic perturbations. This section also explores the dynamics of solitary waves, focusing on their long-term evolution to identify stability properties and the pairwise collisions between stable ones. Finally, § 6 provides conclusions and further remarks.

2. Mathematical formulation

2.1. Governing equations

We consider an axisymmetric column of ferrofluid with constant density $\rho$ and magnetic susceptibility $\chi$, coating a metallic rod of radius $d$. The axisymmetric cylindrical coordinate system $(x,r)$ is introduced, such that the $x$-axis coincides with the rod's central axis, and $r$ is the radial coordinate (see the sketch of the system in figure 1). In the basic configuration, the jet surface is a circular cylinder of radius $R$. An azimuthal magnetic field is induced by the metallic rod carrying a current ${I}$:

(2.1)\begin{equation} \boldsymbol{B}=\frac{\mu_0 {I}}{2{\rm \pi} r}\boldsymbol{e}_{\theta}, \end{equation}

where $\mu _0$ is the magnetic permeability of the vacuum and $\boldsymbol {e}_{\theta }$ is the unit vector in the clockwise azimuthal direction. The assumption of an axisymmetric surface results in the decoupling of the magnetic problem from the hydrodynamic problem, and therefore, the effect of the magnetic field appears in the balance of the normal stress at the surface in an implicit way. The ferrofluid is assumed inviscid and incompressible, which flows irrotationally; therefore, there exists a velocity potential $\phi$ such that the velocity field of fluid can be expressed as $\boldsymbol {\nabla }\phi$, where $\phi$ satisfies the axisymmetric Laplace equation:

(2.2)\begin{equation} \frac{\partial^2 \phi}{\partial r^2}+\frac{1}{r}\frac{\partial\phi}{\partial r} +\frac{\partial^2\phi}{\partial x^2}=0. \end{equation}

We require no normal flow at the axial rod; that is,

(2.3)\begin{equation} \frac{\partial\phi}{\partial r}=0,\quad\ \textrm{at\ $r=d$}.\end{equation}

At the surface of the jet located at $r=S(x,t)$, the kinematic boundary condition demands

(2.4)\begin{equation} \frac{\partial S}{\partial t}=\frac{\partial\phi}{\partial r}-\frac{\partial S}{\partial x}\frac{\partial\phi}{\partial x}.\end{equation}

In the considered case when the effect of gravity can be neglected, the surface tension and magnetic field jointly provide restoring forces for free-surface fluctuations of the ferrofluid. Thus, under the assumption of a linear magnetization law (see § 5.1 of Rosensweig 1985), the dynamic boundary condition arises from the Young–Laplace equation, which reads

(2.5)\begin{equation} \frac{\partial \phi }{\partial t}+\frac12|\boldsymbol{\nabla}\phi|^2+\frac{\sigma}{\rho}\kappa-\frac{\mu_0\chi I^2}{4{\rm \pi}^2\rho}\frac{1}{2S^2}=0,\end{equation}

where $\sigma$ stands for the surface tension coefficient and $\kappa$ represents the mean curvature in the cylindrical coordinates with radial symmetry, namely

(2.6)\begin{equation} \kappa=\frac{1}{S(1+S_x^2)^{1/2}}-\frac{S_{xx}}{(1+S_x^2)^{3/2}}.\end{equation}

Equations (2.2)–(2.5) form a closed system for axisymmetric free-surface motions of a ferrofluid jet under an externally applied magnetic field. All physical variables in the system can be non-dimensionalized as

(2.7a–e)\begin{equation} r=R\bar{r},\quad x=R\bar{x},\quad S=R\bar{S},\quad t=\left(\frac{\rho R^3}{\sigma}\right)^{1/2}\bar{t},\quad\phi=\left(\frac{\sigma R}{\rho}\right)^{1/2}\bar{\phi}. \end{equation}

Rewriting the system under consideration with the scaling (2.7a–e) and then dropping bars for ease of notation, we can express the free surface as the perturbation about $r=1$, namely $S=1+\eta$, where $\eta$ is the displacement away from equilibrium. Therefore, the kinematic and dynamic boundary conditions become

(2.8)\begin{equation} \frac{\partial\eta}{\partial t}=\frac{\partial\phi}{\partial r}-\frac{\partial\eta}{\partial x}\frac{\partial\phi}{\partial x}\end{equation}

and

(2.9)\begin{equation} \frac{\partial\phi }{\partial t}+\frac12|\boldsymbol{\nabla}\phi|^2+\kappa-1-\frac{B}{2}\left[\frac{1}{(1+\eta)^2}-1\right]=0,\end{equation}

where $B={\mu _0\chi I^2}/{4{\rm \pi} ^2\sigma R}$ is called the magnetic Bond number and $\kappa$ can be recast as

(2.10)\begin{equation} \kappa=\frac{1}{(1+\eta)(1+\eta_x^2)^{1/2}}-\frac{\eta_{xx}}{(1+\eta_x^2)^{3/2}}.\end{equation}

Figure 1. Schematic of the problem in a cylindrical coordinate system.

2.2. Linear theory

To obtain the dispersion relation, we first linearize the whole system around the trivial uniform stream solution $\eta =\phi =0$ and then seek solutions in the form $\exp ({\mathrm {i}(kx-\omega t)})$ with wavenumber $k$ and frequency $\omega$. By requiring the solution of the linearized system to be non-trivial, we obtain, after some algebra, the linear dispersion relation:

(2.11)\begin{equation} c^2=\frac{\omega^2}{k^2}=\frac{1}{|k|}\left[\frac{{\rm K}_1(|k|d){\rm I}_1(|k|)-{\rm I}_1(|k|d){\rm K}_1(|k|)}{{\rm K}_1(|k|d){\rm I}_0(|k|)+{\rm I}_1(|k|d){\rm K}_0(|k|)}\right](k^2-1+B),\end{equation}

where ${\rm I}_j$ and ${\rm K}_j$ ($j=0,1$) represent modified Bessel functions of the first and second kind, respectively, and $c$ is called the phase speed. Note that all solutions are stable if $B>1$, which means a magnetic field of sufficient strength can linearly stabilize the Rayleigh capillary mode responsible for jet breakup under normal conditions. In the long-wave limit ($k\rightarrow 0$), (2.11) can be approximated by

(2.12)\begin{equation} c=c_0\left[1+\left(\frac{1}{2(B-1)}+\frac{A(d)}{1-d^2}\right)k^2+O(k^4)\right], \end{equation}

where

(2.13a,b)\begin{equation} c^2_0=\frac{(1-d^2)(B-1)}{2}, \quad A(d) ={-}\frac{1-4d^2+3d^4-4d^4\ln d}{16}. \end{equation}

The coefficient sign of the quadratic nonlinearity determines whether the curve is concave or convex. It is also noted that $c\sim \sqrt {|k|}$ as $k\rightarrow \infty$. Due to these asymptotic behaviours, a critical magnetic Bond number $B_2(d)$ can be defined,

(2.14)\begin{equation} B_2(d)=\frac{9-12d^2+3d^4-4d^4\ln d}{1-4d^2+3d^4-4d^4\ln d}, \end{equation}

such that the phase speed is a monotonically increasing function of the wavenumber for $1< B< B_2(d)$, while the phase speed minimum occurs at a non-zero wavenumber when $B>B_2(d)$. This threshold, also given by Doak & Vanden-Broeck (2019), distinguishes the KdV- and wavepacket-type solitary waves in the weakly nonlinear theory. This phenomenon is similar to capillary-gravity waves, where the Bond number is defined to measure the importance of surface tension to gravity, and its critical value is 1/3 (see Hunter & Vanden-Broeck (1983) for more details). For small $d$, $B_2(d)=9+24d^2+O(d^4\ln d)$, coincident with the results given by Rannacher & Engel (2006) and Blyth & Părău (2014).

2.3. Hamiltonian formulation

The classic water-wave problem can be written in a canonical form in the sense of Zakharov (1968), who showed that the Hamiltonian is simply the total energy with the surface displacement and surface velocity potential being canonical conjugates. In the subsequent analyses, we demonstrate that for the ferrofluid jet problem studied here, the Hamiltonian functional is also the total energy; however, the canonical pair differs from that of the classical water-wave problem. For axisymmetric ferrofluid jets, the total energy reads

(2.15)\begin{align} \mathcal{H}&=\frac12\int\int_d^{1+\eta}|\boldsymbol{\nabla}\phi|^2r\,\mathrm{d}r\,\mathrm{d}{\kern 0.9pt}x+\frac{B}{2}\int\left[\frac12(2\eta+\eta^2)-\ln(1+\eta)\right]{\rm d}{\kern 0.9pt}x\nonumber\\ &\quad+\int\left[(1+\eta)\sqrt{1+\eta_x^2}-\frac12(2+2\eta+\eta^2)\right]{\rm d}{\kern 0.9pt}x, \end{align}

where the first term on the right-hand side of the equals sign is the kinetic energy, and the next two terms are responsible for the potential energy due to the magnetic field and surface tension, respectively. Next, we will show that the two canonical variables are $\zeta ={(1+\eta )^2}/{2}$ and $\xi =\phi (x,1+\eta,t)$. Further details regarding the rationale behind our selection of canonical variables can be found in Appendix A. We perturb the Hamiltonian by $\delta \xi \neq 0$ and $\delta \zeta =0$. It follows from

(2.16a,b)\begin{equation} \left.\delta\xi=\delta\phi\right|_{r=1+\eta}+\frac{\partial\phi}{\partial r}\delta\eta\quad\text{and}\quad0=\delta\zeta=(1+\eta)\delta\eta \end{equation}

that $\delta \xi =\delta \phi |_{r=1+\eta }$. Obviously, the variation in $\phi$ contributes only to the variation of the kinetic energy, denoted by $E_K$. The variation of $E_{K}$ with respect to $\xi$ is then calculated as

(2.17)\begin{align} \delta E_K&=\delta\left(\frac12\int\int_d^{1+\eta}\boldsymbol{\nabla}\phi\boldsymbol{\cdot} r\boldsymbol{\nabla}\phi\,{\rm d}r\,{\rm d}{\kern 0.9pt}x\right)\nonumber\\ &=\int\left(r\delta\phi\frac{\partial\phi}{\partial\boldsymbol{n}}\right)_{r=1+\eta}\sqrt{1+\eta_x^2}\,{\rm d}{\kern 0.9pt}x=\int\delta\xi(1+\eta)\left(\frac{\partial\phi}{\partial\boldsymbol{n}}\right)_{r=1+\eta}\sqrt{1+\eta_x^2}\,{\rm d}{\kern 0.9pt}x, \end{align}

where $\boldsymbol {n}=(-\eta _x,1)^\top /\sqrt {1+\eta _x^2}$ is the exterior unit normal vector to the free surface and the divergence theorem in the cylindrical coordinate system has been used. It follows immediately that, at $r=1+\eta$,

(2.18)\begin{equation} \frac{\partial\phi}{\partial\boldsymbol{n}}\sqrt{1+\eta_x^2}=\frac{\partial\phi}{\partial r}-\frac{\partial\eta}{\partial x}\frac{\partial\phi}{\partial x}.\end{equation}

Hence,

(2.19)\begin{equation} \frac{\delta \mathcal{H}}{\delta\xi}=\frac{\delta E_K}{\delta\xi}=(1+\eta)\left(\frac{\partial\phi}{\partial r}-\frac{\partial\eta}{\partial x}\frac{\partial\phi}{\partial x}\right)=\zeta_t,\end{equation}

equivalent to the kinematic boundary condition at the free surface. Then, we fix $\xi$ and vary $\zeta$. The relation between $\zeta$ and $\eta$ gives $\delta \zeta =(1+\eta )\delta \eta$. Meanwhile, it can be noticed that

(2.20)\begin{equation} 0=\left.\delta\xi=\delta\phi \right|_{r=1+\eta}+\frac{\partial\phi}{\partial r}\delta\eta\quad\Rightarrow\quad \left.\delta\phi\right|_{r=1+\eta}={-}\frac{\partial\phi}{\partial r}\delta\eta.\end{equation}

Thus, the variation of the kinetic energy reads

(2.21)\begin{align} \delta E_K&=\frac12\int\left.\left(r|\boldsymbol{\nabla}\phi|^2\right)\right|_{r=1+\eta}\delta\eta\,{{\rm d}{\kern 0.9pt}x} +\int\left.\left(\delta\phi r\frac{\partial\phi}{\partial\boldsymbol{n}}\right)\right|_{r=1+\eta}\sqrt{1+\eta_x^2}\,{{\rm d}{\kern 0.9pt}x}\nonumber\\ &=\frac12\int\left(|\boldsymbol{\nabla}\phi|^2\right)_{r=1+\eta}\delta\zeta\,{{\rm d}{\kern 0.9pt}x}+\int\left(\delta\phi \right)_{r=1+\eta}(1+\eta)\eta_t\,{{\rm d}{\kern 0.9pt}x}\nonumber\\ &=\int\left(\frac12|\boldsymbol{\nabla}\phi|^2-\frac{\partial\phi}{\partial r}\eta_t\right)_{r=1+\eta}\delta\zeta\,{{\rm d}{\kern 0.9pt}x}. \end{align}

However, the variation of the potential energy (denoted by $E_P$) can be expressed as

(2.22)\begin{align} \delta E_P&=\frac{B}{2}\int\left[(1+\eta)\delta\eta-\frac{\delta\eta}{1+\eta}\right]{\rm d}{\kern 0.9pt}x + \int\left[{\sqrt{1+\eta_x^2}}-\left(\frac{(1+\eta)\eta_x}{\sqrt{1+\eta_x^2}}\right)_x-(1+\eta)\right] \delta\eta\,{\rm d}{\kern 0.9pt}x\nonumber\\ &=\frac{B}{2}\int\left[1-\frac{1}{(1+\eta)^2}\right]\delta\zeta\,{\rm d}{\kern 0.9pt}x+\int\left[\frac{1}{(1+\eta)\sqrt{1+\eta_x^2}}-\frac{\eta_{xx}}{(1+\eta_x^2)^{3/2}}-1\right]\delta\zeta\,{\rm d}{\kern 0.9pt}x. \end{align}

Using the chain rule when necessary, a direct calculation yields

(2.23)\begin{align} -\frac{\delta \mathcal{H}}{\delta\zeta}&={-}\frac12|\boldsymbol{\nabla}\phi|^2+\frac{\partial\phi}{\partial r}\eta_t+\frac{B}{2}\left[\frac{1}{(1+\eta)^2}-1\right]\nonumber\\ &\quad -\left[\frac{1}{(1+\eta)\sqrt{1+\eta_x^2}}-\frac{\eta_{xx}}{(1+\eta_x^2)^{3/2}}-1\right] =\xi_t. \end{align}

Equations (2.19) and (2.23) confirm that $\zeta$ and $\xi$ are the canonical variables satisfying Hamilton's equations:

(2.24a,b)\begin{equation} \zeta_t=\frac{\delta\mathcal{H}}{\delta\xi},\quad\xi_t={-}\frac{\delta\mathcal{H}}{\delta\zeta},\end{equation}

equivalent to the kinematic and dynamic boundary conditions on the free surface.

2.4. Dirichlet–Neumann operator

The Dirichlet–Neuamnn operator (DNO), which connects the Dirichlet data and Neumann data on the boundary, plays an essential role in boundary problems of potential theory. It was initially proved by Coifman & Meyer (1985) that if the Lipschitz-norm of the surface displacement is smaller than a certain constant, then the DNO is an analytic operator with a convergent Taylor series. The Taylor expansion was obtained by Craig & Sulem (1993) from a recursion formula for the two-dimensional water-wave problem and extended to the three-dimensional case by Craig, Schanz & Sulem (1997). This description of the DNO, involving only the variables on the free surface, achieves a dimension reduction and, hence, is numerically efficient for spatially periodic water waves.

Guyenne & Părău (2016) first introduced the DNO in the problem of axisymmetric ferrofluid jets, assuming a thin current-carrying rod ($d=0$) to simplify the derivation of its Taylor series expansion. They defined the operator by directly using the kinematic boundary condition. We should point out that in the axisymmetric case, the traditional definition of the DNO is not a self-adjoint operator. Indeed, this definition is acceptable for numerical simulations, especially when we directly substitute the expansion into the kinematic and dynamic boundary conditions; however, it may cause inconvenience while using the Hamiltonian structure of the problem for investigations. This can be understood from the kinetic energy in (2.15), where the determinant introduced by the polar coordinate transformation needs to be carefully considered to ensure some ‘symmetry’ of the operator. In the present paper, we derive the Taylor expansion of the DNO for a more general case, namely $d\neq 0$, and therefore, both the first and second kinds of Bessel functions need to be invoked. However, more importantly, we define the DNO in a self-adjoint manner, differing from that of Guyenne & Părău (2016).

In the axisymmetric geometry, the DNO for the velocity potential, denoted by $G$, can be defined by solving the following boundary value problem:

(2.25)\begin{equation} \left.\begin{array}{@{}ll@{}}\displaystyle\frac{\partial^2\phi}{\partial r^2}+\frac{1}{r}\frac{\partial\phi}{\partial r}+\frac{\partial^2\phi}{\partial x^2}=0 & \textrm{for\ $d< r<1+\eta$},\\[10pt] \displaystyle\frac{\partial\phi}{\partial r}=0, & \textrm{at\ $r=d$},\\[9pt] \phi=\xi, & \textrm{at\ $r=1+\eta$}, \end{array}\right\} \end{equation}

and then returning

(2.26)\begin{equation} G(\eta)\xi=r(\phi_{r}-\eta_x\phi_{x})\end{equation}

at the surface $r=1+\eta$. Applying the separation of variables to the radial Laplace equation with the form $\phi (x,r)=f(r){\rm e}^{\mathrm {i}kx}$ yields

(2.27)\begin{equation} \bar{k}^2f''+\bar{k} f'-\bar{k}^2f=0,\end{equation}

where $\bar {k}=|k|r$ and the primes represent the derivations with respect to $\bar {k}$. Equation (2.27) is the modified Bessel differential equation of order zero, whose solution is the linear combination of the modified Bessel functions of the first and second kinds, namely

(2.28)\begin{equation} f(|k|r)=c_1{\rm I}_0(|k|r)+c_2{\rm K}_0(|k|r)\quad\textrm{for\ $k\neq0$}, \end{equation}

where $c_1$ and $c_2$ are arbitrary constants determined by boundary conditions. It is noted that for $k=0$, the Bessel equation reduces to $f''+f'/r=0$, and the solution reads $f(r)=c_1+c_2\ln r$. The impermeability condition at $r=d$ gives

(2.29)\begin{equation} f(|k|r)=c_3\left[{\rm K}_1(|k|d){\rm I}_0(|k|r)+{\rm I}_1(|k|d){\rm K}_0(|k|r)\right],\end{equation}

where $c_3$ is an arbitrary constant, and ${\rm I}_0'={\rm I}_1$ and ${\rm K}_0'=-{\rm K}_1$ have been used in the calculation. We assume the DNO can be expanded as a convergent Taylor series

(2.30)\begin{equation} G(\eta,d)=\sum_{j=0}^\infty G_{j}(\eta,d), \end{equation}

where $G_{j}$ is homogeneous of degree $j$ in $\eta$. We suppress the dependence of $G$ and $G_j$ on $d$ for simplicity of notation. Ignoring the constant in (2.29) and recalling the definition of the pseudo-differential operator, one obtains, at the free surface,

(2.31)\begin{align} &G(\eta)\left\{\left[{\rm K}_1(|k|d){\rm I}_0(|k|r)+{\rm I}_1(|k|d){\rm K}_0(|k|r)\right]{\rm e}^{\mathrm{i}kx}\right\}\nonumber\\ &\quad =r\left(\partial_r-\eta_x\partial_x\right)\left\{\left[{\rm K}_1(|k|d){\rm I}_0(|k|r)+{\rm I}_1(|k|d){\rm K}_0(|k|r)\right]{\rm e}^{\mathrm{i}kx}\right\}. \end{align}

Substituting $r=1+\eta$ into identity (2.31) and expanding it around $\eta =0$ yields

(2.32)\begin{align} &\left(\sum_{j=0}^\infty G_{j}(\eta)\right)\left\{\sum_{n=0}^\infty\frac{(|k|\eta)^n}{n!}\left[{\rm K}_1(|k|d){\rm I}_0^{(n)}(|k|)+{\rm I}_1(|k|d){\rm K}_0^{(n)}(|k|)\right]\right\}{\rm e}^{\mathrm{i}kx}\nonumber\\ &\quad =(1+\eta)\left\{\sum_{n=0}^\infty\frac{(|k|\eta)^n|k|}{n!}\left[{\rm K}_1(|k|d){\rm I}_0^{(n+1)}(|k|)+{\rm I}_1(|k|d){\rm K}_0^{(n+1)}(|k|)\right]\right.\nonumber\\ &\qquad \left.-\sum_{n=0}^\infty\frac{\mathrm{i}\eta_x(|k|\eta)^nk}{n!}\left[{\rm K}_1(|k|d){\rm I}_0^{(n)}(|k|)+{\rm I}_1(|k|d){\rm K}_0^{(n)}(|k|)\right]\right\}{\rm e}^{\mathrm{i}kx}, \end{align}

where the superscript with parentheses stands for the derivative. By identifying terms of the same degree in $\eta$, we can obtain an expansion of the operator $G(\eta )$, with each term in the series defined recursively. The first term reads

(2.33)\begin{equation} G_{0}{\rm e}^{\mathrm{i}kx}=\frac{|k|\left[{\rm K}_1(|k|d){\rm I}_1(|k|)-{\rm I}_1(|k|d){\rm K}_1(|k|)\right]}{{\rm K}_1(|k|d){\rm I}_0(|k|)+{\rm I}_1(|k|d){\rm K}_0(|k|)}\,{\rm e}^{\mathrm{i}kx},\end{equation}

and therefore, we can write $G_{0}$ as a pseudo-differential operator as

(2.34)\begin{equation} G_{0}=\frac{|P|\left[{\rm K}_1(|P|d){\rm I}_1(|P|)-{\rm I}_1(|P|d){\rm K}_1(|P|)\right]}{{\rm K}_1(|P|d){\rm I}_0(|P|)+{\rm I}_1(|P|d){\rm K}_0(|P|)},\end{equation}

where $P=-\mathrm {i}\partial _x$. Based on the properties of the modified Bessel functions of the first and second kinds (see Abramowitz & Stegun 1964), it is not difficult to show that (2.33) is well defined for arbitrary wavenumber $k$. It is remarked that if we consider a very thin current-carrying rod as the limiting configuration (i.e. $d\rightarrow 0$), then $G_0$ reduces to

(2.35)\begin{equation} G_{0}=\frac{|P|{\rm I}_1(|P|)}{{\rm I}_0(|P|)}\,, \end{equation}

due to the fact that ${\rm I}_1(|k|d)\sim {|k|d}/{2}$ for small $d$, which is exactly the same as derived by Guyenne & Părău (2016). For $j=1$,

(2.36)\begin{equation} G_{1}(\eta)+G_{0}\eta G_{0}=\eta|k|^2\frac{{\rm K}_1(|k|d)I^{(2)}_0(|k|)+{\rm I}_1(|k|d)K^{(2)}_0(|k|)}{{\rm K}_1(|k|d){\rm I}_0(|k|)+{\rm I}_1(|k|d){\rm K}_0(|k|)}+\eta G_{0}-\mathrm{i}\eta_xk. \end{equation}

By using the fundamental relations of the modified Bessel functions, one obtains

(2.37)\begin{equation} \left.\begin{gathered} {\rm I}_0^{(2)}(|k|)={\rm I}_1^{(1)}(|k|)={\rm I}_0(|k|)-|k|^{{-}1}{\rm I}_1(|k|),\\ {\rm K}_0^{(2)}(|k|)={-}{\rm K}_1^{(1)}(|k|)={\rm K}_0(|k|)+|k|^{{-}1}{\rm K}_1(|k|). \end{gathered}\right\} \end{equation}

Therefore, (2.36) can be simplified to

(2.38)\begin{equation} G_{1}(\eta)={-}G_{0}\eta G_{0}-\partial_x\eta\partial_x, \end{equation}

which is of the same form as that in the classical two-dimensional water-wave problem. For $j\geqslant 2$, $G_{j}$ can be obtained by a recursive formula as follows:

(2.39)\begin{align} G_{j}(\eta)&={-}\sum_{n=0}^{j-1}G_{n}(\eta)\frac{(|k|\eta)^{j-n}}{(j-n)!}M^{(j-n)}+ \frac{(|k|\eta)^j|k|}{j!}M^{(j+1)}+\frac{(|k|\eta)^j}{(j-1)!}M^{(j)}\nonumber\\ &\quad -\frac{(\mathrm{i}\eta_x)(|k|\eta)^{j-1}k}{(j-1)!}M^{(j-1)}-\frac{(\mathrm{i}\eta\eta_x) (|k|\eta)^{j-2}k}{(j-2)!}M^{(j-2)}, \end{align}

where $M^{(j)}$ is defined as

(2.40)\begin{equation} M^{(j)}=\frac{{\rm K}_1(|k|d){\rm I}_0^{(j)}(|k|)+{\rm I}_1(|k|d){\rm K}_0^{(j)}(|k|)}{{\rm K}_1(|k|d){\rm I}_0(|k|)+{\rm I}_1(|k|d){\rm K}_0(|k|)}, \end{equation}

and the dependence of $M^{(j)}$ on $|k|$ and $d$ has been suppressed for simplicity. Direct corollaries include

(2.41a–c)\begin{equation} M^{(0)}=1,\quad M^{(1)}=|k|^{{-}1}G_0,\quad M^{(2)}=1-|k|^{{-}1}M^{(1)}=1-|k|^{{-}2}G_0. \end{equation}

After rearrangement, (2.39) can be rewritten, in terms of $P$, as

(2.42)\begin{align} G_{j}(\eta)&=\frac{1}{j!}P\eta^jP|P|^{j-1}M^{(j-1)}+\frac{\eta^j}{j!}|P|^j\left[|P|M^{(j+1)}-|P|M^{(j-1)}+jM^{(j)}\right]\nonumber\\ &\quad -\frac{\eta^{j-1}\eta_x}{(j-2)!}\partial_x|P|^{j-2}M^{(j-2)}-\sum_{n=0}^{j-1}\frac{1}{(j-n)!}G_{n}(\eta)\eta^{j-n}|P|^{j-n}M^{(j-n)}. \end{align}

For the convenience of later discussion, we show the expression of $G_2(\eta )$ as follows. Since

(2.43)\begin{equation} \left.\begin{gathered} {\rm I}_1^{(2)}(|k|)=\left(1+2|k|^{{-}2}\right){\rm I}_1(|k|)-|k|^{{-}1}{\rm I}_0(|k|),\\ {\rm K}_1^{(2)}(|k|)=\left(1+2|k|^{{-}2}\right){\rm K}_1(|k|)+|k|^{{-}1}{\rm K}_0(|k|), \end{gathered}\right\} \end{equation}

one immediately obtains

(2.44)\begin{equation} M^{(3)}=\frac{{\rm K}_1(|k|d){\rm I}_1^{(2)}(|k|)-{\rm I}_1(|k|d){\rm K}_1^{(2)}(|k|)}{{\rm K}_1(|k|d){\rm I}_0(|k|)+{\rm I}_1(|k|d){\rm K}_0(|k|)}=\left(1+2|k|^{{-}2}\right)M^{(1)}-|k|^{{-}1}. \end{equation}

A direct calculation yields

(2.45)\begin{equation} \frac{\eta^2}{2}|k|^2\left[|k|M^{(3)}-|k|M^{(1)}+2M^{(2)}\right]=\frac{\eta^2}{2}|k|^2, \end{equation}

and then it follows that

(2.46)\begin{equation} G_{2}(\eta)=\tfrac12G_{0}\eta^2\partial_{xx}+\tfrac12\partial_{xx}\eta^2G_{0}+G_{0}\eta G_{0}\eta G_{0}+\tfrac12G_{0}\eta^2G_{0}-\tfrac12\partial_x\eta^2\partial_x.\end{equation}

The expression of the pseudo-differential operator $G_2(\eta )$ is formally different from that of Craig & Sulem (1993), and the difference lies in the last two items on the right-hand side of (2.46), which are not present in the classical water-wave problem. Nevertheless, $G_1$ and $G_2$ have more symmetrical structures compared with those of Guyenne & Părău (2016). For higher-order terms ($G_j$ for $j\geqslant 3$), the procedure continues to be consistent and straightforward, although more tedious calculations have been omitted for brevity.

Using the DNO, the Hamiltonian formulation (2.15) can be rewritten explicitly with the canonical variables as

(2.47)\begin{align} \mathcal{H}[\zeta,\xi]&=\frac12\int\xi G\xi\,{{\rm d}{\kern 0.9pt}x}+\frac{B}{2}\int\left[\frac12(2\eta+\eta^2)-\ln(1+\eta)\right]{\rm d}{\kern 0.9pt}x\nonumber\\ &\quad +\int\left[(1+\eta)\sqrt{1+\eta_x^2}-\frac12(2+2\eta+\eta^2)\right]{\rm d}{\kern 0.9pt}x. \end{align}

The system has a physical conserved quantity, mass, which can be expressed as

(2.48)\begin{equation} M=\int_{-\infty}^\infty\left(\eta^2+2\eta\right){\rm d}{\kern 0.9pt}x. \end{equation}

This conservation law can be verified as follows:

(2.49)\begin{equation} \frac{{\rm d}M}{{\rm d}t}=\int_{-\infty}^{\infty}\partial_t\left(\eta^2+2\eta\right){\rm d}{\kern 0.9pt}x=2\int_{-\infty}^{\infty} G(\eta)\xi\,{{\rm d}{\kern 0.9pt}x}=2\int_{-\infty}^{\infty} \xi G(\eta)1\,{{\rm d}{\kern 0.9pt}x}=0, \end{equation}

where the definitions of the DNO and its self-adjoint property have been used during the derivation. The conserved quantities, mass (2.48) and energy (2.47), can be further used to monitor the global accuracy of numerical computations.

3. Asymptotic models

In this section, a systematic procedure for deriving various nonlinear dispersive models is presented by taking advantage of the Hamiltonian structure of the problem and the analyticity property of the DNO. First of all, it is noticed that Hamilton's equations (2.24a,b) realize an extremum of the action $\int \mathbb {L}\,\textrm {d}t$ with the Lagrangian,

(3.1)\begin{align} \mathbb{L}&=\int\left\{\xi\zeta_t-\frac12\xi G\xi - \frac{B}{2}\left[\frac12\left(2\eta+\eta^2\right)-\ln(1+\eta)\right]\right\}{{\rm d}{\kern 0.9pt}x}\nonumber\\ &\quad -\int\left[(1+\eta)\sqrt{1+\eta_x^2}-\frac12\left(2+2\eta+\eta^2\right)\right]{\rm d}{\kern 0.9pt}x. \end{align}

We can construct reduced models by applying the variational principle to an approximated Lagrangian, called the Lagrangian perturbation method. This can be achieved by truncating either the Taylor series expansion of the DNO due to its analytical property or the Lagrangian expansion about a small parameter induced by a multi-scale analysis. We should note that the Lagrangian perturbation method differs slightly from the Hamiltonian perturbation method detailed in, e.g. Craig, Guyenne & Kalisch (2005); as will be shown below, it directly accommodates the use of non-canonical variables when taking variational derivatives.

3.1. Long-wave modelling

To showcase the operation routine of the suggested approach, we first consider the classical situation, addressing the problem with $d\neq 0$ within the Boussinesq regime. The Boussinesq scaling indicates that the free-surface wave varies on temporal and horizontal spatial scales of $1/\epsilon$, and the amplitude of the surface displacement scales by $\epsilon ^2$, where $\epsilon \ll 1$ is the ratio of the mean depth of fluid to typical wavelength, a small parameter characterizing the long-wave approximation. The magnetic Bond number is assumed to be of order $\epsilon ^0$ to enable the force exerted by the magnetic field to compete with surface tension. In summary, the following scales are chosen:

(3.2a–d)\begin{equation} t=\frac{\bar{t}}{\epsilon},\quad x=\frac{\bar{x}}{\epsilon},\quad \eta=\epsilon^2\bar{\eta},\quad \xi=\epsilon\bar{\xi}. \end{equation}

After changing variables and dropping bars, the pseudo-differential operator $G$ can be expanded as

(3.3)\begin{equation} G={-}\frac{\epsilon^2}{2}\left(1-d^2\right)\partial_x^2+\epsilon^4\left[A(d)\partial_x^4 -\partial_x\eta\partial_x\right]+O(\epsilon^6). \end{equation}

Upon substituting the asymptotic expansion of the DNO into the Lagrangian (3.1) and retaining terms valid up to $O(\epsilon ^5)$, one then obtains

(3.4)\begin{align} \tilde{\mathbb{L}}&=\epsilon^3\int\left[\xi\eta_t+\dfrac{1-d^2}{4}\xi\partial_x^2\xi -\dfrac{B-1}{2}\eta^2 \right]{\rm d}{\kern 0.9pt}x\nonumber\\ &\quad+\epsilon^5\int\left[\xi\eta\eta_t - \dfrac{1}{2}A(d)\xi\partial_x^4\xi+\dfrac{1}{2}\xi\partial_x\eta\partial_x\xi + \frac{B}{3}\eta^3 -\dfrac{1}{2}\eta_x^2\right]{\rm d}{\kern 0.9pt}x. \end{align}

Minimizing the action of the approximate Lagrangian $\int \tilde {\mathbb {L}}\,\textrm {d}t$ yields two evolution equations for $\eta$ and $\xi$:

(3.5)\begin{align} \eta_t&=\frac{1}{1+\epsilon^2\eta}\left[-\frac{1-d^2}{2}\xi_{xx}+\epsilon^2\left(A(d)\xi_{xxxx}-\partial_x\eta\partial_x\xi\right)\right]\nonumber\\ &={-}\frac{1-d^2}{2}\xi_{xx}+\epsilon^2\left[\frac{1-d^2}{2}\eta\xi_{xx}+A(d)\xi_{xxxx}-\partial_x \eta\partial_x\xi\right]+O(\epsilon^4) \end{align}

and

(3.6)\begin{align} \xi_{t}&={-}\frac{1}{1+\epsilon^2\eta}\left[(B-1)\eta+\epsilon^2\left(\frac{1}{2}\xi_{x}^2-\frac{B}{2}\eta^2-\eta_{xx}\right)\right]\nonumber\\ &=(1-B)\eta+\epsilon^2\left[\left(\frac{3B}{2}-1\right)\eta^2-\frac{1}{2}\xi_{x}^2+\eta_{xx}\right]+O(\epsilon^4). \end{align}

To proceed, taking the derivative of (3.6) with respect to $t$ yields

(3.7)\begin{equation} \xi_{tt}=(1-B)\eta_t+\epsilon^2\left[\eta_{txx} -\xi_{x}\xi_{xt}+\left( 3B-2\right)\eta\eta_t\right], \end{equation}

where the terms of $O(\epsilon ^4)$ and higher have been neglected. Replacing $\eta _t$ with (3.5) and retaining terms valid up to $O(\epsilon ^2)$, one obtains a single evolution equation of $\xi$:

(3.8)\begin{align} \xi_{tt}&=\frac{(B-1)(1-d^2)}{2}\xi_{xx}+\epsilon^2\left\{-\left(\xi^2_{x}\right)_t+ \left[\frac{2B-\dfrac{3}{2}}{B-1}(1-d^2)-1\right]\xi_{t}\xi_{xx}\right.\nonumber\\ &\quad\left.-\left[ \frac{1-d^2}{2}+(B-1)A(d)\right]\xi_{xxxx}\vphantom{\left[\frac{2B-\dfrac{3}{2}}{B-1}(1-d^2)-1\right]}\right\}, \end{align}

where the relation $\eta \sim \xi _t/(1-B)$ arising from the leading order of (3.6) has been used. Furthermore, by assuming the unidirectional propagation of waves, the classic KdV equation can be obtained,

(3.9)\begin{equation} \eta_t +c_0\eta_x+\alpha\eta_{xxx}+\beta\eta\eta_x=0,\end{equation}

with the dispersive and nonlinear coefficients given by

(3.10)$$\begin{gather} \alpha={-}c_0\left[ \frac{1}{2(B-1)}+\frac{A(d)}{1-d^2} \right], \end{gather}$$

(3.11)$$\begin{gather}\beta=\frac{B-1}{2c_0} \left[ 3-\frac{1-d^2}{2}-\frac{ \left(3B-2\right)(1-d^2) }{2(B-1)}\right], \end{gather}$$

and the long wave speed $c_0$ defined in (2.13a,b). It is worth noting that we have dropped $\epsilon$ here by reverting to the original variables. Since the derivation from the bidirectional model to the unidirectional is standard, we state the resultant equation and refer interested readers to Whitham (1974) or Guan & Wang (2022) for details. As $d\rightarrow 0$, (3.9) is precisely the same as derived by Rannacher & Engel (2006). The celebrated (single) soliton solution to (3.9), corresponding to a locally confined travelling wave, is given by

(3.12)\begin{equation} \eta(x,t)=\frac{3\delta}{\beta}\mathrm{sech}^{2}\left(\sqrt{\frac{\delta}{4\alpha}}\left(x-(\delta+c_0)t\right)\right), \end{equation}

where $\delta \ll 1$ is a free constant with the same sign as $\alpha$. A critical parameter, $B_1(d)$, can be found to identify elevation solitons for $1< B< B_1(d)$ and depression solitons for $B_1(d)< B< B_2(d)$, both bifurcating from zero amplitude; this critical value is expressed by

(3.13)\begin{equation} B_1(d)=\frac{3(1+d^2)}{2(1+2d^2)}. \end{equation}

As is well known, soliton formation arises from a dynamic balance between dispersive and nonlinear effects. However, the dispersive and nonlinear coefficients in (3.9), $\alpha$ and $\beta$, may be close to zero for certain parameter sets; hence, rescaling is required in the asymptotic analysis. For the first case, namely $\beta \ll 1$, a higher-order nonlinearity must be introduced to balance the dispersion, resulting in the modified KdV equation. To proceed, we choose a set of new scales

(3.14a–e)\begin{equation} t=\frac{\bar{t}}{\epsilon},\quad x=\frac{\bar{x}}{\epsilon},\quad\eta=\epsilon\bar{\eta},\quad \xi=\bar{\xi},\quad \beta=\epsilon\bar{\beta}.\end{equation}

Following the same procedure, symbolic calculations yield

(3.15)\begin{equation} G ={-}\epsilon^2\frac{1-d^2}{2}\partial_x^2 -\epsilon^3\partial_x\eta\partial_x+\epsilon^4A(d)\partial_x^4 +O(\epsilon^5).\end{equation}

Substituting (3.14a–e) and (3.15) into (3.1) and truncating the expansion at $O(\epsilon ^4)$, an approximate Lagrangian can be obtained as

(3.16)\begin{align} \tilde{\mathbb{L}}&=\int\left\{\epsilon\left[\xi\eta_t-\frac{1-d^2}{4}\xi_x^2-\frac{A(d)}{2}-\dfrac{B-1}{2}\eta^2\right]\right.\nonumber\\ &\quad \left.+\;\epsilon^2\left[\xi\eta\eta_t-\dfrac{1}{2}\eta\xi_x^2+\dfrac{B}{6}\eta^3\right]-\epsilon^3\left[\dfrac{B}{8}\eta^4+\dfrac{1}{2}\eta_x\right] \right\}{\rm d}{\kern 0.9pt}x. \end{align}

Minimizing the action of the approximate Lagrangian yields

(3.17)\begin{gather} \begin{aligned}[b] \eta_t &={-}\frac{1-d^2}{2}\xi_{xx}+\epsilon\left(\dfrac{1-d^2}{2}\eta\xi_{xx}-\partial_x\eta\partial_x\xi\right)\nonumber\\ &\quad +\epsilon^2\left[A(d)\xi^4_{x}+\eta\partial_x\eta\partial_x\xi-\dfrac{1-d^2}{2}\eta^2\xi_{xx}\right], \end{aligned}\end{gather}

(3.18)\begin{gather} \xi_t ={-}(B-1)\eta+\epsilon\left[\left(\tfrac{3}{2}B-1\right)\eta^2-\tfrac{1}{2}\xi_x^2\right] +\epsilon^2\left[\tfrac{1}{2}\eta\xi_x^2+\left(1-2B\right)\eta^3+\eta_{xx}\right], \end{gather}

where the variational principle is retained up to $O(\epsilon ^2)$. Following the standard argument, after some algebra, a modified KdV equation can be derived to govern the free surface dynamics

(3.19)\begin{equation} \eta_t+c_0\eta_x+\alpha\eta_{xxx}+\beta\eta\eta_x+\beta_1\eta^2\eta_x=0,\end{equation}

where $\epsilon$ has been removed by returning to the original variables, and the coefficient of the cubic nonlinearity reads

(3.20)\begin{equation} \beta_1=\frac{(B-1)^2}{2c_0}\left[\frac{6+d^2}{2(1-B)}-\frac{2(3B-2)}{(B-1)^2}+\frac{(1-d^2)(9B-5)}{2(B-1)^2}\right]. \end{equation}

Of note is the balance between the two effects, which can be achieved for larger amplitude waves indicated by the scales of (3.14a–e), although the temporal and spatial scales for solitons remain the same. It is well known that the modified KdV equation (3.19) has a one-parameter family of soliton solutions:

(3.21)\begin{equation} \eta(x,t)=\frac{6\alpha}{\beta\varXi^2+\sqrt{\beta^2\varXi^4+6\alpha\beta_1\varXi^2}\cosh\left[(x-c_0t-\alpha t/\varXi^2)/\varXi\right]}, \end{equation}

where $\varXi$ is the free parameter.

For the second case, $\alpha \ll 1$, a fifth-order derivative term needs to be introduced to balance the quadratic nonlinearity, resulting in the fifth-order KdV equation. For this purpose, we choose the scales

(3.22a–e)\begin{equation} t=\frac{\bar{t}}{\epsilon},\quad x=\frac{\bar{x}}{\epsilon},\quad\eta=\epsilon^4\bar{\eta},\quad \xi=\epsilon^3\bar{\xi},\quad \alpha=\epsilon^2\bar{\alpha}. \end{equation}

By assuming one-way wave propagation along a primary direction (right-going waves, say) and following a standard argument, one finally, after some tedious calculation, obtains a fifth-order KdV equation,

(3.23)\begin{equation} \eta_t+c_0\eta_x+\alpha\eta_{xxx}+\beta\eta\eta_x+\alpha_1\eta_{xxxxx}=0, \end{equation}

in terms of the original variables. The fifth-order dispersive coefficient reads

(3.24)\begin{equation} \alpha_1=\frac{1}{2c_0}\left[ \frac{1-d^2}{2(1-B)}+(B-1)A_1(d)\right] \end{equation}

with

(3.25)\begin{equation} A_1(d)=\tfrac{1}{192}\left[2-15d^2+30d^4-17d^6+12d^4({-}2+3d^2)\ln d-24d^6(\ln d)^2\right]. \end{equation}

The coefficients of the fifth-order dispersive term can also be obtained by directly expanding the dispersion relation (2.12) up to terms of $O(\epsilon ^4)$. It is also highlighted that Groves & Nilsson (2018) provided a rigorous derivation of the KdV and fifth-order KdV equations within this context.

In the water-wave community, in addition to weakly nonlinear dispersive models, regularized strongly nonlinear models are usually required to describe moderate-amplitude wave phenomena as well as to afford a significant simplification over primitive equations. For ferrofluid jets, it was Cornish (2018) who first proposed a strongly nonlinear and weakly dispersive model by extending the model of Papageorgiou & Orellana (1998) via including the magnetic effect. We propose in the subsequent analyses a new strongly nonlinear model using the Lagrangian perturbation method and show that the model documented by Cornish (2018) can also be obtained based on our method. First, we introduce a new set of non-dimensional scales

(3.26a–e)\begin{equation} x=l\bar{x},\quad r=R\bar{r},\quad S=R\bar{S},\quad t=\sqrt{\frac{\rho R l^2}{\sigma}}\bar{t},\quad \phi=\sqrt{\frac{\sigma l^2}{\rho R}}\bar{\phi}, \end{equation}

where $l$ is the typical wavelength and $\eta$ is assumed to be $O(1)$ (recalling that $S=1+\eta$). After changing variables and dropping bars, the Lagrangian $\mathbb {L}$ becomes

(3.27)\begin{equation} \mathbb{L}=\int\left[\xi \zeta_t-\frac12\int_d^S\left(\phi_x^2+\frac{1}{\epsilon^2}\phi_{r}^2\right) {\rm d}r+\left(\frac{B}{2}\ln S-S\sqrt{1+\epsilon^2 S_x^2}\right)\right]{\rm d}{\kern 0.9pt}x, \end{equation}

where $\epsilon =R/l$ is a small parameter reflecting the long-wave approximation. We denote by $\tilde {G}$ the rescaled DNO, involving the small parameter $\epsilon$. Since

(3.28)\begin{equation} \tilde{G}_0={-}\frac{\epsilon^2(1-d^2)}{2}\partial_x^2+\epsilon^4A(d)\partial_x^4+O(\epsilon^6), \end{equation}

the pseudo-differential operator $\tilde {G}$ can be expanded, in terms of $\epsilon$, as

(3.29)\begin{align} \tilde{G}&=\epsilon^2\left(-\frac{1-d^2}{2}\partial_x^2-\partial_x\eta\partial_x-\frac12\partial_x\eta^2\partial_x\right)\nonumber\\ &\quad+\epsilon^4 \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}

The approximate Lagrangian is constructed by inserting (3.29) into (3.1) and retaining terms valid up to $O(\epsilon ^2)$, resulting in

(3.30)\begin{align} \tilde{\mathbb{L}}&=\int\left[\xi\zeta_{t}+\frac{1}{2}\xi\left(\frac{1-d^2}{2}\partial_x^2+\partial_x\eta\partial_x+\frac{1}{2}\partial_x\eta^2\partial_x\right)\xi+\frac{B}{2}\ln(1+\eta)-(1+\eta)\right]{\rm d}{\kern 0.9pt}x\nonumber\\ &\quad-\frac{\epsilon^2}{2}\int\left\{A(d)\xi\partial_x^4\xi-\frac{(1-d^2)^2}{4}\xi\partial_{x}^2\eta\partial_x^2\xi\right.\nonumber\\ &\quad\left.+\left[\frac{(1-d^2)^2}{8}-\frac{1-d^2}{2}\right]\xi\partial_x^2\eta^2\partial_x^2\xi+(1+\eta)\eta_x^2\right\}{\rm d}{\kern 0.9pt}x. \end{align}

Minimizing the approximate action $\int \tilde {{\mathbb {L}}}\,\textrm {d}t$ yields

(3.31)\begin{align} 0&=\eta_t+u\eta_x+\frac{(1+\eta)^2-d^2}{2(1+\eta)}u_x-\frac{\epsilon^2}{1+\eta}\left[\vphantom{\frac{(1-d^2)^2}{4}}A(d)u_{xxx}\right.\nonumber\\ &\quad\left.-\frac{(1-d^2)^2}{4}\partial_x^2(\eta u_x)+\frac{d^4+2d^2-3}{8}\partial_x^2(\eta^2u_x)\right] \end{align}

and

(3.32)\begin{align} 0&=u_t+uu_x-\frac{(1+\eta-B)\eta_x}{(1+\eta)^3}-\epsilon^2\left[\frac{(1-d^2)^2(\eta-1)}{8(1+\eta)}u_x^2\right.\nonumber\\ &\quad\left.-\frac{1-d^2}{2(1+\eta)}\eta u_x^2\right]_x-\epsilon^2\left[ \frac{\eta_x^2}{2(1+\eta)}+\eta_{xx} \right]_x, \end{align}

where the governing system has been rewritten in a more revealing form by letting $u=\xi _x$. The small parameter $\epsilon$ can be removed from the problem by stretching variables: $x\rightarrow \epsilon x$, $t\rightarrow \epsilon t$ and $u\rightarrow u$.

It is worth mentioning that, although based on the same scaling, the strongly nonlinear model established by Cornish (2018), which does not have a clear Hamiltonian conservation form, is different from (3.31)–(3.32). This is not only because a different unknown variable was adopted, but also due to a certain special treatment on the kinematic boundary condition. Next, we will show that the Cornish system can also be obtained using the Lagrangian perturbation method. To proceed the derivation, we rewrite the Laplace equation with the rescaled velocity potential as

(3.33)\begin{equation} \phi_{rr}+\frac{1}{r}\phi_{r}+\epsilon^2\phi_{xx}=0,\quad\text{in}\ d< r< S,\end{equation}

subject to the impermeability boundary condition $\phi _r=0$ at $r=d$. The asymptotic expansion of $\phi$ can be written as $\phi =\phi ^{(0)}+\epsilon ^2\phi ^{(1)} +O(\epsilon ^4)$. The substitution of the above expansion into (3.33) gives the leading and next-to-leading order terms:

(3.34a,b)\begin{equation} \phi^{(0)}=\varPhi(x,t),\quad\phi^{(1)}=\varPhi_{xx}\left(-\frac{r^2}{4}+\frac{d^2}{2}\ln r\right). \end{equation}

The approximate Lagrangian is constructed by substituting (3.34a,b) into (3.27) and retaining terms up to $O(\epsilon ^2)$, resulting in

(3.35)\begin{align} \tilde{\mathbb{L}}&=\int\left[\varPhi\zeta_{t}-\frac{S^2-d^2}{4}\varPhi_x^2 +\frac{B}{2}\ln S-S\right]{\rm d}{\kern 0.9pt}x+\epsilon^2\int\varPhi_{xx}\left(\frac{d^2}{2}\ln S-\frac{S^2}{4}\right)\zeta_t\,{\rm d}{\kern 0.9pt}x\nonumber\\ &\quad-\frac{\epsilon^2}{2}\left\{\int \varPhi_x\varPhi_{xxx}\left[\frac{3d^4-S^4-2S^2d^2}{8}+\frac{1}{2}d^2\left( S^2\ln S-d^2\ln d\right)\right]{\rm d}{\kern 0.9pt}x\right.\nonumber\\ &\quad\left.+\int\left[\varPhi_{xx}^2\left(\frac{S^4+3d^4-4d^2S^2}{16}+\frac{1}{4}d^2\ln\frac{S}{d} \right)+SS_x^2\right]{\rm d}{\kern 0.9pt}x\right\}. \end{align}

Again, we obtain the governing equations by taking the variational derivatives of the approximate action $\int \tilde {{\mathbb {L}}}\,\textrm {d}t$:

(3.36)\begin{align} &S_{t}+S_{x}u+\frac{S^2-d^2}{2S}u_x+\frac{\epsilon^2}{S}\left[\left(\frac{d^2}{2}\ln S-\frac{S^2}{4}\right)SS_t\right]_{xx}\nonumber\\ &\quad +\frac{\epsilon^2}{2S}\left\{u_{xx}\left[\frac{3d^4-S^4-2d^2S^2}{8}+\frac{d^2}{2}\left(S^2\ln S-d^2\ln d\right)\right]\right\}_x\nonumber\\ &\quad +\frac{\epsilon^2}{2S}\left\{u\left[\frac{3d^4-S^4-2d^2S^2}{8}+\frac{d^2}{2}\left(S^2\ln S-d^2\ln d\right)\right]\right\}_{xxx}\nonumber\\ &\quad-\frac{\epsilon^2}{S}\left[u_x\left(\frac{S^4+3d^4-4d^2S^2}{16}+\frac{d^2}{4}\ln\frac{S}{d}\right)\right]_{xx}=0\,, \end{align}

(3.37)\begin{align} &u_t+uu_x+\left(\frac{1}{S}-\frac{B}{2S^2}\right)_x+\frac{\epsilon^2}{4}\left[\left(2d^2\ln S-S^2\right)\left( u_{xt}+uu_{xx} \right)\right]_x\nonumber\\ &\quad +\frac{\epsilon^2}{2}\left[u_x^2\left(\frac{d^2-S^2}{2S}\right)^2\right]_x-\epsilon^2\left(\frac{S_x^2}{2S}+S_{xx}\right)_x=0. \end{align}

Following the idea of Cornish (2018) to retain only the leading order terms in the kinematic boundary condition, i.e. neglect the terms of $O(\epsilon ^2)$ in (3.36), one arrives at the same governing system; see (4.50) and (4.51) of Cornish (2018).

3.2. A fully dispersive model

Hitherto, we derive various models in the long-wave approximation. Alternatively, these models only describe wave behaviours near $k=0$ and are suitable for $1< B\leqslant B_2(d)$. However, when the phase speed minimum occurs at $k\neq 0$ or, equivalently, $B>B_2(d)$, the narrow band approximation and normal form analysis are commonly conducted near this minimum, resulting in the nonlinear Schrödinger type equation (see Groves & Nilsson (2018) for the case of ferrofluid jets). Next, we expect to propose a model that can capture the essential features near the phase speed minimum for an arbitrary magnetic Bond number. Therefore, the fully dispersive model, which accurately describes linear waves across the whole spectrum of wavelengths, is a good candidate.

In this part, the main goal is to derive a cubic truncation model, which incorporates the complete linear dispersion relation of the full potential flow equations, and keeps the full nonlinearity involving the surface tension and magnetic field. A similar counterpart has been proposed for capillary-gravity waves on water of infinite depth in the classical water-wave system (the interested readers are referred to Wang & Milewski (2012) for more details). The non-dimensionalization is made in the same way as (2.7a–e). First, in the sense explained by Craig & Sulem (1993), the DNO can be expanded as a Taylor power series (see § 2.3). The first three terms, expressed by (2.34), (2.38) and (2.46), are used to approximate the DNO, namely $G\approx G_0+G_1+G_2$. In this set-up, the Lagrangian (3.1) can be approximated by

(3.38)\begin{align} \tilde{\mathbb{L}}&\triangleq\int\left[\xi(1+\eta)\eta_t-\frac{1}{2}\xi \left(G_0+G_1+G_2 \right)\xi-\frac{B}{2}\left(\frac{1}{2}(2\eta+\eta^2)-\ln(1+\eta)\right)\right]{\rm d}{\kern 0.9pt}x\nonumber\\ &\quad -\frac{1}{2}\int\left[ (1+\eta)\sqrt{1+\eta_x^2}-\frac{1}{2}\left( 2+2\eta+\eta^2 \right)\right]{\rm d}{\kern 0.9pt}x. \end{align}

Following the same procedure, by ignoring higher order terms and then taking the variational derivatives of the truncated action, $\int \tilde {\mathbb {L}}\,\textrm {d}t$, with respect to $\xi$ and $\eta$, we obtain the cubic full-dispersion model:

(3.39)\begin{gather} \eta_t=\frac{1}{1+\eta}\left(G_0+G_1+G_2\right)\xi,\end{gather}

(3.40)\begin{gather} \begin{aligned}[b] \xi_t &={-}\frac{1}{2(1+\eta)}\left[ \left(G_0\xi \right)\left({-}G_0\xi+2\eta\xi_{xx}+\eta G_0\xi+2G_0\eta G_0\xi\right)\right]\\ &\quad -\frac12\xi_x^2+\frac{B}{2}\left[\frac{1}{\left(1+\eta \right)^2}-1\right]-(\kappa-1).\end{aligned} \end{gather}

This is more than just a weakly nonlinear model. First, as Craig & Sulem (1993) pointed out, the DNO expansion in Cartesian coordinates does not assume a smallness of wave amplitude and surface potential, but instead uses its analytical property, namely the Taylor expansion of the DNO converges under certain regularity conditions on the surface displacement, differing from other high-order spectral methods for water waves. It is reasonable to assume that the analytical property of the operator also holds in the axisymmetric cylindrical configuration, though we leave the rigorous analyses outside the scope of this paper. Second, the restoring forces exerted by the magnetic field and surface tension are retained precisely.

Following the same procedure, it is natural to derive other truncation models, such as the quadratic, quartic and quintic truncation models, corresponding to the third-, fifth- and sixth-order kinetic energy truncations. For conciseness, we omit the detailed derivations and governing equations for other truncation models. We finally remark that a limitation of the cubic full-dispersion model is that overturning waves with a multi-valued profile are excluded. In the subsequent sections, we will conduct theoretical and numerical research for the newly proposed model (3.39)–(3.40). We will demonstrate the effectiveness of this fully dispersive system from the perspective of solitary waves, including the bifurcation mechanism, stability and dynamics, as well as define the model's applicability.

4. Solitary waves

4.1. Numerical method

Despite the dimension reduction compared with the full Euler equations, finding theoretical solutions in the newly developed model is still overly complicated. To examine the qualitative similarities between the cubic full-dispersion model and the full Euler equations, and to assess the quantitative validity of the reduced model, we seek help from numerical tools. The pseudo-spectral method is applied to numerically solve the nonlinear dispersive equations (3.39)–(3.40), where the integrating factors technique is incorporated to handle the stiffness problem. For convenience, we recast the system to a single evolution equation. Numerical experiments are implemented on a double-periodic domain, where the Fourier transform is applicable in spatial variables, and hence

(4.1)\begin{equation} \begin{pmatrix} \hat{\eta}\\ \hat{\xi} \end{pmatrix}_t+\begin{bmatrix} 0 & -\hat{G}_0\\ k^2-1+B & 0 \end{bmatrix} \begin{pmatrix} \hat{\eta}\\ \hat{\xi} \end{pmatrix}=\begin{pmatrix} \widehat{\mathcal{N}_1}\\ \widehat{\mathcal{N}_2} \end{pmatrix},\end{equation}

where the hat indicates the Fourier transform and $k$ is the wavenumber in the Fourier space. The nonlinear terms $\mathcal {N}_1$ and $\mathcal {N}_2$ read

(4.2)\begin{equation} \left.\begin{gathered} \mathcal{N}_1=\frac{1}{1+\eta}\left(G_1\xi+G_2\xi\right),\\ \mathcal{N}_2={-}\frac{1}{2(1+\eta)}\left[\left(G_0\xi\right)\left({-}G_0\xi+2\eta\xi_{xx}+\eta G_0\xi+2G_0\eta G_0\xi\right)\right]-\frac12\xi_x^2\\ \quad\quad+\frac{B}{2}\frac{3\eta^2+2\eta^3}{(1+\eta)^2}-\frac{\eta^2\left(1+\eta_x^2+\sqrt{1+\eta_x^2}\right)-\eta_x^2}{(1+\eta)\left(1+\eta_x^2+\sqrt{1+\eta_x^2}\right)}-\left[\frac{\eta_x^3}{1+\eta_x^2+\sqrt{1+\eta_x^2}}\right]_x. \end{gathered}\right\} \end{equation}

We then multiply (4.1) by a matrix to diagonalize the system:

(4.3)\begin{equation} \begin{bmatrix} 1 & \dfrac{\mathrm{i}\hat{G}_0}{w}\\[10pt] 1 & -\dfrac{\mathrm{i}\hat{G}_0}{w} \end{bmatrix} \begin{pmatrix} \hat{\eta}\\ \hat{\xi} \end{pmatrix}_t+ \begin{bmatrix} \mathrm{i}w & 0 \\ 0 & -\mathrm{i}w \end{bmatrix} \begin{bmatrix} 1 & \dfrac{\mathrm{i}\hat{G}_0}{w}\\[10pt] 1 & -\dfrac{\mathrm{i}\hat{G}_0}{w} \end{bmatrix} \begin{pmatrix} \hat{\eta}\\ \hat{\xi} \end{pmatrix}= \begin{bmatrix} 1 & \dfrac{\mathrm{i}\hat{G}_0}{w} \\[10pt] 1 & -\dfrac{\mathrm{i}\hat{G}_0}{w} \end{bmatrix} \begin{pmatrix} \hat{\mathcal{N}}_1\\ \hat{\mathcal{N}}_2 \end{pmatrix}. \end{equation}

Introducing notation $\hat {p}=\hat {\eta }+({\mathrm {i}\hat {G}_0}/{w})\hat {\xi }$ and $\hat {q}=\hat {\eta }-({\mathrm {i}\hat {G}_0}/{w})\hat {\xi }$, we can rewrite the system as

(4.4)$$\begin{gather} \hat{p}_t+\mathrm{i}w\hat{p}=\hat{\mathcal{N}}_1+\frac{\mathrm{i}\hat{G}_0}{w}\hat{\mathcal{N}}_2, \end{gather}$$

(4.5)$$\begin{gather}\hat{q}_t-\mathrm{i}w\hat{q}=\hat{\mathcal{N}}_1-\frac{\mathrm{i}\hat{G}_0}{w}\hat{\mathcal{N}}_2. \end{gather}$$

Using the fact that $\xi$ and $\eta$ are both real, these two equations are essentially equivalent, and $\xi$ and $\eta$ can be recovered from $p$ alone with

(4.6)$$\begin{gather} \hat{\eta}(k)=\frac{1}{2}\left[\hat{p}(k)+\hat{p}({-}k)^* \right], \end{gather}$$

(4.7)$$\begin{gather}\hat{\xi}(k)=\frac{1}{2\mathrm{i}}\frac{w}{\hat{G}_0}\left[\hat{p}(k)-\hat{p}({-}k)^*\right], \end{gather}$$

where the asterisks represent complex conjugation. Thus, the problem is reduced to solving (4.4), a single complex evolution equation. It is worth mentioning that $\hat {\xi }(0,t)$ requires special treatment since the right-hand side of (4.7) becomes singular for $k=0$. This can be generally overcome by directly integrating the second equation of (4.1) for $k=0$. Nevertheless, upon noticing that on the right-hand sides of (3.39)–(3.40), $\xi$ always appears in the form of spatial derivatives or being acted on by the pseudo-differential operator $G_0$, $\hat {\xi }(0)$ has no substantial contribution to other Fourier coefficients and hence there is no need to update at every time step.

To compute a travelling wave translating with a constant speed $c$, we set $\widehat {p_t}=-\mathrm {i}ck\hat {p}$ and solve the resulting nonlinear algebraic system for the Fourier coefficients using the Newton–Raphson method. The initial guess for the Newton algorithm is obtained by using the vanishing pressure method, which was first presented by Vanden-Broeck & Dias (1992). Time integration of the system is accomplished with a classic fourth-order Runge–Kutta method combined with an integrating factor. All the computations are de-aliased with a doubling of Fourier modes, and no filters are used. For solitary waves of wavepacket type, small-amplitude solutions often present significant challenges in terms of accurate computation due to their slower spatial decay. Consequently, it is necessary to gradually expand the computational domain as the amplitude becomes smaller (or equivalently, the solutions are in the vicinity of the bifurcation point).

4.2. Bifurcation of solitary waves and model accuracy

The numerical investigation of axisymmetric solitary-wave solutions in the fully nonlinear regime for arbitrary values of $d$ was performed by Blyth & Părău (2014) for the one-layer scenario and Doak & Vanden-Broeck (2019) for the two-layer counterpart via using a finite difference method. They solved the full Euler equations using the hodograph transformation that exchanges the dependent and independent variables. Based on the fact that the Stokes stream function $\psi$ on the free surface and the rigid bottom is constant for a solitary wave in the co-moving frame of reference, the hodograph transformation transforms the computational domain into a rectangular stripe in the $(\phi,\psi )$ plane. Then, the rewritten equations and boundary conditions for $x(\phi,\psi )$ and $r(\phi,\psi )$ can be discretized over a rectangular grid and solved using the finite difference method. We omit the complete description of the numerical scheme, and the interested readers are referred to Blyth & Părău (2014) and Doak & Vanden-Broeck (2019) for more details.

In this subsection, we focus on justifying that the cubic full-dispersion model is in excellent quantitative agreement with the full Euler equations across a range of physical parameter scenarios. For both systems, solitary waves were computed by Newton's method in the appropriate ranges of $B$ and $d$, and the bifurcation branches were found by straightforward numerical continuation in a chosen parameter (speed or amplitude). A solution was considered to have converged when the $l^\infty$-norm of the residual error was less than $10^{-10}$. It should be pointed out that for the Euler equations, the computational domain in the $\phi$ direction, originally unbounded, was truncated, and an extra condition $r_\phi =0$ was imposed at the endpoints to numerically approximate solitary waves ascribed to their decaying nature in the far field. However, solitary waves in the cubic full-dispersion model were approximated by long periodic waves with flat tails, and the period and mesh point number were both chosen sufficiently large so that, to the numerical accuracy we used, the solution did not change when they were further enlarged.

It is well known that monotonic elevation and depression solitary waves were identified in the regimes of $1< B< B_1$ and $B_1< B< B_2$, respectively, both bifurcating from the uniform stream at the global minimum of the phase speed (denoted by $c_0$), which can be understood from the KdV theory since these minima occur at $k=0$. We should clarify that only symmetric solitary waves are considered, and half of the solution is ‘monotonic’. We set $d=5/11$, which is the value of $d$ used in the experiments of Bourdin et al. (2010), so that $B_1\approx 1.28$ and $B_2\approx 15.55$. For demonstration purposes, two typical magnetic Bond numbers, $B=1.25$ and $B=1.75$, were chosen so that the respective elevation and depression solitary waves predicted by the associated KdV equation can be expected in the model equations. Comparisons of the speed-amplitude bifurcation diagram of the cubic full-dispersion model to those of the Euler equations and the theoretical predictions of the KdV equation are plotted in figure 2. The quantitative agreement between the newly proposed model and the full Euler equations is reasonably good, extending far beyond the narrow validity range of the KdV region. Typical profiles shown in figure 3 illustrate that the cubic full-dispersion model is remarkably accurate at relatively large amplitudes; specifically, the relative difference in profiles between the Euler equations and the model is approximately of order $10^{-4}$ on this scale.

Figure 2. Speed-amplitude bifurcation curves of axisymmetric ferrofluid solitary waves obtained with $d=5/11$ and two magnetic Bond numbers: $B=1.25$ (on the left-hand side) and $B=1.75$ (on the right-hand side), emerging from bifurcation points $c_0=0.3149,0.5455$, respectively. The solutions were obtained with different models: the full Euler equations (black dashed lines), the cubic full-dispersion model (blue solid lines) and the KdV equation (red dash-dotted lines). The inset shows the bifurcation curve of the elevation branch for $B=1.75$ computed with the full Euler equations.

Figure 3. ($a$) Typical wave profiles for $B=1.25$ and $S(0)=0.6$, $1.10$, $1.50$, $1.90$. ($b$) Typical wave profiles for the depression branch for $B=1.75$ and $S(0)=0.7$. Solutions were obtained with the cubic full-dispersion model (blue solid lines), the full Euler equations (black dashed lines) and the KdV theory (red dash-dotted line).

A striking phenomenon found in ferrofluid jets is the existence of another type of axisymmetric solitary waves of the opposite polarity in the Euler equations for a fixed magnetic Bond number satisfying $1< B<2$, as shown by Blyth & Părău (2014). The weakly nonlinear theory cannot capture these new solutions since they do not bifurcate from zero amplitude. For $B=1.25$, in contrast to elevation solitary waves predicted by the KdV equation, depression waves were also found in both the cubic full-dispersion model and the Euler equations, and excellent agreement was obtained for the bifurcation diagram (figure 2) and wave profile (figure 3a), confirming the effectiveness of the new model. However, it should be noted that the cubic full-dispersion model is impractical to capture the branch bifurcating from a considerably large finite amplitude. This limitation is demonstrated by the model's failure to converge to a solution at a magnetic Bond number of $B=1.75$. However, at this specific value, the numerical results from the full Euler equations reveal that an elevation branch actually bifurcates from a finite amplitude of $S(0)\approx 2.98$ at $c=c_0$, as illustrated in the embedded diagram in figure 2. Taking this limitation into account, we keep $d=5/11$ unchanged while selecting a magnetic field of $B=1.30$ in the vicinity of the critical parameter $B_1$. This choice is based on the fact that, for a fixed $d$, the solution amplitude at the bifurcation point decreases as $B$ approaches $B_1$. As shown in figure 4, the finite-amplitude elevation branch bifurcating with $S(0)\approx 1.1$ at $c=c_0$ is identified by using the cubic full-dispersion model. As the parameter $B$ continues to increase – though not depicted in a figure – the divergence between the speed-amplitude bifurcation curve of the finite-amplitude elevation branch derived from the cubic model and that of the full Euler equations becomes increasingly pronounced. This growing discrepancy ultimately leads to the failure of the cubic model to converge to a solution.

Figure 4. ($a$) Speed-amplitude bifurcation diagrams for $B=1.30>B_1$ and $d=5/11$ obtained with the cubic full-dispersion model (blue solid line), the Euler equations (black dashed line) and the KdV theory (red dash-dotted line). ($b$) Typical wave profiles for $S(0)=0.6,1.6,1.9$, and the line styles are the same as in panel ($a$).

Given that the nonlinear coefficient $\beta$ in (3.9) is close to zero, achievable by selecting $B\approx 3(1+d^3)/(2+4d^2)$, and thus cannot balance the dispersion, a rescaling is invoked, resulting in the modified KdV equation (3.19), which can be used to describe the small- and moderate-amplitude solutions. The following four sets of parameters are chosen: $d=0.5,0.4,0.3,0.2$, and the associated magnetic Bond numbers $B=1.26,1.33,1.40,1.48$. To highlight the superiority of the cubic full-dispersion model in the fully nonlinear regime, solitary waves were also computed in the strongly nonlinear model (3.31)–(3.32) and participated in comparison. As a locally confined travelling wave translating with speed $c$, the system can be integrated once, yielding

(4.8)\begin{align} {-}\frac{c}{2}\!\left(\eta^2+2\eta\right)+\frac{(1+\eta)^2-d^2}{2}u =\epsilon^2\left[\!A(d)u_x-\frac{(1-d^2)^2}{4}\eta u_x+\frac{d^4+2d^2-3}{8}\eta^2 u_x\!\right]_x \end{align}

and

(4.9)\begin{align} &-cu +\frac{u^2}{2}-\frac{\eta}{1+\eta}+\frac{B}{2}\frac{\eta^2+2\eta}{(1+\eta)^2}\nonumber\\ &\quad=\epsilon^2\left[\frac{(1-d^2)^2(\eta-1)}{8(1+\eta)}u_x^2-\frac{1-d^2}{2(1+\eta)}\eta u_x^2+\frac{\eta_x^2}{2(1+\eta)}+\eta_{xx}\right], \end{align}

which is then solved by the Newton iteration method. Depression solitary waves can be found for all four parameter sets we tested, and the speed-amplitude curves from four types of models are plotted in figure 5(a). The modified KdV equation offers qualitative descriptions for the small- and moderate-amplitude solitary waves, and the cubic full-dispersion model and the strongly nonlinear model show overall consistency with the magnetized Euler equations. Typical profiles for depression waves with $S(0)=0.8$ and $S(0)=0.6$ are shown in figures 5(b) and 5(c) for $d=0.3$ and $B=1.40$. The cubic full-dispersion model outperforms other reduced models when compared with the full Euler equations, as the two wave profiles overlap significantly, making the difference nearly indistinguishable. Additionally, the situation of $\alpha \ll 1$ is also examined. In figure 6(a), we show the dispersion relations of the fifth-order KdV equation and the full Euler equations with $d=0.3$ and $B=10.2$. The solitary-wave solutions to (3.23) were obtained through an extension of the Petviashvili method, as described by Guan & Wang (2022), with the speed-amplitude bifurcation curves and typical profiles shown in figure 6(b). As expected, the cubic full-dispersion model still performs excellently. However, as is well known, unlike the modified KdV equation, the fifth-order KdV equation is only valid in the small-amplitude limit considering the scale balance in the asymptotic analysis.

Figure 5. ($a$) From left to right, it shows the speed-amplitude bifurcation curves for $d=0.5,0.4,0.3,0.2$ (the corresponding values of $B$ are $1.26$, $1.33$, $1.40$ and $1.48$) obtained with different models: the cubic full-dispersion model (blue solid lines), the full Euler equations (downward-pointing black triangles), the modified KdV equation (dash-dotted lines) and the strongly nonlinear model (red asterisks). ($b$) Comparison of typical wave profiles for $d=0.3$ and $S(0)=0.8$: the cubic full-dispersion model (blue solid line), the Euler equations (black dashed line), the modified KdV equation (black dash-dotted line) and the strongly nonlinear model (red dotted line). ($c$) Comparison of typical wave profiles for $d=0.3$ and $S(0)=0.6$, and the line styles are the same as in panel ($b$).

Figure 6. ($a$) Comparison of the linear dispersion relation between the fifth-order KdV equation (red solid line) and the Euler equations (black dashed line) for $d=0.3$ and $B=10.2$. ($b$) Speed-amplitude bifurcation curves obtained from the fifth-order KdV equation (red solid line), the cubic full-dispersion model (blue solid line) and the Euler equations (black dashed line). The inset shows the typical profiles for $S(0)=0.9$, and the line styles are the same as the bifurcation curves.

As $B>B_2(d)$, the linear dispersive relation of the full Euler equations has a global minimum at a non-zero wavenumber, from which wavepacket solitary waves can bifurcate. A typical example with $d=5/11$ and $B=28$, including the bifurcation diagrams and typical profiles, is shown in figure 7. It is worth mentioning that the solutions within this region are highly oscillatory and very broad near the bifurcation point, requiring more grid points than monotonic solitary waves for $B< B_2(d)$ to resolve solutions satisfactorily, for both the Euler equations and reduced models. To further validate the rationality of the cubic truncation, we made a comparison between solutions found using the different truncations mentioned in § 3.2 (quadratic, cubic, quartic and quintic) and those of the full potential flow equations for depression wavepacket solitary waves. As shown in figure 7 for the speed-amplitude bifurcation diagrams, it is clear that, over an extensive range of wave speeds corresponding to a large range of amplitudes, the cubic full-dispersion model exhibits better performance than the other truncated models, confirming its accuracy. Regarding the bifurcation diagrams of the cubic model, we found two branches of symmetric wavepacket solitary waves, elevation (with a positive free-surface displacement at the centre) and depression (with a negative free-surface displacement at the centre), tending to meet at the point $S(0)=1$ and $c=3.021$ (which is precisely equal to the phase speed minimum). This numerical evidence confirms that the two branches simultaneously bifurcate from the infinitesimal periodic waves at this minimum point (denoted by $c_m$ hereafter). These wavepacket solitary waves are governed at small amplitude by a nonlinear Schrödinger equation, recently derived by Groves & Nilsson (2018) under the assumption of $d=0$. As for the small deviation in bifurcation curves between the cubic full-dispersion model and the Euler equations, in addition to the truncation error of the DNO, it is also attributed to the under-resolved potential flow results. This is especially noticeable near the bifurcation point $c_m$, where the waves are of small amplitude, highly oscillatory and spatially extended. There is a numerical clue that increasing the number of grid points in the $\psi$-direction from 50 to 150 makes the solution branch of the Euler equations closer to that of the cubic full-dispersion model. However, the calculation of the Euler equations is actually two-dimensional, and thus limited by computational power, so further increasing the number of grid points becomes challenging.

Figure 7. A comparison of speed-amplitude bifurcation diagram between the full Euler equations (black dashed line) and various truncation models: the quadratic model (yellow solid line), the cubic full-dispersion model (blue solid lines), the quartic model (green solid line) and the quintic model (red solid line). The branches of elevation and depression wavepacket solitary waves were computed with $d=5/11$ and $B=28$, emerging from the bifurcation point $c_m=3.021$.

In this part, all numerical evidence on solitary waves basically justifies the assertion that the cubic full-dispersion model is an appropriate reduced model to approximate the primitive Euler equations. Since the unsteady simulation of the full Euler equations is still challenging, in what follows, we use (3.39)–(3.40) as the time-dependent model to conduct a series of numerical experiments on wave dynamics, including the stability of solitary waves and collisions between two stable ones.

5. Stability and dynamics of solitary waves

5.1. Stability of solitary waves – theory

In this subsection, using the cubic full-dispersion model, we aim to establish a criterion for the stability exchange of axisymmetric solitary waves when subjected to longitudinal perturbations. To achieve this, we employ the asymptotic expansion method near the neutral mode, following a similar approach to that of Kataoka (2006).

To start, we introduce a steady solution to the model (3.39)–(3.40) in the frame of reference moving with a constant speed $c$ along the longitudinal direction $x$, which we denote as $\bar {\eta }(x-ct)$ and $\bar {\xi }(x-ct)$. We refer to $(\bar {\eta },\bar {\xi })$ as a solitary-wave solution. The existence of this solution has been confirmed through numerical computations for various sets of parameters $d$ and $B$, as demonstrated in § 4.2. Thus, $\bar {\eta }$ and $\bar {\xi }$ satisfy the following equations:

(5.1)\begin{align} -c\bar{\eta}_x&=\left[G_0\bar{\xi}-\partial_x\bar{\eta}\partial_x\bar{\xi}-G_0(\bar{\eta}G_0\bar{\xi})+\frac12G_0(\bar{\eta}^2\bar{\xi}_{xx})+ \frac12\partial_{xx}(\bar{\eta}^2G_0\bar{\xi})\right.\nonumber\\ &\quad\left.+\;\frac12G_0(\bar{\eta}^2G_0\bar{\xi})+G_0(\bar{\eta}G_0\bar{\eta}G_0\bar{\xi})-\frac{1}{2} \partial_x(\bar{\eta}^2\bar{\xi}_x)\right]\frac{1}{1+\bar{\eta}}, \end{align}

(5.2)\begin{align} -c\bar{\xi}_x&={-}\frac{1}{2(1+\bar{\eta})}\left[(G_0\bar{\xi})\left({-}G_0\bar{\xi}+2\bar{\eta}\bar{\xi}_{xx}+\bar{\eta}G_0\bar{\xi}+2G_0\bar{\eta}G_0\bar{\xi}\right)\right]-\frac12\bar{\xi}_x^2\nonumber\\ &\quad+\frac{B}{2}\left[\frac{1}{(1+\bar{\eta})^2}-1\right]-\left[\frac{1}{(1+\bar{\eta})(1+\bar{\eta}_x^2)}-\frac{\bar{\eta}_{xx}}{\left(1+\bar{\eta}_x^2\right)^{3/2}}-1\right]. \end{align}

Consider perturbations of the form

(5.3a,b)\begin{equation} \eta=\bar{\eta}+\hat{\eta}{\rm e}^{\lambda t},\quad \xi=\bar{\xi}+\hat{\xi}{\rm e}^{\lambda t},\end{equation}

where the magnitudes of $\hat {\eta }$ and $\hat {\xi }$ are much smaller than the undisturbed solitary-wave solution, and $\lambda$ is an eigenvalue. Substituting (5.3a,b) into (3.39)–(3.40) and linearizing around the fundamental state $(\bar {\eta },\bar {\xi })$, we obtain the following set of linear equations:

(5.4)\begin{equation} \left.\begin{gathered} -\lambda\hat{\eta}=\mathcal{L}_1\left[\hat{\eta},\hat{\xi}\right]=\mathcal{L}_{11}[\hat{\eta}]+\mathcal{L}_{12}[\hat{\xi}]\,,\\ -\lambda\hat{\xi}=\mathcal{L}_2\left[\hat{\eta},\hat{\xi}\right]=\mathcal{L}_{21}[\hat{\eta}]+\mathcal{L}_{22}[\hat{\xi}]\,, \end{gathered}\right\}\end{equation}

where the linearized operators are defined as follows:

(5.5)\begin{align} \mathcal{L}_{11}[\hat{\eta}]&={-}c\hat{\eta}+\frac{\hat{\eta}}{\left(1+\bar{\eta}\right)^2} \left[G_0\xi-G_0\bar{\eta}G_0\bar{\xi}-\partial_x\bar{\eta}\partial_x\bar{\xi}+\frac12G_0\left(\bar{\eta}^2\bar{\xi}_{xx} \right)\right.\nonumber\\ &\quad\left.+\;\frac12\partial_{xx}\left(\bar{\eta}^2G_0\bar{\xi} \right)+\frac12G_0\bar{\eta}^2G_0\bar{\xi} +G_0\bar{\eta}G_0\bar{\eta}G_0\bar{\xi}-\frac12\partial_x(\bar{\eta}^2\bar{\xi}_x)\right]\nonumber\\ &\quad+\frac{1}{1+\bar{\eta}}\left[G_0\hat{\eta}G_0\bar{\xi}+\partial_x(\bar{\eta}\hat{\eta}\bar{\xi}_x) -G_0(\bar{\eta}\hat{\eta}\bar{\xi}_{xx})-\partial_{xx}(\bar{\eta}\hat{\eta}G_0\bar{\xi})\right.\nonumber\\ &\quad\left.-\;G_0(\bar{\eta}\hat{\eta}G_0\bar{\xi})-G_0\hat{\eta}G_0\bar{\eta}G_0\bar{\xi}- G_0\bar{\eta}G_0\hat{\eta}G_0\bar{\xi}+\partial_x\hat{\eta}\partial_x{\bar{\xi}}\right], \end{align}

(5.6)\begin{align} \mathcal{L}_{12}[\hat{\xi}]&=\frac{1}{1+\bar{\eta}}\left[{-}G_0\hat{\xi}+G_0\bar{\eta}G_0\hat{\xi}+\partial_x\bar{\eta}\partial_x\hat{\xi}-\frac{1}{2}G_0(\bar{\eta}^2\hat{\xi}_{xx})\right.\nonumber\\ &\quad\left.-\;\frac12\partial_{xx}(\bar{\eta}^2G_0\hat{\xi})-\frac12G_0(\bar{\eta}^2G_0\hat{\xi}) -G_0\bar{\eta}G_0\bar{\eta}G_0\hat{\xi}+\frac{1}{2} \partial_x(\bar{\eta}^2\hat{\xi}_x)\right], \end{align}

(5.7)\begin{align} \mathcal{L}_{21}[\hat{\eta}]&=\frac{\hat{\eta}}{(1+\bar{\eta})^2}\left[\frac{1}{2}(G_0\bar{\xi})^2(1-\bar{\eta})-(G_0\bar{\xi})(\bar{\eta}\bar{\xi}_{xx}+G_0\bar{\eta}G_0\bar{\xi})\right]\nonumber\\ &\quad +\frac{1}{1+\bar{\eta}}\left[(G_0\bar{\xi})\left( \hat{\eta}\bar{\xi}_{xx}+\frac12\hat{\eta}G_0\bar{\xi}+G_0\hat{\eta}G_0\bar{\xi}\right)\right]+B_0\frac{\hat{\eta}}{(1+\bar{\eta})^3}\nonumber\\ &\quad-\left[\frac{\hat{\eta}}{(1+\bar{\eta})^2(1+\bar{\eta}_x^2)^{1/2}}+\frac{\bar{\eta}_x\hat{\eta}_x}{(1+\bar{\eta})(1+\bar{\eta}_x^2)^{3/2}}\right]-\left[\frac{\hat{\eta}_x}{(1+\bar{\eta}_x^2)^{3/2}}\right]_x, \end{align}

(5.8)\begin{align} \mathcal{L}_{22}[\hat{\xi}]&={-}c\hat{\xi}_x+\frac{1}{1+\bar{\eta}}\left[(G_0\hat{\xi})\left({-}G_0\bar{\xi}+\bar{\eta}\bar{\xi}_{xx}+\bar{\eta}G_0\bar{\xi}+G_0\bar{\eta}G_0\bar{\xi} \right)\right]+\bar{\xi}_x\hat{\xi}_x \nonumber\\ &\quad+\frac{1}{1+\bar{\eta}}\left[(G_0\bar{\xi})\left(\bar{\eta}\hat{\xi}_{xx}+G_0\bar{\eta}G_0\hat{\xi}\right) \right]. \end{align}

Under the assumption of $|\lambda |\ll 1$, we seek the solution $(\hat {\eta },\hat {\xi })$ in the following power series of $\lambda$:

(5.9a,b)\begin{equation} \hat{\eta}=\hat{\eta}^{(0)}+\lambda\hat{\eta}^{(1)}+\lambda^2\hat{\eta}^{(2)}+\cdots,\quad \hat{\xi}=\hat{\xi}^{(0)}+\lambda\hat{\xi}^{(1)}+\lambda^2\hat{\xi}^{(2)}+\cdots. \end{equation}

Substituting the series (5.9a,b) into (5.4) and collecting powers of $\lambda$, we obtain a series of linear problems with the eigenvalue removed for $(\hat {\eta }^{(n)},\hat {\xi }^{(n)})$:

(5.10a,b)\begin{equation} \mathcal{L}_1\left[\hat{\eta}^{(n)},\hat{\xi}^{(n)}\right]=K^{(n)}, \quad\mathcal{L}_2\left[\hat{\eta}^{(n)},\hat{\xi}^{(n)}\right]=H^{(n)}, \end{equation}

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

(5.11)$$\begin{gather} K^{(n)}=\left\{\begin{array}{ll} 0, & n=0,\\ -\hat{\eta}^{(n-1)}, & n=1,2,\ldots, \end{array}\right. \end{gather}$$

(5.12)$$\begin{gather}H^{(n)}=\left\{\begin{array}{ll} 0, & n=0,\\ -\hat{\xi}^{(n-1)}, & n=1,2,\ldots\,. \end{array}\right. \end{gather}$$

The leading-order equation in this context corresponds to the original linear problem associated with the zero eigenvalue. The solution to this equation is determined to be

(5.13)\begin{equation} \left(\hat{\eta}^{(0)},\hat{\xi}^{(0)}\right)=\left(\bar{\eta}_x,\bar{\xi}_x\right). \end{equation}

At $O(\lambda )$, the first-order perturbation solutions satisfy

(5.14a,b)\begin{equation} \mathcal{L}_1\left[\hat{\eta}^{(1)},\hat{\xi}^{(1)}\right]={-}\hat{\eta}^{(0)},\quad \mathcal{L}_2\left[\hat{\eta}^{(1)},\hat{\xi}^{(1)}\right]={-}\hat{\xi}^{(0)}. \end{equation}

System (5.14a,b) has a particular solution of $(\hat {\eta }^{(1)},\hat {\xi }^{(1)})=(-\bar {\eta }_c,-\bar {\xi }_c)$, which can be verified by taking derivatives of (5.1)–(5.2) with respect to $c$. At $O(\lambda ^2)$, the first-order correction forces the equations for $(\hat {\eta }^{(2)},\hat {\xi }^{(2)})$ as

(5.15a,b)\begin{equation} \mathcal{L}_1\left[\hat{\eta}^{(2)},\hat{\xi}^{(2)}\right]={-}\hat{\eta}^{(1)},\quad \mathcal{L}_2\left[\hat{\eta}^{(2)},\hat{\xi}^{(2)}\right]={-}\hat{\xi}^{(1)}. \end{equation}

In terms of the Fredholm alternative, a solvability condition for (5.15a,b) is that the vector $(\bar {\eta }_c,\bar {\xi }_c)$, composed of the forcing terms, is orthogonal to the solutions of the homogeneous adjoint problem. The adjoint problem is given by

(5.16a,b)\begin{equation} \mathcal{L}^+_1[f,g]=0,\quad\mathcal{L}^+_2[f,g]=0, \end{equation}

where

(5.17a,b)\begin{equation} \mathcal{L}^+_1[f,g]=\mathcal{L}_{21}[g]-\mathcal{L}_{22}[f],\quad \mathcal{L}^+_2[f,g]=\mathcal{L}_{12}[g]-\mathcal{L}_{11}[f], \end{equation}

formally akin to the inverse of the linear operators $\mathcal {L}_1$ and $\mathcal {L}_2$. This adjoint problem is solved by

(5.18a,b)\begin{equation} f={-}\bar{\xi}_x,\quad g=\bar{\eta}_x. \end{equation}

Thus, the corresponding solvability condition reads

(5.19)\begin{equation} \int_{-\infty}^{+\infty}(1+\bar{\eta})\left(\bar{\eta}_c\bar{\xi}_x-\bar{\xi}_x\bar{\eta}_c\right){{\rm d}{\kern 0.9pt}x}=\frac{1}{c}\frac{{\rm d}\mathcal{H}}{{\rm d}c}=0. \end{equation}

A necessary condition for the translational eigenmode to switch between stable and unstable states is that the wave energy, as a function of wave speed, has an extremum. This coincides with the necessary condition for the stability exchange of periodic water waves subjected to superharmonic perturbations (see, e.g. Tanaka 1983; Saffman 1985; Tanaka 1986; Zufiria & Saffman 1986; Kataoka 2006). The detailed proof of (5.19) is provided in Appendix B.

5.2. Stability of solitary waves – numerics

In this subsection, we examine numerically the superharmonic instability of locally confined travelling waves computed in § 4.2. This is achieved by integrating numerically the single evolution equation (4.4) with slightly perturbed solitary waves as the initial condition. Unlike periodic waves, whose stability properties are in some cases primarily affected by subharmonic perturbation, such as the Benjamin–Feir instability in water waves, it is sufficient to check superharmonic instability for solitary waves since the real axis can adapt to arbitrary harmonic perturbations. However, according to the necessary condition for the stability exchange subjected to superharmonic perturbations proved in § 5.1, the stability characteristics of solitary waves on a monotonic segment of the speed-energy bifurcation diagram should be the same. Therefore, it suffices to check the stability property of one typical wave between two successive stationary points on the speed-energy curve. In the subsequent numerical experiments, the mesh spacing and time step are typically chosen as ${\textrm {d}{\kern 0.9pt}x}=0.05$ and $\textrm {d}t=5\times 10^{-4}$ in most computations, and a frame of reference moving with the speed of the unperturbed solitary wave is adopted for convenience.

We begin by considering the solution branches of monotonic solitary waves bifurcating from zero amplitude. As shown in figure 8, each speed-energy bifurcation curve computed with the fully nonlinear Euler equations (black dashed line) features a stationary point marked with an asterisk where the energy reaches its maximum. Specifically, the elevation branch for $B=1.25$ achieves maximum energy of $\mathcal {H}^*=0.4138$ at $c^*=0.1766$, while the depression branch for $B=1.75$ peaks at $\mathcal {H}^*=0.2861$ and $c^*=0.3158$. However, for the depression branch shown in figure 8($b$), the speed-energy curve of the cubic full-dispersion model becomes significantly more deviated from the counterparts obtained from the Euler equations upon $c$ being less than (approximate) $0.3961$, a threshold marked by the red circle, and ultimately fails to capture the stationary point. We speculate that the discrepancy arises from the intense nonlinearity of solutions as they approach the bottom boundary $r=d$. In the spirit of Saffman's theory, we first verify numerically the stability of solitary waves for the elevation branch illustrated in figure 8($a$). For the region $c\leq 0.1766$, where solitary waves have a relatively large amplitude, an examination was carried out to observe the free propagation of a single hump travelling-wave solution with $c=0.1364$ and $S(0)=1.90$ in the presence of a $1\,\%$ amplitude-decreasing perturbation at $t=0$. The snapshots of the time evolution of the perturbed wave are shown in figure 9. The numerical solution was observed to undergo a noticeable transition to a smaller amplitude solitary wave as time evolves. The other elevation waves were also tested within the region, yielding a similar phenomenon. Based on these numerical experiments, we can conclude that monotonic elevation solitary waves falling within the range $0< c\leqslant 0.1766$ are unstable according to the generalized Saffman theory. Within the region $0.1766< c< c_0$, superharmonic perturbations with magnitudes ranging from $1\,\%$ to $5\,\%$ amplitude decrease of initial waves showed no instability. Figure 10(a) displays the time-evolution snapshots of a solitary wave subjected to a $2\,\%$ amplitude-decreasing perturbation. Apart from translation in the $x$-direction, the wave structure and amplitude remain barely unchanged from the initial wave, and the maximum difference between the two profiles is within the order of $O(10^{-3})$ throughout the time-evolution process. Finally, we can draw a conclusion for the elevation branch that the exchange stability does occur at the stationary point. Examining the stability property of the depression branch is analogous and straightforward, but the range $c<0.3961$ is excluded due to the limitation of the cubic full-dispersion model. Based on a series of numerical experiments, it has been verified numerically that solitary waves in the region $0.3961< c<0.546$ are stable. One illustrative example of the time-evolution snapshots, subjected to a $3\,\%$ amplitude-decreasing perturbation, is presented in figure 10(b).

Figure 8. Speed-energy bifurcation curves with fixed $d=5/11$ and varying $B$: ($a$) the elevation branch with $B=1.25$ and ($b$) the depression branch with $B=1.75$. The solution branches computed with the full Euler equations and the cubic full-dispersion model are shown with black dashed lines and blue solid lines, respectively. The two insets demonstrate the stationary points obtained in the full Euler equations.

Figure 9. Time evolution of a perturbed elevation wave with $c=0.1364< c^*=0.1766$. The snapshots (blue solid lines) represent local wave profiles at $t=0$, 400, 800, 1200, 1600 and 2000, and the $t-S$ curve (black dash-dotted line) shows the maximum amplitude variations over time.

Figure 10. Time evolution of perturbed solitary waves with fixed $d=5/11$ and varying $B$. ($a$) Snapshots of an elevation wave with $B=1.25$ and $c=0.2489>c^{*}=0.1766$, subjected to the $2\,\%$ amplitude-decreasing perturbation. ($b$) Snapshots of a depression wave with $B=1.75$ and $c=0.4077$, subjected to the $3\,\%$ amplitude-decreasing perturbation. ($c$) Snapshots of a wavepacket depression wave with $B=28$ and $c=2.603$, subjected to the $2\,\%$ amplitude-decreasing perturbation. The perturbed initial profiles are plotted with dash-dotted lines, and the terminal profiles at $t=3000$, 4000 and 5000 are shown with solid lines from top to bottom, respectively. In each panel, the maximum difference between the two profiles is within the order of $O(10^{-3})$.

Moreover, we carried out numerical stability tests for monotonic solitary waves bifurcating from finite amplitudes. In figures 11(a) and 11(b), we also observed stationary points in the speed-energy bifurcation curves obtained with the full Euler equations. The cubic full-dispersion model behaves better for the elevation branch than the depression branch since, in the former situation, the model captures the stationary point, but it fails in the latter case. The numerical experiments show that the elevation branch computed with the cubic model also experiences an exchange of stability near the stationary point where the maximum energy is $\mathcal {H}^*= 0.7735$ at $c^{*}= 0.1948$, as shown in figure 11($a$). The time-evolution dynamical behaviours of these waves closely resemble those of their counterparts bifurcating from zero amplitude and thus will not be presented again here.

Figure 11. Speed-energy curves computed with fixed $d=5/11$ and varying $B$: ($a$) elevation branch with $B=1.3$; ($b$) depression branch with $B=1.25$ and ($c$) wavepacket branches with $B=28$. Both branches presented in panels ($a$) and ($b$) bifurcate from finite amplitudes. The solution branches computed with the full Euler equations and the cubic full-dispersion model are shown with black dashed lines and blue solid lines, respectively.

Subsequently, we move on to the stability characteristics associated with wavepacket solitary waves, whose speed-energy curves are monotonic, as shown in figure 11($c$). All the computational evidence suggests that the whole elevation branch is unstable and the depression one is stable. Two illustrative examples are shown in figures 10($c$) and 12. The time-evolution snapshots of a stable depression solitary wave with $c=2.603$ and an initial 2 % amplitude-decreasing perturbation are presented in figure 10($c$). Obviously, the snapshot observed at $t=5000$ exhibits a striking similarity to the initial profile, resulting in a difference that is nearly indistinguishable. However, as shown in figure 12, an exact elevation wave evolves rapidly into two depression waves of distinct amplitudes as time goes on due to numerical errors. As demonstrated in figure 12, for the long-time dynamics, the depression wave with the smaller amplitude travels rightward more quickly than the larger one, so an overtaking collision occurs.

Figure 12. Time evolution of a wavepacket elevation wave for $c=2.68$, $d=5/11$ and $B=28$. The snapshots of wave profiles are shown at $t=0$, 30, 60, 90 and 120 from top to bottom.

5.3. Dynamics of solitary waves

In this subsection, we supplement the stability analysis by investigating the interaction between two stable solitary waves. In all the collision experiments, the initial data were constructed by properly shifting the existing travelling waves so as to minimize their overlap and then adding them together. The direction of wave propagation (left or right) can be set initially by simply altering the sign of the velocity potential $\xi$. To quantify the degree of the inelastic interactions, two fundamental quantities can be recorded: the energy loss and the amplitude variation of solitary waves, with the former measured by the ratio of the residual energy $\mathcal {H}_R$ to the total energy $\mathcal {H}_T$, as discussed by Craig et al. (2006). Following Craig et al. (2006), the residual is calculated by matching the well-separated solitary waves after collision with computed travelling-wave solutions and subtracting them from the time-dependent solution. Note that time-evolution simulations presented here extend those reported by Guyenne & Părău (2016) to include more types of colliding solitary waves.

The first set of computational experiments are head-on collisions. Figures 13 and 14 present several illustrative scenarios in which both colliding waves bifurcate from zero amplitude. These numerical results indicate that both waves persist after the collision, and the noticeable residual waves, induced by the interaction and occurring ahead of the two separating pulses, suggest the inelasticity of the collision. Furthermore, we investigate the head-on collision between stable elevation and depression waves, one on the branch bifurcating from zero amplitude and the other on the branch bifurcating from finite amplitude. As shown in figure 15, two numerical experiments, $B=1.25$ and $B=1.3$, were conducted. For $B=1.25$, figure 15($a$) shows an illustrative example of head-on interaction for the speeds of $c=0.2698$ (a right-going wave on the branch bifurcating from zero amplitude) and $c=0.3076$ (a left-going wave on the branch bifurcating from finite amplitude) at $t=0$. The amplitudes of solitary waves before (after) the collision are $0.4$ and $-0.2$ ($0.3989$ and $-0.1622$ at $t=350$). At $t=350$, the total energy was calculated to be $\mathcal {H}_T=0.4147$ and the energy of residual waves to be $\mathcal {H}_R=0.0198$, giving rise to $\mathcal {H}_T/\mathcal {H}_R\approx 4.77\,\%$. Figure 15($b$) depicts an alternative dynamic evolution scenario at the magnetic Bond number $B=1.30$, wherein the bifurcation role of the two solution branches is reversed compared with the former experiment. The time-evolution snapshots of the head-on interaction for the speeds of $c=0.3269$ (travelling rightwards) and $c=0.3376$ (travelling leftwards) at $t=0$ are presented. The amplitudes of the two solitary waves before (after) the collision are $0.3$ and $-0.15$ ($0.2968$ and $-0.1341$ at $t=350$). At $t=350$, the total energy was calculated to be $\mathcal {H}_T=0.4038$, and the energy of residual waves to be $\mathcal {H}_R=0.0139$, leading to $\mathcal {H}_T/\mathcal {H}_R\approx 3.44\,\%$. It can be observed that the two depression waves undergo an apparent decrease in amplitude after the interaction.

Figure 13. ($a$) Asymmetric head-on collision between two monotonic depression solitary waves with the initial speeds $c=0.4935$ (travelling rightwards) and $c=0.3625$ (travelling leftwards), respectively, for $d=5/11$ and $B=1.75$. ($b$) Snapshots of wave profiles at $t=90$, 117 and 140 from top to bottom are presented, and the dotted line in the middle panel denotes the rigid boundary.

Figure 14. Head-on collisions between wavepacket depression solitary waves for $B=28$ and $d=5/11$. ($a$) Snapshots of collision between depression waves of different amplitudes, $S(0)=0.9$ ($c=2.758$) and $S(0)=0.85$ ($c=2.5638$), propagating in opposite directions, at $t=0$, 1, 2 and 3 from top to bottom. ($b$) Snapshots of collision between depression waves of identical amplitude, $S(0)=0.85$ ($c=2.5638$), propagating in opposite directions, at $t=0$, 1, 2 and 3 from top to bottom.

Figure 15. Head-on collisions between monotonic solitary waves of opposite polarities for ($a$) $B=1.25$ and ($b$) $B=1.30$, together with $d=5/11$. $(a)$ Snapshots of collision between an elevation wave on the branch bifurcating from zero amplitude ($S(0)=1.40$ and $c=0.2698$) and a depression wave on the branch bifurcating from finite amplitude ($S(0)=0.8$ and $c=0.3076$), at $t=0$, 100, 150, 200 and 350 from top to bottom. ($b$) Snapshots of collision between an elevation wave (right-going wave on the branch bifurcating from finite amplitude with $S(0)=1.30$ and $c=0.3269$) and a depression wave (left-going wave on the branch bifurcating from zero amplitude with $S(0)=0.85$ and $c=0.3376$), at $t=0$, 100, 150, 200 and 350 from top to bottom.

The second set of computational experiments is overtaking collisions where two waves propagate in the same direction, with a larger and slower wave positioned in front of a faster and smaller wave. Several illustrative cases are shown in figures 16–18. In contrast with the head-on collision, the inelasticity of the overtaking interaction is relatively weak and produces residual waves that are small, even barely noticeable in some instances. A feature of these interactions observed is that after the collision, the larger solitary wave is of slightly larger amplitude than its initial counterpart, which has already been reported in the water-wave problem (Craig et al. 2006). As shown in figure 18($a$), the amplitudes of the two waves change from $0.4$ and $-0.2$ (before the interaction) to $0.4004$ and $-0.1961$ (after the interaction), and a non-zero energy transfer from the incident waves to the residual results in the ratio of the residual to the total energy $\mathcal {H}_R/\mathcal {H}_T$ being approximately $0.41\,\%$ at $t=3500$.

Figure 16. Overtaking collision between monotonic depression solitary waves for $d=5/11$ and $B=1.75$. Two initial waves are characterized by $S(0)=0.8$, $c=0.4935$ and $S(0)=0.50$, $c=0.2405$. ($a$) Time-dependent solution in the $x\unicode{x2013} t$ plane. ($b$) Snapshots at $t=50$, 95 and 150 from top to bottom.

Figure 17. Overtaking collision between wavepacket depression solitary waves for $d=5/11$ and $B=28$. Two initial waves are characterized by $S(0)=0.9$, $c=2.758$ and $S(0)=0.85$, $c=2.5638$. The snapshots at $t =0$, 20, 24, 26, 30 and 40 are shown in a frame of reference moving with speed $c=2.5638$.

Figure 18. Overtaking collision between monotonic solitary waves of different polarities for ($a$) $B=1.25$ and ($b$) $B=1.3$, together with $d=5/11$. $(a)$ Snapshots of collision between a depression wave on the branch bifurcating from finite amplitude ($S(0)=0.8$, $c=0.3076$) and an elevation wave on the branch bifurcating from zero amplitude ($S(0)=1.40$, $c=0.2698$), at $t=0$, 1500, 2000, 2500 and 3500 from top to bottom. ($b$) Snapshots of collision between a depression wave on the branch bifurcating from zero amplitude ($S(0)=0.9$, $c=0.3426$) and an elevation wave on the branch bifurcating from finite amplitude ($S(0)=1.50$, $c=0.3018$), at $t=0$, 1000, 2000, 3000 and 4000 from top to bottom. The snapshots in panels ($a$) and ($b$) are shown in the respective frames of reference moving with speeds $c=0.2698$ and $c=0.3018$.

6. Conclusions

In this paper, nonlinear waves propagating on the surface of an axisymmetric ferrofluid jet under an azimuthal magnetic field have been investigated from asymptotic and numerical perspectives. The problem has been proved to be a Hamiltonian system, and a new definition has been introduced for the crucial pseudo-differential operator involved (i.e. the Dirichlet–Neumann operator in the axisymmetric geometry). Subsequently, multi-scale models with weak or strong nonlinearity in various asymptotic limits have been derived systematically based on the Hamiltonian structure and the expansions of the Dirichlet–Neumann operator. Through a comparison with the full Euler equations and other reduced models (the KdV equation, the fifth-order KdV equation, the modified KdV equation and the strongly nonlinear long-wave model), such as wave profiles and bifurcation diagrams, the cubic full-dispersion model can be considered as the quantitatively accurate and optimally efficient reduced model to accommodate both monotonic and wavepacket solitary waves found previously in the ferrofluid jet problem.

The stability and dynamics of solitary-wave solutions have also been explored theoretically and numerically based on the cubic full-dispersion model. First, the stability characteristics of solitary waves subjected to longitudinal superharmonic perturbations have been discussed analytically through an asymptotic procedure, which broadens the scope of Saffman's theory valid for axisymmetric solitary waves. Subsequently, a time-dependent numerical scheme has been implemented using the accurate pseudo-spectral method coupled with a fourth-order Runge–Kutta scheme to simulate the long-time propagation and pairwise collisions of solitary waves. All the numerical experiments we tested reveal that the monotonic elevation waves on a branch bifurcating from zero amplitude (or finite amplitude) experience stability exchange near the stationary point on the speed-energy curve (denoted by $c^*$). When $c>c^*$, the elevation waves remain stable, but for $c< c^*$, they become unstable and undergo a transition to different amplitudes during the long-time propagation. Monotonic solitary waves on the depression branch appear stable as they are in the region of $c>c^*$. However, when $c< c^*$, the cubic full-dispersion model fails to depict the phenomenon of stability exchange due to the complicated nonlinearity while approaching the rigid bottom. Concurrently, wavepacket depression waves have been observed to exhibit stability, while wavepacket elevation waves have displayed instability and evolved into two wavepacket depression waves. Finally, in all numerical experiments that we have considered, head-on and overtaking collisions are both inelastic, generating residual waves ahead of the separating pulses.

Our primary results from the cubic full-dispersion model motivate further investigations with the full Euler equations. First, it would be interesting to employ an asymptotic analysis for small eigenvalues and explicitly give their growth rates near the critical point where an exchange of stability occurs. Second, in terms of the limitation of the cubic model, it is also necessary to develop a new time-dependent numerical strategy for the Euler equations to thoroughly investigate the time-evolution characteristics of solitary waves, especially involving highly nonlinear wave motions, such as the stability and dynamics of overhanging waves. Finally, the bifurcation mechanism for monotonic solitary waves of ferrofluid jets is extraordinary, which cannot always be captured by the cubic model (especially when the magnetic Bond number is considerably large). This bifurcation mechanism merits a thorough study, and proposing a suitable reduced model is particularly interesting.