2 Formulation

We consider an axisymmetric column of ferrofluid with constant density $\unicode[STIX]{x1D70C}_{1}$ and magnetic susceptibility $\unicode[STIX]{x1D712}_{1}$ , coating a copper rod of radius $d$ . We choose the cylindrical coordinate system $(x,r)$ such that $x$ points along the rod, and $r$ is the radial coordinate. The ferrofluid is surrounded by a non-magnetisable fluid ( $\unicode[STIX]{x1D712}_{2}=0$ ) of density $\unicode[STIX]{x1D70C}_{2}\leqslant \unicode[STIX]{x1D70C}_{1}$ (we do not consider values of $\unicode[STIX]{x1D70C}_{2}>\unicode[STIX]{x1D70C}_{1}$ to avoid the Rayleigh–Taylor instability). The interface is given by $r=\unicode[STIX]{x1D702}(x,t)$ , the mean radius of which is denoted $R$ . Denote the velocity fields in the ferrofluid and surrounding fluid as $\boldsymbol{u}_{1}=(u_{1},v_{1})$ and $\boldsymbol{u}_{2}=(u_{2},v_{2})$ in $(x,r)$ respectively. The system is contained inside a fixed cylindrical container of radius $D$ (see figures 1 and 2). We note that in the experiments of Arkhipenko et al. (Reference Arkhipenko, Barkov, Bashtovoi and Krakov1980) and Bourdin et al. (Reference Bourdin, Bacri and Falcon2010), the fluids were contained in a rectangular box. However, since axisymmetric interfaces were witnessed, the box must have been of a sufficient size to not destroy the axisymmetry of the problem. Therefore, comparison between the experiments and the model in this paper can be made by considering large values of $D$ .

A current $I$ is passed through the copper wire. This induces a purely azimuthal magnetic field, given by

(2.1) $$\begin{eqnarray}\boldsymbol{H}=\unicode[STIX]{x1D707}_{0}I/(2\unicode[STIX]{x03C0}r)\boldsymbol{e}_{\unicode[STIX]{x1D703}},\end{eqnarray}$$

where $\unicode[STIX]{x1D707}_{0}$ is the magnetic permeability in a vacuum, and $e_{\unicode[STIX]{x1D703}}$ is the unit vector in the clockwise azimuthal direction. We assume the linear magnetisation law, such that the magnetisation $\boldsymbol{M}$ satisfies $\boldsymbol{M}=\unicode[STIX]{x1D712}\boldsymbol{H}$ . The assumption of an axisymmetric interface results in the decoupling of the magnetic problem from the hydrodynamic problem (Rannacher & Engel Reference Rannacher and Engel2006). Continuity of pressure (Rosensweig Reference Rosensweig1985, § 5.1) on the interface is given by

(2.2) $$\begin{eqnarray}P_{1}=P_{2}+T\unicode[STIX]{x1D705}-\frac{\unicode[STIX]{x1D707}_{0}}{2}\unicode[STIX]{x1D712}^{2}(\boldsymbol{H}\boldsymbol{\cdot }\hat{\boldsymbol{n}}).\end{eqnarray}$$

Here, $P_{1}$ and $P_{2}$ are the pressures in the ferrofluid and outer fluid respectively, $T$ the surface tension, and $\unicode[STIX]{x1D705}$ the mean curvature, given by

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D705}=\frac{1}{\unicode[STIX]{x1D702}}\left(1+\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}{\unicode[STIX]{x2202}x}\right)^{2}\right)^{-1/2}-\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D702}}{\unicode[STIX]{x2202}x^{2}}\left(1+\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}{\unicode[STIX]{x2202}x}\right)^{2}\right)^{-3/2}.\end{eqnarray}$$

Note that since $\boldsymbol{H}$ is azimuthal, the pressure jump associated with the magnetic field is zero.

Figure 1.

Formulation in cylindrical coordinates.

Figure 2.

Three-dimensional visualisation of the problem.

We consider a wave of unchanging form with wavelength $\unicode[STIX]{x1D706}$ and celerity $c$ . Under the assumption that the flows in either region are irrotational and incompressible, both velocity fields can be written in terms of a velocity potential $\boldsymbol{u}_{1,2}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{1,2}$ , where $\unicode[STIX]{x1D719}_{1}$ and $\unicode[STIX]{x1D719}_{2}$ satisfy the axisymmetric Laplace equation, given by

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}_{i}=\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{i}}{\unicode[STIX]{x2202}r^{2}}+\frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{i}}{\unicode[STIX]{x2202}r}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{i}}{\unicode[STIX]{x2202}z^{2}}=0,\quad i=1,2,\end{eqnarray}$$

in their respective flow domains. We assume the wave is symmetric about the point $\unicode[STIX]{x1D719}_{1}=\unicode[STIX]{x1D719}_{2}=0$ . We require no normal flow through the rod and outer cylinder, that is

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}r}=0\quad \text{for}~r=d, & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}r}=0\quad \text{for}~r=D. & \displaystyle\end{eqnarray}$$

The Bernoulli principle (Rosensweig Reference Rosensweig1985, § 5.2) satisfied on the interface gives

(2.7) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D70C}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}t}+\frac{1}{2}(q_{1}^{2}-\unicode[STIX]{x1D70C}q_{2}^{2})+\frac{1}{\unicode[STIX]{x1D70C}_{1}}(P_{1}-P_{2})-\frac{\unicode[STIX]{x1D707}_{0}\unicode[STIX]{x1D712}_{1}I^{2}}{4\unicode[STIX]{x03C0}^{2}\unicode[STIX]{x1D70C}_{1}}\frac{1}{2r^{2}}={\hat{C}},\end{eqnarray}$$

where $q_{i}=|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{i}|$ , $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{2}/\unicode[STIX]{x1D70C}_{1}$ and ${\hat{C}}$ is the Bernoulli constant. We take $R$ as the reference length and $\sqrt{T/(R\unicode[STIX]{x1D70C}_{1})}$ as the reference velocity. Making use of (2.2), we find that the non-dimensionalised Bernoulli equation is

(2.8) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D70C}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}t}+\frac{1}{2}(q_{1}^{2}-\unicode[STIX]{x1D70C}q_{2}^{2})+\unicode[STIX]{x1D705}-\frac{B}{2r^{2}}=C,\end{eqnarray}$$

where the magnetic Bond number $B$ is defined as

(2.9) $$\begin{eqnarray}B=\frac{\unicode[STIX]{x1D707}_{0}\unicode[STIX]{x1D712}_{1}I^{2}}{4\unicode[STIX]{x03C0}^{2}RT}.\end{eqnarray}$$

The magnetic Bond number is a ratio of magnetic to surface tension forces. It is shown in § 3 that the stability of linear perturbations is determined by $B$ . Finally, the kinematic boundary condition on the interface is given by

(2.10) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{i}}{\unicode[STIX]{x2202}x}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{i}}{\unicode[STIX]{x2202}r},\quad i=1,2.\end{eqnarray}$$

Note that for solitary waves with a flat far field, instead of fixing the mean of $\unicode[STIX]{x1D702}$ to unity, we fix $\unicode[STIX]{x1D702}$ in the far field to be unity. This choice of scaling gives rise to the far-field condition

(2.11) $$\begin{eqnarray}\unicode[STIX]{x1D702}\rightarrow 1,\quad \text{as}~x\rightarrow \pm \infty .\end{eqnarray}$$

It is left to solve the governing equation (2.4) for $\unicode[STIX]{x1D719}_{1}$ and $\unicode[STIX]{x1D719}_{2}$ in their respective flow domains, subject to boundary conditions, (2.5), (2.6), (2.8) and (2.10). We consider two values of $\unicode[STIX]{x1D70C}$ , that is $\unicode[STIX]{x1D70C}=0$ (one-layer model) and $\unicode[STIX]{x1D70C}=1$ (two-layer model). For the one-layer model, we ignore the outer boundary $r=D$ . This removes the requirement to solve for $\unicode[STIX]{x1D719}_{2}$ , since the equations concerning just $\unicode[STIX]{x1D719}_{1}$ form a closed system (no $\unicode[STIX]{x1D719}_{2}$ terms are present in (2.8) when $\unicode[STIX]{x1D70C}=0$ ).

In the following section, we derive the linear dispersion relation for the system.

3 Linear theory

Consider a small perturbation to the uniform stream of the form

(3.1a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{1}=\unicode[STIX]{x1D716}\mathop{\sum }_{m=1}^{\infty }F_{m}(r)\text{e}^{\text{i}mk(x-ct)},\quad \unicode[STIX]{x1D719}_{2}=\unicode[STIX]{x1D716}\mathop{\sum }_{m=1}^{\infty }G_{m}(r)\text{e}^{\text{i}mk(x-ct)},\end{eqnarray}$$

where $|\unicode[STIX]{x1D716}|\ll 1$ , and $F_{m}$ and $G_{m}$ are unknown functions of $r$ . Note that if $c^{2}>0$ , the solution is stable, while if $c^{2}<0$ , the amplitude grows exponentially in time and the solution is unstable. Ignoring terms of $O(\unicode[STIX]{x1D716}^{2})$ , and solving the linearised system, one finds the equation for the free surface,

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D702}=1+C_{1}\unicode[STIX]{x1D716}\left(\text{I}_{1}(k)-\frac{\text{I}_{1}(kd)}{\text{K}_{1}(kd)}\text{K}_{1}(k)\right)\text{e}^{\text{i}k(x-ct)}.\end{eqnarray}$$

Here, $\text{I}_{n}$ and $\text{K}_{n}$ are the modified Bessel functions of the first and second kind of order $n$ , and $C_{1}$ is an arbitrary constant. Equation (3.2) is a linear perturbation of wavenumber $k$ , travelling at speed $c$ . Furthermore, we recover the linear dispersion relation

(3.3) $$\begin{eqnarray}c^{2}=\frac{1}{\displaystyle k\left(\frac{m_{2}^{d}}{m_{1}^{d}}-\unicode[STIX]{x1D70C}\frac{m_{2}^{D}}{m_{1}^{D}}\right)}(k^{2}-1+B),\end{eqnarray}$$

where

(3.4a,b ) $$\begin{eqnarray}m_{1}^{d}=\text{I}_{1}(k)\text{K}_{1}(kd)-\text{K}_{1}(k)\text{I}_{1}(kd),\quad m_{2}^{d}=\text{I}_{0}(k)\text{K}_{1}(kd)+\text{I}_{1}(kd)\text{K}_{0}(k).\end{eqnarray}$$

Replacing all instances of $d$ with $D$ in the above equations gives $m_{1}^{D}$ and $m_{2}^{D}$ . If it is the case that $c(k)=c(nk)$ for some positive integer $n\geqslant 2$ , then the leading-order solution is given by

(3.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702} & = & \displaystyle 1+\unicode[STIX]{x1D716}C_{1}\left(\text{I}_{1}(k)-\frac{\text{I}_{1}(kd)}{\text{K}_{1}(kd)}\text{K}_{1}(k)\right)\text{e}^{\text{i}k(x-ct)}\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D716}C_{n}\left(\text{I}_{1}(nk)-\frac{\text{I}_{1}(nkd)}{\text{K}_{1}(nkd)}\text{K}_{1}(nk)\right)\text{e}^{\text{i}nk(x-ct)}.\end{eqnarray}$$

This phenomenon is called Wilton ripples, and is only possible when a minimum occurs in the dispersion relation (corresponding to $B>B_{2}$ ). A similar property is seen in gravity–capillary waves, and was originally derived by Wilton (Reference Wilton1915). One finds at higher order a solvability condition for $C_{n}$ . However, the algebra quickly becomes complicated, and instead these solutions are recovered via fully nonlinear computations, as seen in § 6.

Since $d<1<D$ and $0\leqslant \unicode[STIX]{x1D70C}\leqslant 1$ , the denominator in (3.3) is always positive, meaning that the stability of the solution depends on $k$ and $B$ . We find that solutions with wavenumber $k$ are stable if

(3.6) $$\begin{eqnarray}k^{2}>1-B.\end{eqnarray}$$

This is true for all $k$ if $B>1$ . Note that we recover the stability condition found by Rayleigh (Reference Rayleigh1878) by taking $B=0$ (that all solutions with $k<1$ are unstable).

The right-hand side of (3.3) tends to infinity as $k\rightarrow \infty$ . Hence, whether the dispersion curve has a minimum or not can be determined by considering the gradient of $c^{2}$ for small $k$ . A negative gradient for small $k$ corresponds to the existence of a minimum. Denoting the dispersion relation when $\unicode[STIX]{x1D70C}=0$ as $c_{\unicode[STIX]{x1D70C}}$ , we find that

(3.7) $$\begin{eqnarray}c_{\unicode[STIX]{x1D70C}}^{2}=\frac{1}{k}\left(\frac{m_{1}^{d}}{m_{2}^{d}}\right)(k^{2}-1+B).\end{eqnarray}$$

Taking a small $k$ expansion of the above equation, and differentiating with respect to $k$ , one gets

(3.8) $$\begin{eqnarray}2c_{\unicode[STIX]{x1D70C}}\frac{\text{d}c_{\unicode[STIX]{x1D70C}}}{\text{d}k}\approx \frac{1}{8}[(-1+4d^{2}-3d^{2}+4d^{4}\log d)(B-1)+8(1-d^{2})]k+O(k^{3}).\end{eqnarray}$$

Hence, there exists a minimum in $c_{\unicode[STIX]{x1D70C}}(k)$ given that the coefficient of $k$ in the above equation is negative. This is the case if $B>B_{2}$ , where $B_{2}$ has the following dependence on $d$ :

(3.9) $$\begin{eqnarray}B_{2}(d)=1+\frac{8(1-d^{2})}{1-4d^{2}+3d^{4}-4d^{4}\log d}.\end{eqnarray}$$

This expression is in agreement with (3.5) in the paper of BP. We see in § 6 that the characteristics of the solution space changes upon the existence of a minimum.

When $\unicode[STIX]{x1D70C}=1$ , we now expect $B_{2}$ to have dependence on both $d$ and $D$ . In the case when $D\rightarrow \infty$ , BP demonstrated that $B_{2}=1$ . For a finite value of $D$ , one can follow the same argument given above for $c_{\unicode[STIX]{x1D70C}}$ to find that

(3.10) $$\begin{eqnarray}B_{2}(d,D)=1-\frac{E}{F},\end{eqnarray}$$

where

(3.11) $$\begin{eqnarray}\displaystyle E & = & \displaystyle \frac{\frac{1}{2}(1-d^{2})(D^{2}-1)}{(D^{2}-1)+(1-d^{2})},\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle F & = & \displaystyle \frac{1}{8(D^{2}-d^{2})^{2}} [(D^{2}-1)(1-d^{2})(D-d)(D+d)+2d^{4}(D^{2}-1)^{2}\log d\nonumber\\ \displaystyle & & \displaystyle -\,2D^{4}(d^{2}-1)^{2}\log D].\end{eqnarray}$$

We will find it useful to denote the value of $c$ at $k=0$ as $c_{0}$ , and the minimum value of $c$ occurring at $k=k_{m}$ to be denoted $c_{m}$ . When $B<B_{2}$ , $c_{0}=c_{m}$ . In the following section, we describe the numerical method used to solve the fully nonlinear problem.

4 Numerical scheme

We consider a wave of wavelength $\unicode[STIX]{x1D706}$ travelling with unchanging form at a constant speed $c$ . We remove time dependence by taking a frame of reference travelling with the wave. We will use a finite difference scheme originally proposed by Woods (Reference Woods1951) for axisymmetric flows, and later independently formulated by Jeppson (Reference Jeppson1970). The method has since been used by a variety of authors for axisymmetric capillary waves (Vanden-Broeck et al. Reference Vanden-Broeck, Miloh and Spivack1998) under the effect of electric (Grandison et al. Reference Grandison, Vanden-Broeck, Papageorgiou, Miloh and Spivak2008) or magnetic (BP) fields, and the rise of Taylor bubbles in a tube (Doak & Vanden-Broeck Reference Doak and Vanden-Broeck2018). We will first describe the method used to find solutions to the two-layer model. This involves adapting the finite difference scheme to allow for two computational domains, as described below. Following this, we state the simplifications made to the method to solve the one-layer problem.

The idea is to solve the problem by finding the physical variable $r$ in the two potential spaces $(\unicode[STIX]{x1D719}_{1},\unicode[STIX]{x1D713}_{1})$ and $(\unicode[STIX]{x1D719}_{2},\unicode[STIX]{x1D713}_{2})$ , where $\unicode[STIX]{x1D713}_{1}$ and $\unicode[STIX]{x1D713}_{2}$ are the Stokes streamfunctions, defined by

(4.1a,b ) $$\begin{eqnarray}u_{i}=\frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}}{\unicode[STIX]{x2202}r},\quad v_{i}=-\frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}}{\unicode[STIX]{x2202}x},\quad i=1,2.\end{eqnarray}$$

Lines given by $\unicode[STIX]{x1D713}_{i}=$ constant are everywhere parallel to the velocity vector $\boldsymbol{u}_{i}$ , and are orthogonal to lines of constant $\unicode[STIX]{x1D719}_{i}$ . However, it is interesting to note that the Stokes streamfunction, unlike it’s two-dimensional counterpart, does not satisfy the Laplace equation. Therefore, the powerful tools of complex analysis are unavailable to us here, since the mapping from the $(x,r)$ space to the $(\unicode[STIX]{x1D719},\unicode[STIX]{x1D713})$ space is not conformal. Without loss of generality, we choose to define $\unicode[STIX]{x1D713}_{1}=d^{2}c/2=Q_{d}$ on $r=d$ , $\unicode[STIX]{x1D713}_{1}=\unicode[STIX]{x1D713}_{2}=Q$ on the interface, and $\unicode[STIX]{x1D713}_{2}=Q_{D}$ on $r=D$ . We note that, in the case of a flat free surface (uniform stream solution), $Q=c/2$ and $Q_{D}=cD^{2}/2$ . Integrating (4.1) with $u_{i}=c$ and $v_{i}=0$ , the uniform stream solution is found to be

(4.2) $$\begin{eqnarray}r=\left\{\begin{array}{@{}ll@{}}\displaystyle \sqrt{\frac{2\unicode[STIX]{x1D713}_{1}}{c}},\quad & \text{if }Q_{d}\leqslant \unicode[STIX]{x1D713}_{1}\leqslant Q,\\ \displaystyle \sqrt{\frac{2\unicode[STIX]{x1D713}_{2}}{c}},\quad & \text{if }Q<\unicode[STIX]{x1D713}_{2}\leqslant Q_{D}.\end{array}\right.\end{eqnarray}$$

This encourages the coordinate transformations $\unicode[STIX]{x1D713}_{1}=t^{2}$ and $\unicode[STIX]{x1D713}_{2}=s^{2}$ to better distribute streamlines between the interface and the boundaries. This choice of transformation means that taking equally spaced points in the discretisation of $t$ and $s$ results in equally spaced streamlines in the computation of the uniform stream solution. Seeking a periodic wave of wavelength $\unicode[STIX]{x1D706}$ , symmetric about $\unicode[STIX]{x1D719}_{1}=\unicode[STIX]{x1D719}_{2}=0$ , the ferrofluid and surrounding fluid flow domains are mapped onto the rectangular domains $\unicode[STIX]{x1D6FA}_{1}$ and $\unicode[STIX]{x1D6FA}_{2}$ respectively, where

(4.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FA}_{1} & = & \displaystyle \{\unicode[STIX]{x1D719}_{1}\in [-c\unicode[STIX]{x1D706}/2,0],t\in [Q_{d}^{1/2},Q^{1/2}]\},\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FA}_{2} & = & \displaystyle \{\unicode[STIX]{x1D719}_{2}\in [-c\unicode[STIX]{x1D706}/2,0],s\in [Q^{1/2},Q_{D}^{1/2}]\}.\end{eqnarray}$$

Figure 3.

The two flow domains in potential space. The interface between the two fluids is in bold, and corresponds to the same streamline in physical space.

Here, we only consider the flow domains over half a wavelength, making use of the assumed symmetry. The flow domain in the potential space is shown in figure 3. Seeking $r$ as a function of the independent variables $(\unicode[STIX]{x1D719}_{1},\unicode[STIX]{x1D713}_{1})$ in $\unicode[STIX]{x1D6FA}_{1}$ and $(\unicode[STIX]{x1D719}_{2},\unicode[STIX]{x1D713}_{2})$ in $\unicode[STIX]{x1D6FA}_{2}$ , we find that (2.4) under the mapping becomes

(4.5) $$\begin{eqnarray}r^{3}\frac{\unicode[STIX]{x2202}^{2}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}^{2}}+r\frac{\unicode[STIX]{x2202}^{2}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{i}^{2}}+r^{2}\left(\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}}\right)^{2}-\left(\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{i}}\right)^{2}=0.\end{eqnarray}$$

Furthermore, one can express $q_{i}=|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{i}|$ and the mean curvature $\unicode[STIX]{x1D705}$ evaluated on the interface as functions of $\unicode[STIX]{x1D719}_{i}$ , using the identities

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle q_{i}(\unicode[STIX]{x1D719}_{i})=(u_{i}^{2}+v_{i}^{2})^{1/2}=\left(\left(\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{i}}\right)^{2}+r^{2}\left(\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}}\right)^{2}\right)^{1/2}, & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D705}_{i}(\unicode[STIX]{x1D719}_{i})=-q_{i}^{3}\left(r\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}}\frac{\unicode[STIX]{x2202}^{2}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{i}^{2}}-\left(\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{i}}\right)^{2}\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}}-r\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{i}}\frac{\unicode[STIX]{x2202}^{2}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{i}\unicode[STIX]{x1D713}_{i}}\right)+\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}}q_{i}. & \displaystyle\end{eqnarray}$$

Note that $\unicode[STIX]{x1D705}_{1}$ here denotes the mean curvature as a function of $\unicode[STIX]{x1D719}_{1}$ , and likewise for $\unicode[STIX]{x1D705}_{2}$ . These functions correspond to the same curve in physical space (the interface), and hence have the same value at given points along the interface, but are different functions due to the discontinuity in tangential velocities across the interface. We discretise $\unicode[STIX]{x1D6FA}_{1}$ and $\unicode[STIX]{x1D6FA}_{2}$ into equidistant points with $M$ points in $\unicode[STIX]{x1D719}_{1}$ and $\unicode[STIX]{x1D719}_{2}$ , $N$ points in $t$ and $P$ points in $s$ as follows

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{1_{i}}=\unicode[STIX]{x1D719}_{2_{i}}=-\frac{c\unicode[STIX]{x1D706}}{2(M-1)}(M-i)\quad i=1,\ldots ,M, & \displaystyle\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle t_{j}=Q_{d}^{1/2}+(Q^{1/2}-Q_{d}^{1/2})\frac{j-1}{N-1}\quad j=1,\ldots ,N, & \displaystyle\end{eqnarray}$$

(4.10) $$\begin{eqnarray}\displaystyle & \displaystyle s_{j}=Q^{1/2}+(Q_{D}^{1/2}-Q^{1/2})\frac{j-1}{P-1}\quad j=1,\ldots ,P. & \displaystyle\end{eqnarray}$$

We satisfy the governing equation (4.5) at the interior nodes of $\unicode[STIX]{x1D6FA}_{1}$ and $\unicode[STIX]{x1D6FA}_{2}$ , finding the values of derivatives with finite difference approximations. We use second-order central differences, making use of the symmetry by imposing $\unicode[STIX]{x2202}r/\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{i}=0$ at $\unicode[STIX]{x1D719}_{i}=0$ and $\unicode[STIX]{x1D719}_{i}=-c\unicode[STIX]{x1D706}/2$ (for $i=1,2$ ). On the interface, we use second-order backwards differences to compute derivatives with respect to $t$ , and forward differences for derivatives with respect to $s$ . Derivatives with respect to $\unicode[STIX]{x1D713}_{1}$ are given in terms of derivatives with respect to $t$ via the identities

(4.11) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}} & = & \displaystyle \frac{1}{2t}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t},\end{eqnarray}$$

(4.12) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}^{2}} & = & \displaystyle \frac{1}{4t^{2}}\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}t^{2}}-\frac{1}{t}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\right).\end{eqnarray}$$

The same is done for $\unicode[STIX]{x1D713}_{2}$ and $s$ . Equations (2.5) and (2.6) can be written as

(4.13a,b ) $$\begin{eqnarray}r(\unicode[STIX]{x1D719}_{1},Q_{d})=d,\quad r(\unicode[STIX]{x1D719}_{2},Q_{D})=D,\end{eqnarray}$$

respectively. Finally, we satisfy the dynamic boundary condition (2.8) on the interface in both $\unicode[STIX]{x1D6FA}_{1}$ and $\unicode[STIX]{x1D6FA}_{2}$ . For example, consider (2.8) satisfied in $\unicode[STIX]{x1D6FA}_{1}$ . Making use of (4.6), this gives

(4.14) $$\begin{eqnarray}\frac{1}{2}\left(\left[\left(\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}\right)^{2}+r^{2}\left(\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}\right)^{2}\right]-\unicode[STIX]{x1D70C}\left[\left(\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}\right)^{2}+r^{2}\left(\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{2}}\right)^{2}\right]\right)+\unicode[STIX]{x1D705}_{1}-\frac{B}{2r^{2}}=C,\end{eqnarray}$$

where $\unicode[STIX]{x1D705}_{1}$ is computed using (4.7). Note that the time-dependant term is removed due to the moving frame of reference. We see that we require $\unicode[STIX]{x2202}r/\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}$ and $\unicode[STIX]{x2202}r/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{2}$ as functions of $\unicode[STIX]{x1D719}_{1}$ on the interface to solve this equation in $\unicode[STIX]{x1D6FA}_{1}$ . Similarly, we require $\unicode[STIX]{x2202}r/\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}$ and $\unicode[STIX]{x2202}r/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}$ as functions of $\unicode[STIX]{x1D719}_{2}$ to solve it in $\unicode[STIX]{x1D6FA}_{2}$ . This is done by integrating the identities

(4.15) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}x}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{i}}=r\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}},\quad i=1,2,\end{eqnarray}$$

on the interface to find $x$ as a function of $\unicode[STIX]{x1D719}_{1}$ in $\unicode[STIX]{x1D6FA}_{1}$ , and $x$ as a function of $\unicode[STIX]{x1D719}_{2}$ in $\unicode[STIX]{x1D6FA}_{2}$ . We then interpolate in $x$ to find $\unicode[STIX]{x1D719}_{2}$ as a function of $\unicode[STIX]{x1D719}_{1}$ , since the interface is the same in either domain. An unfortunate consequence of the interpolation procedure is that it requires the interface $\unicode[STIX]{x1D702}$ to be a single valued function of $x$ , meaning the method will not work for overhanging waves.

Fixing a value of $B$ , the system above provides $M(P+N)$ equations for $M(P+N)+4$ unknowns ( $r$ at each mesh point, $C,c,Q$ and $Q_{D}$ ). We obtain three additional equations by fixing the amplitude $A$ of the wave,

(4.16) $$\begin{eqnarray}A=r(0,Q)-r(-c\unicode[STIX]{x1D706}/2,Q),\end{eqnarray}$$

and the wavelength $\unicode[STIX]{x1D706}$ ,

(4.17a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D706}=\left[\int _{-c\unicode[STIX]{x1D706}}^{0}rr_{\unicode[STIX]{x1D713}_{1}}\,\text{d}\unicode[STIX]{x1D719}_{1}\right]_{\unicode[STIX]{x1D713}_{1}=Q},\quad \unicode[STIX]{x1D706}=\left[\int _{-c\unicode[STIX]{x1D706}}^{0}rr_{\unicode[STIX]{x1D713}_{2}}\,\text{d}\unicode[STIX]{x1D719}_{2}\right]_{\unicode[STIX]{x1D713}_{2}=Q}.\end{eqnarray}$$

Finally, we fix the mean displacement of the interface ( $R=1$ ) by writing

(4.18) $$\begin{eqnarray}\left[\int _{-\unicode[STIX]{x1D706}c/2}^{0}(r-1)rr_{\unicode[STIX]{x1D713}_{1}}\,\text{d}\unicode[STIX]{x1D719}_{1}\right]_{\unicode[STIX]{x1D713}_{1}=Q}=0.\end{eqnarray}$$

In some instances, it is convenient to fix instead the speed $c$ and allow the amplitude $A$ to be an unknown. The discrete system of $M(P+N)+4$ equations for $M(P+N)+4$ unknowns can be solved numerically via Newton’s method. We terminate the iterations in Newton’s method once the $L^{\infty }$ -norm of the residuals is of order $10^{-11}$ .

When considering pure solitary waves, the far-field condition (2.11) is equivalent to demanding  $r$ tends to the uniform stream solution (4.2) as $\unicode[STIX]{x1D719}_{i}\rightarrow \pm \infty$ . Furthermore, the far-field condition fixes the Bernoulli constant $C=(1-\unicode[STIX]{x1D70C})c^{2}/2+1-B/2$ (see (2.8)) and the fluxes $Q=c/2$ and $Q_{D}=cD^{2}/2$ . In such circumstances, we replace the governing equation (4.5) with (4.2) at the mesh points $\unicode[STIX]{x1D719}_{1_{1}}=\unicode[STIX]{x1D719}_{2_{1}}=-c\unicode[STIX]{x1D706}/2$ . Again, we obtain $M(N+P)$ equations from the field equation and boundary conditions. We obtain an additional equation by fixing the amplitude of the wave, which for solitary waves we choose to be the value of $r$ on the interface at the point of symmetry. This results in $M(P+N)+1$ equations for the $M(P+N)+1$ unknowns ( $r$ at each mesh point, and $c$ ). We must take $\unicode[STIX]{x1D706}$ large enough such that the solution becomes identical within graphical accuracy to further increase in $\unicode[STIX]{x1D706}$ . This is common practise when computing solitary waves (for example, see Byatt-Smith & Longuet-Higgins (Reference Byatt-Smith and Longuet-Higgins1976)), since computationally we cannot solve for infinitely large domains.

The numerical scheme described above is used to find solutions for the two-layer model. When finding solutions for the one-layer problem, we do not need to solve for values of $r$ in the domain $\unicode[STIX]{x1D6FA}_{2}$ , or the value $Q_{D}$ . For example, for one-layer periodic waves, there are $MN+3$ unknowns ( $r$ at each mesh point in $\unicode[STIX]{x1D6FA}_{1}$ , $C,c$ , and $Q$ ). We solve the field equation (4.5) at interior nodes of $\unicode[STIX]{x1D6FA}_{1}$ . Furthermore, we satisfy (4.14) with $\unicode[STIX]{x1D70C}=0$ on $\unicode[STIX]{x1D713}_{1}=Q$ , as well as (4.13a ), (4.16), (4.17a ) and (4.18). This results in a closed discrete system of $MN+3$ equations for $MN+3$ unknowns. Furthermore, since we do not require values from $\unicode[STIX]{x1D6FA}_{2}$ to solve (4.14) in $\unicode[STIX]{x1D6FA}_{1}$ , we no longer need to interpolate values in $x$ , as is done in the two-layer problem. This allows us to compute overhanging solutions for the one-layer model.

Typical mesh sizes for periodic waves are $M=200$ , and $N$ and $P$ are chosen such that differences in $t$ are approximately equal to differences in $s$ . For example, with $d=0.5$ and $D=2$ , we took $N=30$ and $P=60$ . For solitary waves, larger values of $M$ are considered. Meshes of this size are possible due to the sparsity of the Jacobian matrix. Furthermore, for more extreme profiles, it can be useful to perform the coordinate transforms $\unicode[STIX]{x1D719}=-c\unicode[STIX]{x1D706}(1-\unicode[STIX]{x1D6FC}^{2})/2$ or $\unicode[STIX]{x1D719}=-c\unicode[STIX]{x1D706}\unicode[STIX]{x1D6FC}^{2}/2$ on either $\unicode[STIX]{x1D719}_{1}$ or $\unicode[STIX]{x1D719}_{2}$ (or both), and then take equally spaced points in $\unicode[STIX]{x1D6FC}\in [0,1]$ . The first transformation condenses points close to $\unicode[STIX]{x1D719}=-c\unicode[STIX]{x1D706}/2$ , while the second condenses points near $\unicode[STIX]{x1D719}=0$ . The transformation is chosen such that the distribution of points is more uniform. There are less points in areas of small velocities if equally spaced points in $\unicode[STIX]{x1D719}$ are used.

In the following section, we discuss the possible static configurations of the problem.

5 Static profiles

It is helpful to discuss static configurations of this problem ( $c=0$ ), since many of the dynamic solution branches terminate on static profiles. Setting all time derivatives and velocities to zero in (2.8), it is left to find $\unicode[STIX]{x1D702}$ that satisfies

(5.1) $$\begin{eqnarray}\unicode[STIX]{x1D705}-\frac{B}{2\unicode[STIX]{x1D702}^{2}}=C,\end{eqnarray}$$

where $\unicode[STIX]{x1D705}$ is the mean curvature. BP solved (5.1) by parameterising the problem in terms of arclength $s$ , and expressing it as a two-dimensional system for the unknowns $\unicode[STIX]{x1D702}$ and $\unicode[STIX]{x1D6FC}$ , where $\unicode[STIX]{x1D6FC}=\tan \unicode[STIX]{x1D702}_{x}$ . We repeat their findings below for the sake of completeness. They found that the energy, $E$ , given by

(5.2) $$\begin{eqnarray}E=\unicode[STIX]{x1D702}\cos \unicode[STIX]{x1D6FC}-\frac{C}{2}\unicode[STIX]{x1D702}^{2}-\frac{B}{2}\log \unicode[STIX]{x1D702},\end{eqnarray}$$

is a conserved quantity. Curves of constant $E$ correspond to trajectories in the $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D702})$ plane. Full details can be found in § 4 of BP. There are four possible fixed points of the system, given by

(5.3a-d ) $$\begin{eqnarray}(2n\unicode[STIX]{x03C0},\unicode[STIX]{x1D6FD}_{+}),\quad (2n\unicode[STIX]{x03C0},\unicode[STIX]{x1D6FD}_{-}),\quad ((2n+1)\unicode[STIX]{x03C0},\unicode[STIX]{x1D6FE}_{+}),\quad ((2n+1)\unicode[STIX]{x03C0},\unicode[STIX]{x1D6FE}_{-}),\end{eqnarray}$$

where

(5.4a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}_{\pm }=\frac{1\pm \sqrt{1-2CB}}{2C},\quad \unicode[STIX]{x1D6FE}_{\pm }=\frac{-1\pm \sqrt{1-2CB}}{2C}.\end{eqnarray}$$

Since we only consider solutions with $\unicode[STIX]{x1D702}\geqslant 0$ , assuming $B>0$ , the existence of these fixed points can be broken down into three cases. In the first case, when $C<0$ , we find that the fixed point $(2n\unicode[STIX]{x03C0},\unicode[STIX]{x1D6FD}_{-})$ is a saddle point, and the fixed point $((2n+1)\unicode[STIX]{x03C0},\unicode[STIX]{x1D6FE}_{-})$ is a centre. The other two fixed points are unphysical, and are ignored. Figure 4(a) shows trajectories in the $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D702})$ space. The heteroclinic orbits (solid lines) connecting the saddle points at $(2n\unicode[STIX]{x03C0},\unicode[STIX]{x1D6FD}_{-})$ and $(2(n+1)\unicode[STIX]{x03C0},\unicode[STIX]{x1D6FD}_{-})$ correspond to solitary waves with radial displacement $\unicode[STIX]{x1D6FD}_{-}$ in the far field. As expected, when $C=1-B/2$ , this value is unity. These solutions self-intersect, and are hence unphysical. The circular orbits (dotted lines) contained inside the heteroclinic orbits correspond to smooth periodic profiles, while the $2\unicode[STIX]{x03C0}$ periodic curves (dashed lines) correspond to self-intersecting periodic profiles.

Figure 4.

Curves of constant $E$ in the $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D702})$ plane. Both panels are for $B=1.25$ . (a) Has $C=-1$ , while (b) has $C=0.2$ . The critical points are labelled with crosses. The solid lines correspond to heteroclinic and homoclinic orbits, the dotted lines circular orbits and dashed lines $2\unicode[STIX]{x03C0}$ periodic curves.

Next, when $0<C<1/2B$ , we find that the fixed point $(2n\unicode[STIX]{x03C0},\unicode[STIX]{x1D6FD}_{-})$ is again a saddle point, and the fixed point $(2n\unicode[STIX]{x03C0},\unicode[STIX]{x1D6FD}_{+})$ is a centre. Figure 4(b) is an example of the $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D702})$ space. The homoclinic orbit connecting the saddle point to itself corresponds to a smooth elevation solitary wave profile. The heteroclinic orbits connecting the saddle points are again self-intersecting solitary waves. Circular orbits correspond to smooth periodic profiles, and $2\unicode[STIX]{x03C0}$ periodic curves are self-intersecting periodic solutions. Finally, when $C>1/(2B)$ , there are no physical fixed points, and all trajectories are $2\unicode[STIX]{x03C0}$ periodic curves in the $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D702})$ space, corresponding to self intersecting periodic profiles.

All of the solutions described above are one-dimensional profiles that satisfy (5.1) and $\unicode[STIX]{x1D702}\geqslant 0$ with $c=0$ . We integrate for values of $x$ along the curve of constant $E$ via the integral

(5.5) $$\begin{eqnarray}x=\int _{E=\text{const.}}\cot \unicode[STIX]{x1D6FC}\,\text{d}\unicode[STIX]{x1D702},\end{eqnarray}$$

to obtain the profile in the $(x,\unicode[STIX]{x1D702})$ space. This integral is evaluated numerically using the trapezoidal rule.

The solutions do not take into consideration any boundaries at $r=d$ and $r=D$ . It is of interest to note that we can interpret the profiles even if they intersect a boundary: the solutions can be seen as profiles which touch a boundary. We then consider the profile up to the point of contact, where the solution is reflected. This is demonstrated in figure 5, where two examples of a static profile (dashed curves) crossing a boundary (dotted curve) from above and below (a and b respectively) are interpreted this way. We only consider the portion of the profile satisfying $d\leqslant \unicode[STIX]{x1D702}\leqslant D$ . In (a), the dashed profile self-intersects, and is hence not physical without the inclusion of a boundary. The dashed profile of (b) is a static elevation solitary wave, and is a valid solution without the boundary. We note that these modified solutions disregard the complicated physical properties of contact angles (for example, see Batchelor Reference Batchelor1994, § 1.9). Despite this, the solutions are still of importance to consider, since many dynamic solution branches approach such static limiting configurations, as shown in § 6.

Figure 5.

The dashed curves are profiles of static configurations. In (a), this static solution corresponds to a $2\unicode[STIX]{x03C0}$ periodic curve in the $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D702})$ space, while in (b), it corresponds to a homoclinic orbit. The dotted curves we take as boundaries below (a) or above (b) the profile. The black curves show the modified solution, taken by reflecting the relevant part of the dashed profile.

In the next section, we discuss the results of the numerical procedure discussed in § 4 for non-static solutions, and how they relate to the static solutions discussed above.