3.1. Formulation of the problem

We consider a column of ferrofluid, fluid 1, with magnetic susceptibility $\chi =\chi _1$, centred on a rigid wire with radius $a$. We choose the cylindrical system $(r,\theta,z)$ such that $r$ and $\theta$ are the radial and azimuthal coordinates, and $z$ points along the wire. In the stationary state, fluid 1 is in the region $a< r< R$ and is surrounded by another ferrofluid, fluid 2, whose domain is unbounded, with susceptibility $\chi =\chi _2$. Both fluids are incompressible and isothermal, with constant density $\rho$, viscosity $\eta$, and permeability $\mu _0$. A steady electric current ${\boldsymbol {J}}=J_0{\boldsymbol {e}_z}$ runs through the wire, producing an azimuthal magnetic field ${\boldsymbol {H}}=J_0/2{\rm \pi} r\boldsymbol {e}_{\theta }$, where $\boldsymbol {e}_z$ is the unit vector in the $z$ direction, and $\boldsymbol {e}_{\theta }$ is the unit vector in the anticlockwise, azimuthal direction. The set-up is shown in figure 1.

Figure 1.

Schematic of the two-fluid system.

Take $R$ as the length scale, and define $a=a_*R$. We non-dimensionalise pressure as $p=\sigma p_*/R$, where $\sigma$ is the surface tension, and non-dimensionalise the field as $\boldsymbol {H}=J_0\boldsymbol {H_*}/2{\rm \pi} R$, and pick the time scale, $T$, to be $T=\eta R/\sigma$, such that the velocity is non-dimensionalised as $\boldsymbol {u}=\sigma /\eta \boldsymbol {u}_*$. The starred variables are dimensionless, but we drop the stars from here on. Equation  (2.10) becomes

(3.1)\begin{equation} {{Re}}\,\frac{{\rm D}\boldsymbol{u}}{{\rm D}t}+\boldsymbol{\nabla}p={\nabla}^2\boldsymbol{u}+B\boldsymbol{f}, \end{equation}

where

(3.2)\begin{equation} {{Re}}=\frac{\rho \sigma R}{\eta^2} \end{equation}

is the Reynolds number, and

(3.3)\begin{equation} B=\frac{\mu_0J_0^2}{4{\rm \pi}^2 R\sigma} \end{equation}

is the magnetic Bond number. Since

(3.4)\begin{equation} \left. \begin{array}{ll@{}} \displaystyle \chi=0, & {\rm for}\ r \leq a,\\ \displaystyle \chi=\chi_1, & {\rm for}\ a< r < R,\\ \displaystyle \chi=\chi_2, & {\rm for}\ r \geq R, \end{array}\right\} \end{equation}

$\chi$ is always constant, resulting in $\boldsymbol {f}=\boldsymbol {0}$ in (3.1), and magnetic effects are felt only at the interface.

Initially, both fluids are at rest and the interface between the two fluids is located at $r=1$. We consider perturbations to the surface such that the surface is located at

(3.5)\begin{equation} r=1+\epsilon\,{\rm Re}(\hat{S}\zeta), \end{equation}

where $\zeta =\exp \left ({\textrm {i}(kz+m\theta )+st}\right )$, $\epsilon \ll 1$, $k,m$ are real and positive wavenumbers, $\hat {S}$ may be a complex constant (or it could be unity), and $s$ is the growth rate of the disturbance and could be complex. In (3.5), the real part of the perturbation is taken, and this is done from here on for the other variables, but it is not written explicitly. The normal vector to the surface becomes

(3.6)\begin{equation} \boldsymbol{{n}}=\left(1,-\frac{\epsilon{\rm i}m\hat{S}\zeta}{r},-\epsilon{\rm i}k\hat{S}\zeta\right)^{\rm T}, \end{equation}

and the tangential vectors are

(3.7a,b)\begin{equation} {\boldsymbol{{\tau}_1}}=(\epsilon{\rm i}k\hat{S}\zeta,0,1)^{\rm T}, \quad \boldsymbol{{\tau}_2}=\left(\frac{\epsilon{\rm i}m\hat{S}\zeta}{r},1,0\right)^{\rm T}. \end{equation}

The perturbed pressure $p^{(i)}$ and velocity $\boldsymbol {u}^{(i)}$, where $i=1,2$ for fluids $1$ and $2$, are

(3.8a,b)\begin{equation} p^{(i)}=p_0+\epsilon\,\hat{p}^{(i)}(r)\,\zeta \quad {\rm and} \quad \boldsymbol{u}^{(i)}=\epsilon\,\hat{\boldsymbol{u}}^{(i)}(r)\,\zeta, \end{equation}

where $p_0$ is constant. The perturbations satisfy

(3.9a,b)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \hat{\boldsymbol{u}}=0 \quad {\rm and} \quad s\,{{Re}}\,\hat{\boldsymbol{u}}+\boldsymbol{\nabla}\hat{p}={\nabla}^2\hat{\boldsymbol{u}}. \end{equation}

In component form we have

(3.10)\begin{gather} (r\hat{u})'+{\rm i}m\hat{v}+{\rm i}kr\hat{w}=0, \end{gather}

(3.11)\begin{gather} r^2\hat{p}'={-}2{\rm i}m\hat{v}-\hat{u}+r^2\hat{u}''+r\hat{u}'-(m^2+\bar{k}^2r^2)\hat{u}, \end{gather}

(3.12)\begin{gather} {\rm i}mr\hat{p}={-}\hat{v}+2{\rm i}m\hat{u}+r^2\hat{v}''+r\hat{v}'-(m^2+\bar{k}^2r^2)\hat{v}, \end{gather}

(3.13)\begin{gather} {\rm i}kr^2\hat{p}=r^2\hat{w''}+rw'-(m^2+\bar{k}^2r^2)\hat{w}, \end{gather}

where $'$ denotes the first derivative with respect to $r$, and $\bar {k}=\sqrt {k^2+s\,{{Re}}}$. The general solution is in terms of the modified Bessel functions of the first and second kind, $\textrm {I}_n(z)$ and $\textrm {K}_n(z)$, respectively, where we write $\textrm {I}_n, \textrm {K}_n$ when $z=kr$, and $\bar {\textrm {I}}_n, \bar {\textrm {K}}_n$ when $z=\bar {k}r$, but give the argument otherwise. The general solution of (3.9a,b) is modified from Saville (Reference Saville1971) and Mestel (Reference Mestel1996) to account for the inner wire at $r=a$, and is found to be

(3.14a)\begin{align} \hat{p}^{(i)} & = c_1^{(i)}{\rm I}_m + c_2^{(i)}{\rm K}_{m}, \end{align}

(3.14b)\begin{align} \hat{u}^{(i)} & ={-}\frac {1}{(s\,{{Re}})^{2}r} \left(c_1^{(i)}(kr{\rm I}_{m+1}+m{\rm I}_m)+c_2^{(i)}(m{\rm K}_m-kr{\rm K}_{m+1})\right)\nonumber\\ &\quad -\frac{{\rm i}k}{\bar{k}}\left(c_3^{(i)}\bar{{\rm I}}_{m+1}+c_4^{(i)}\bar{{\rm K}}_{m+1}\right)+\frac{2m}{\bar{k}r}\left(c_5^{(i)}\bar{{\rm I}}_m+c_6^{(i)}\bar{{\rm K}}_{m}\right),\end{align}

(3.14c)\begin{align} \hat{v}^{(i)}& = \left( {{\frac {-c_3^{(i)}k}{\bar{k}}}}+2{\rm i}c_5^{(i)} \right) \bar{{\rm I}}_{m+1}+{\frac {2{\rm i}m}{\bar{k}r}}\left(c_5^{(i)}\bar{{\rm I}}_m+c_6^{(i)}\bar{{\rm K}}_{m}\right)-\left({{\frac {c_4^{(i)}k}{\bar{k}}}}+2{\rm i}c_6^{(i)} \right) \bar{{\rm K}}_{m+1}\nonumber\\ &\quad -\frac{{\rm i}m}{{(s\,{{Re}})}^{2}r}\left(c_2^{(i)}{\rm K}_{m}+c_1^{(i)}{\rm I}_m\right),\end{align}

(3.14b)\begin{align} \hat{w}^{(i)}& = {\frac {-{\rm i}k}{{(s\,{{Re}})}^{2}}}\left(c_1^{(i)}{\rm I}_m+c_2^{(i)}{\rm K}_{m}\right)+c_3^{(i)}\bar{{\rm I}}_m-c_4^{(i)}\bar{{\rm K}}_{m}, \end{align}

for constants $c_1^{(i)},\ldots, c_6^{(i)}$. To satisfy $u^{(2)} \rightarrow 0$ as $r \rightarrow \infty$, $\mathrm {Re}(\bar {k})>0$ and $c_1^{(2)}=c_3^{(2)}=c_5^{(2)}=0$.

We perturb the magnetic potential such that

(3.15)\begin{equation} \phi^{(l)}=\theta+\epsilon\,\hat{\phi}^{(l)}(r)\,\zeta, \end{equation}

where $l=0,1,2$ for the wire, inner fluid and outer fluid, respectively. Equation (2.4) gives

(3.16)\begin{equation} r^2\hat{\phi''}^{(l)}+\hat{\phi'}^{(l)}-(m^2+k^2r^2)\hat{\phi}^{(l)}=0, \end{equation}

with general solution

(3.17)\begin{equation} \hat{\phi}(r)^{(l)}=q_1^{(l)}{\rm I}_m+q_2^{(l)}{\rm K}_m, \end{equation}

for constants $q_1^{(l)},q_2^{(l)}$. For $\hat {\phi }^{(0)}$ regular at $r=0$,

(3.18)\begin{equation} \hat{\phi}^{(0)}=q_1^{(0)}{\rm I}_m, \end{equation}

and imposing $\phi ^{(2)}\rightarrow 0$, as $r\rightarrow \infty$ gives

(3.19)\begin{equation} \hat{\phi}^{(2)}=q_2^{(2)}{\rm K}_m. \end{equation}

Equations (2.5) and (2.6) give

(3.20a,b)\begin{equation} (1+\chi_2)\hat{\phi}'^{(2)}-(1+\chi_1)\hat{\phi}'^{(1)}+im(\chi_1-\chi_2)\hat{S}=0 , \quad \hat{\phi}^{(1)}=\hat{\phi}^{(2)} \end{equation}

at $r=1$, and

(3.21a,b)\begin{equation} (1+\chi_1)\hat{\phi}'^{(1)}-\hat{\phi}'^{(0)}=0, \quad \hat{\phi}^{(1)}=\hat{\phi}^{(0)} \end{equation}

at $r=a$, determining the constants $q_1^{(0)},q_1^{(1)},q_2^{(1)},q_2^{(2)}$ given in (A1)–(A4) in Appendix A.

At $r=a$,

(3.22)\begin{equation} \hat{\boldsymbol{u}}^{(1)}=0, \end{equation}

and at $r=1$,

(3.23)\begin{equation} \hat{\boldsymbol{u}}^{(1)}=\hat{\boldsymbol{u}}^{(2)}. \end{equation}

At the interface of the fluids, there is a normal stress condition

(3.24) \begin{equation} [\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{T}\boldsymbol{\cdot} \boldsymbol{n}]=\sigma\, \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{n}, \end{equation}

and two tangential stress conditions

(3.25)\begin{gather} {[}\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{T}\boldsymbol{\cdot} \boldsymbol{\tau_1}]=0, \end{gather}

(3.26)\begin{gather} {[}\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{T}\boldsymbol{\cdot} \boldsymbol{\tau_2}]=0. \end{gather}

Non-dimensionalising (3.24), substituting the perturbed variables, linearising, and invoking continuity of the normal component of $\boldsymbol {B}$, we obtain

(3.27) \begin{align} (m^2+k^2-1+B(\chi_1-\chi_2))\hat{S}={\rm i}mB(\chi_1\hat{\phi}^{(1)} -\chi_2\hat{\phi}^{(2)})+\hat{p}^{(1)}- \hat{p}^{(2)}-2\hat{u}'^{(1)}+2\hat{u}'^{(2)} \end{align}

at $r=1$. Similarly, (3.25) and (3.26) give

(3.28)\begin{equation} \hat{v}'^{(1)}-\hat{v}'^{(2)}=0 \end{equation}

and

(3.29)\begin{equation} \hat{w}'^{(1)}-\hat{w}'^{(2)}=0 \end{equation}

at $r=1$. Equations (3.22), (3.23) (3.27)–(3.29) determine the constants $c_1^{(i)},\ldots, c_6^{(i)}$ given in (A13)–(A21).

The growth rate appears in the kinematic condition at $r=1$:

(3.30)\begin{equation} s\hat{S}-\hat{u}^{(i)}=0. \end{equation}

Consequently, substituting $\hat {u}^{(i)}$ into (3.30), we obtain

(3.31) \begin{equation} s={-}gf(f_1(k^2+m^2-1+{B}(\chi_1-\chi_2))+f_2{B}(\chi_1-\chi_2)^2m^2), \end{equation}

where $g, f_1, f_2>0$. Here, $g, f_1, f_2$ are functions of $m,k,a,\chi _1,\chi _2$, and $f$ is a function of $\bar {k}, m, k, a,\chi _1,\chi _2$, all given in (A10)–(A12). Since $\bar {k}$ is a function of $s$, $f$ is a function of $s$, and therefore (3.31) is an implicit relation that must be solved numerically.

3.2. Highly viscous and inviscid limits

In the limit ${{Re}} \rightarrow 0$, $f\rightarrow f_v$, where $f_v>0$ and $f_v$ is no longer a function of $s$, given in (A22). The growth rate in the highly viscous regime, $s_v$, is expressed explicitly as

(3.32)\begin{equation} s_v={-}gf_v(f_1(k^2+m^2-1+{B}(\chi_1-\chi_2))+f_2{B}(\chi_1-\chi_2)^2m^2). \end{equation}

In the inviscid limit, $\eta \rightarrow 0$, and a more appropriate scaling for time, $T_I$, is $T_I=\sqrt {R^3\rho /\sigma }$. Since $T_I=\sqrt {{{Re}}}\,T$, we substitute $s=s_I/\sqrt {{{Re}}}$ into (3.31), where $s_I$ is the inviscid growth rate. Taking the limit as ${{Re}} \rightarrow \infty$ gives $f\rightarrow f_I/s_I$, where $f_I>0$ and $f_I$ is no longer a function of $s$, given in (A23). Consequently,

(3.33)\begin{equation} s_I^2={-}gf_I(f_1(k^2+m^2-1+{B}(\chi_1-\chi_2))+f_2{B}(\chi_1-\chi_2)^2m^2). \end{equation}

More simply, (3.33) can be obtained by taking the limit $\eta \rightarrow 0$ in the governing equations from the outset, and applying the boundary conditions for an inviscid system ($\eta =0$), namely $\hat {u}^{(1)}=0$ at the wire, and $\hat {u}^{(1)}=\hat {u}^{(2)}$, (3.27), (3.29) at the interface.

Since $f_I,f_v>0$, the system is stable (or neutrally stable), in both the inviscid and highly viscous regimes, if and only if

(3.34)\begin{equation} f_1(k^2+m^2-1+{B}(\chi_1-\chi_2))+f_2{B}(\chi_1-\chi_2)^2m^2\geq 0. \end{equation}

Note that in the inviscid regime, if (3.34) holds, then the system is neutrally stable as the growth rate is imaginary. If the inner fluid has a higher susceptibility than the outer fluid, then the system can be unstable only as a result of capillary forces. Only axisymmetric modes can be unstable when $k<1$, and increasing the current in the wire will stabilise the system provided that $B(\chi _1-\chi _2)>1$. Figures 2 and 3 show the growth rate of the modes being dampened as the current in the wire is increased for the viscous and inviscid regimes, respectively.

Figure 2.

Viscous growth rate, $a=0.1$, $\chi _1=5$, $\chi _2=1$: (a) $B=0$, (b) $B=0.1$, (c) $B=0.25$, (d) $B=4$.

Figure 3.

Inviscid growth rate, $a=0.1$, $\chi _1=5$, $\chi _2=1$: (a) $B=0$, (b) $B=0.1$, (c) $B=0.25$, (d) $B=4$.

On the other hand, when the outer fluid has a higher susceptibility, capillary and magnetic forces may be destabilising. Increasing the current, thereby increasing the magnetic forcing at the interface, renders non-axisymmetric modes unstable, as well as axisymmetric modes, and for sufficiently large $B$, all modes $m$ can be rendered unstable. Figures 4–6 show that increasing $B$ results in an increase in unstable modes, and increases the magnitude of their growth rates. Moreover, in the highly viscous regime, figure 5 shows $m=1$, $k\rightarrow 0$, is the most unstable mode for $a=0.1$, and in fact $s\rightarrow \infty$ as $k\rightarrow 0$, whereas when $a=0.5$ in figure 6, $m=0$ remains the most unstable mode. Performing a series expansion on the viscous growth rate as $a,k\rightarrow 0$, we find $s\sim \ln (ka)$ when $m=1$, thus $s\rightarrow \infty$ in the limit, but $s$ converges to a constant when $m=0$ or $m>1$, a result seen in the context of electro-hydrodynamics too (Saville Reference Saville1971; Mestel Reference Mestel1996). It is important to note that when comparing the inviscid regime with the viscous regime, $s_I$ and $s_v$ are on different time scales.

Figure 4.

Inviscid growth rate, $a=0.1$, $\chi _1=1$, $\chi _2=5$: (a) $B=0$, (b) $B=1$, (c) $B=10$.

Figure 5.

Viscous growth rate, $a=0.1$, $\chi _1=1$, $\chi _2=5$: (a) $B=0$, (b) $B=1$, (c) $B=10$.

Figure 6.

Viscous growth rate, $a=0.5$, $\chi _1=1$, $\chi _2=5$: (a) $B=0$, (b) $B=1$, (c) $B=10$.

3.3. Arbitrary Reynolds number

A numerical root solver in the program Maple is used on (3.31) for specific values of $k,a,m,\chi _1,\chi _2,B,{{Re}}$, to find the associated growth rate of the mode. The stability condition (3.34) appears to hold for all Reynolds numbers. Figure 7 is the growth rate plotted when $a=0.1$, $k=0.5$, $\chi _1=1$, $\chi _2=5$, $B=0.1$, for a range of Reynolds numbers, showing two stable branches when $m=1$ and an unstable branch for $m=0$. Given a stable mode, as ${{Re}} \rightarrow 0$, there are two branches, both real: one branch tends to $s_v$, and the other tends to $-\infty$, the latter a result of ${{Re}} \rightarrow 0$ for the chosen time scale. As ${{Re}}$ is increased, the two branches meet and then split, becoming complex conjugates of each other, tending towards $\pm s_I$, where $s_I$ is purely imaginary for a stable mode. Given an unstable mode, we get one branch, starting at $s_v$ and tending to $|s_I|$. There exists a branch that tends to the negative inviscid root, $-|s_I|$, but this is invalid for finite ${{Re}}$, since the boundary conditions require $\mathrm {Re}({\bar {k}})\geq 0$. When $\chi _2>\chi _1$, and $a,k$ are sufficiently small, $m=1$ is more unstable than $m=0$. Yet for all $k$, $m=1$ is the most unstable mode only for sufficiently small ${{Re}}$. This is shown in figure 8, where the growth rate is plotted against $k$ for different Reynolds numbers, showing that when ${{Re}}=0.001,0.1$, $m=1$ is the most unstable mode, but for the other Reynolds numbers shown, it is not. We find that increasing the current in the wire stabilises the system if $\chi _1>\chi _2$ for all Reynolds numbers. Figures 7 and 9, where $\chi _2>\chi _1$, show that increasing $B$ from $B=0.1$ to $B=0.5$ renders the mode $m=1$ unstable, and we find that for all Reynolds numbers, increasing the current does not stabilise the system if $\chi _2>\chi _1$, but renders more modes unstable.

Figure 7.

Growth rate plotted when $a=0.1$, $k=0.5$, $\chi _1=1$, $\chi _2=5$ and $B=0.1$ for arbitrary ${{Re}}$. Panels (a) and (b) plot, respectively, the real and imaginary parts of $s$. The $m=0$ branch is purely real, but the $m=1$ branch starts as real for low ${{Re}}$ and then becomes a complex conjugate pair.

Figure 8.

Growth rate plotted against $k$, for arbitrary ${{Re}}$, when $a=0.1$, $\chi _1=1$, $\chi _2=5$, $B=10$: (a) ${{Re}}=0.0001$, (b) ${{Re}}=0.1$, (c) ${{Re}}=0.5$, (d) ${{Re}}=1$, (e) ${{Re}}=10$, ( f) ${{Re}}=100$.

Figure 9.

Growth rate plotted when $a=0.1$, $k=0.5$, $\chi _1=1$, $\chi _2=5$ and $B=0.5$ for arbitrary ${{Re}}$.

3.4. Stabilisation with an axial field

To stabilise the system, irrespective of whether $\chi _2>\chi _1$ or $\chi _1>\chi _2$, we consider $\boldsymbol {H}_0=(0,1/r,Z)^\textrm {T}$, $Z$ constant, thereby adding an axial field. It follows that

(3.35)\begin{equation} \phi_0=\theta +Zz, \end{equation}

and we perform analysis analogous to that in § 3.1. The general solutions (3.14), (3.17) still hold, and $\hat {\phi }^{(0)},\hat {\phi }^{(2)}$ are still given by (3.18) and (3.19), respectively. At $r=a$, we apply (3.21) and (3.22). At $r=1$, we apply (3.23), (3.28), (3.29), but (2.5), (2.6), (3.24) now give

(3.36a,b)\begin{equation} (1+\chi_2)\hat{\phi}'^{(2)}-(1+\chi_1)\hat{\phi}'^{(1)}+{\rm i}(m+kZ)(\chi_1-\chi_2)\hat{S}=0 , \quad \hat{\phi}^{(1)}=\hat{\phi}^{(2)}, \end{equation}

and

(3.37) \begin{align} (m^2+k^2-1+B(\chi_1-\chi_2))\hat{S}&={\rm i}(m+kZ)B(\chi_1\hat{\phi}^{(1)} -\chi_2\hat{\phi}^{(2)})\nonumber\\&\quad +\hat{p}^{(1)}- \hat{p}^{(2)}-2\hat{u}'^{(1)}+2\hat{u}'^{(2)}. \end{align}

We now obtain

(3.38)\begin{gather} s={-}gf(f_1(k^2+m^2-1+{B}(\chi_1-\chi_2))+f_2{B}(\chi_1-\chi_2)^2(Zk+m)^2), \end{gather}

(3.39)\begin{gather} s_v={-}gf_v(f_1(k^2+m^2-1+{B}(\chi_1-\chi_2))+f_2{B}(\chi_1-\chi_2)^2(Zk+m)^2), \end{gather}

(3.40)\begin{gather} s_I^2={-}gf_I(f_1(k^2+m^2-1+{B}(\chi_1-\chi_2))+f_2{B}(\chi_1-\chi_2)^2(Zk+m)^2), \end{gather}

analogously to (3.31), (3.32) and (3.33). Equations (3.39) and (3.40) show that a sufficiently large $kZ$ will stabilise all modes in the inviscid and highly viscous regimes, irrespective of the sign of $(\chi _1-\chi _2)$, provided that $B\neq 0$, and this result appears to hold for all Reynolds numbers. Although extremely long waves in the $z$ direction, $k\rightarrow 0$, would remain unstable, by physical restrictions of the system, $k$ is bounded from zero.

4.1. Axisymmetric disturbances

For axisymmetric disturbances, (4.11) becomes

(4.18a)$$\begin{gather} (s\,{{Re}}-\mathcal{L}_0)\hat{\omega}_r = 0, \end{gather}$$

(4.18b)$$\begin{gather}(s\,{{Re}}-\mathcal{L}_0)\hat{\omega}_{\theta} ={-}\frac{{\rm i}kB\hat{\chi}}{r^3}, \end{gather}$$

(4.18c)$$\begin{gather}(s\,{{Re}}-\mathcal{D}_0)\hat{\omega}_z =0, \end{gather}$$

where

(4.19)\begin{gather} \mathcal{L}_0=\frac{{\rm d}^2}{{\rm d}r^2}+\frac{1}{r}\,\frac{{\rm d}}{{\rm d}r}-k^2-\frac{1}{r^2}, \end{gather}

(4.20)\begin{gather} \mathcal{D}_0=\frac{{\rm d}^2}{{\rm d}r^2}+\frac{1}{r}\,\frac{{\rm d}}{{\rm d}r}-k^2. \end{gather}

Define a stream function $\varPsi$ such that $\hat {\boldsymbol {u}}=\boldsymbol {\nabla }\times (0,\varPsi /r,0),$ and use the change of variables $\varPsi =r\psi$ to give

(4.21)\begin{equation} \hat{\boldsymbol{\omega}}={-}\hat{\boldsymbol{\theta}}\mathcal{L}_0\hat{\psi}, \quad {\nabla}^2 \hat{\boldsymbol{\omega}}={-}\hat{\boldsymbol{\theta}}\mathcal{L}_0^2\hat{\psi}, \end{equation}

where $\psi =\hat {\psi }(r)\,\textrm {e}^{\textrm {i}kz+st}$. It follows from the boundary conditions for $\boldsymbol {u}$ that $\hat {\psi },\hat {\psi }'=0$ at $r=a$ and as $r \rightarrow \infty$. Equations (4.15) and (4.18) give the eigenvalue equation

(4.22)\begin{equation} (s^2\,{{Re}}\,\mathcal{L}_0-s\mathcal{L}_0^2)\hat{\psi}={-}\frac{k^2{B}\chi_0'\hat{\psi}}{r^3}. \end{equation}

Rather than find the eigenvalues of (4.22) numerically, we prove a stability condition. Multiply (4.22) by $r\hat {\psi }^*$, where $\hat {\psi }^*$ is the complex conjugate of $\hat {\psi }$. Integrate over the domain, use integration by parts, and invoke the boundary conditions, to obtain

(4.23)\begin{align} s^2\,{{Re}}\int_{a}^{\infty} \left(|\hat{\psi}'|^2+\left(\frac{1}{r^{2}}+k^2\right)|\hat{\psi}|^2\right)r\,{\rm d}r+s\int_a^{\infty}|\mathcal{L}_0\hat{\psi}|^2r\,{\rm d}r\,{-}\,k^2{B}\int_a^{\infty}\frac{\chi_0'|\hat{\psi}|^2}{r^{2}}\,{\rm d}r\,{=}\,0, \end{align}

an equation for $s$ of the form $as^2+bs+c=0$, where $a,b,c$ all depend on $\hat {\psi }$ and therefore $s$ too, but $a,b,c$ are real, as well as $a,b>0$, bounded away from zero. Thus if $\chi _0'>0$, then $c<0$ and there exists a root with $s>0$. Note that in this case, $s$ and $\hat {\psi }$ are real, whereas when $\chi _0'<0$, $c>0$, $\mathrm {Re}(s)<0$ and $s,\hat {\psi }$ could be complex. We conclude that if $\chi _0'>0$ everywhere, then there exists an unstable mode, while if $\chi _0' \leq 0$ everywhere, then the axisymmetric modes are stable.

We prove a stronger stability condition using variational methods. Crucially, we have shown that if the flow is unstable, then $s$ must be real and therefore $\hat {\psi }$ is real, thus it suffices to consider only $y$ real in the following argument. Consider the functional

(4.24)\begin{align} I(y)=\frac{-\textstyle\int_a^{\infty}(\mathcal{L}_0 y)^2r\,{\rm d}r \textstyle+\sqrt{\left(\int_a^{\infty}(\mathcal{L}_0 y)^2r\,{\rm d}r\right)^2+ W_0}}{2\,{{Re}}\int_{a}^{\infty}\left((y')^2+\left(\frac{1}{r^{2}}+k^2\right)y^2\right)r\,{\rm d}r},\end{align}

where

(4.25)\begin{align} W_{0}= 4k^2\,{{Re}}\,{B}\int_{a}^{\infty}\left((y')^2+\left(\frac{1}{r^{2}}+k^2\right)y^2\right)r\, {\rm d}r\int_a^{\infty}\frac{\chi_0' y^2}{r^{2}}\,{\rm d}r,\end{align}

for all real functions $y(r)$ satisfying $y,y'=0$ at $r=a$, $r \rightarrow \infty$, and proceed in a manner similar to the Rayleigh–Ritz argument. It can be shown that $I(y)$ is bounded above and has a maximum. Suppose that $y= y_0$ is a stationary point of $I$, and $I( y_0)=I_0$. Consider $y={ y}_0+\epsilon y_1$, where $\epsilon \ll 1$, and $y_1$ satisfies the boundary conditions for $y$. Since $y=y_0$ is a stationary point of $I(y)$, the first variation is zero, and thus after Taylor expanding $I(y_0+\epsilon y_1)$, it follows that

(4.26)\begin{align} -{{Re}}\,I_0^2\int_{a}^{\infty}\left( y_0' y_1'+\left(\frac{1}{r^{2}}+k^2\right) y_0 y_1\right)r\,{\rm d}r- I_0\int_a^{\infty}\mathcal{L}_0 y_0\mathcal{L}_0 y_1r\,{\rm d}r=k^2{B}\int_a^{\infty}\frac{\chi_0' y_0 y_1}{r^{2}}\,{\rm d}r. \end{align}

Invoking the self-adjoint property of $\mathcal {L}_0$, using integration by parts, and the boundary conditions for $y_0, y_1$, we can write (4.26) as

(4.27)\begin{equation} {{Re}}\,I_0^2\int_{a}^{\infty}y_1r\mathcal{L}_0 y_0\,{\rm d}r-I_0\int_a^{\infty}r y_1\mathcal{L}_0^2 y_0\,{\rm d}r=k^2{B}\int_a^{\infty}\frac{\chi_0' y_0 y_1}{r^{2}}\,{\rm d}r. \end{equation}

Equation (4.27) is valid for any $y_1$, and therefore

(4.28)\begin{equation} {{Re}}\,I_0^2\mathcal{L}_0 y_0-I_0\mathcal{L}_0^2 y_0=\frac{k^2{B}\chi_0'y_0}{r^{3}}. \end{equation}

It follows that the stationary points of $I(y)$ satisfy (4.22) with real eigenvalues $s=I_0$, and therefore the stationary points of $I(y)$ correspond to the real eigenvalues of (4.22).

Suppose that

(4.29)\begin{equation} \chi_0'>0, \quad {\rm for}\ r_1\leq r \leq r_2, \end{equation}

and pick an arbitrary real function $\hat {y}(r)$ that satisfies the boundary conditions of $y$ such that

(4.30) \begin{equation} \left. \begin{array}{ll@{}} \displaystyle \hat{y}(r)\neq0, & {\rm for }\ r_1\leq r \leq r_2,\\[3pt] \displaystyle \hat{y}(r)=0, & {\rm for }r\notin [r_1,r_2]. \end{array}\right\} \end{equation}

Substitute $y=\hat {y}(r)$ into (4.24) to give $I(\hat {y})=\eta _M$. It follows that $\eta _M>0$, and either $\eta _M$ is a stationary point of $I(y)$, and therefore a positive, real eigenvalue of (4.22), or there exists a positive stationary point of $I(y)$ greater than $\eta _M$, since $I$ is bounded above. Thus there exists a positive eigenvalue of (4.22), resulting in an unstable mode. We conclude that if and only if $\chi _0'>0$ anywhere in the domain, every axisymmetric mode is unstable.

4.2. Two-dimensional modes

By considering two-dimensional modes $k=0$ in (3.10), and (4.11)–(4.16), we obtain an eigenvalue equation,

(4.31)\begin{equation} (s^2\,{{Re}}\,\mathcal{L}_m-s\mathcal{L}_m^2)\mathcal{M}_m\hat{\phi}={-}m^2{B}\left(\frac{\chi_0'\mathcal{M}_m\hat{\phi}}{r^5}+\frac{m^2\chi_0'\hat{\phi}}{r^3}\right), \end{equation}

where

(4.32)\begin{equation} {\mathcal{L}_m}=\frac{1}{r}\,\frac{{\rm d}}{{\rm d}r}\left(r\,\frac{{\rm d}}{{\rm d}r}\right)-\frac{m^2}{r^2}, \quad \mathcal{M}_m=\frac{r^3}{\chi_0'}\left[(1+\chi_0)\mathcal{L}_m+\chi_0'\,\frac{{\rm d}}{{\rm d}r}\right], \end{equation}

and $\phi,\phi '\rightarrow 0$ as $r\rightarrow \infty$. Taking the limit $k\rightarrow 0$ in (4.17) gives

(4.33)\begin{equation} \hat{\phi'}=\frac{m}{r(1+\chi_0)}\,\hat{\phi}(r) \end{equation}

at $r=a$.

We multiply (4.31) by $r\mathcal {M}_m\hat {\phi }^*$, where $\hat {\phi }^*$ is the complex conjugate of $\hat {\phi }$, and integrate over the domain to obtain

(4.34)\begin{align} \int_a^{\infty}(\mathcal{M}_m\hat{\phi}^*(s^2\,{{Re}}\,\mathcal{L}_m-s\mathcal{L}_m^2)\mathcal{M}_m\hat{\phi})r\,{\rm d}r={-}{B}m^2\int_a^{\infty}\left(\frac{\chi_0'\,|\mathcal{M}_m\hat{\phi}|^2}{r^4}+\frac{m^2\chi_0'}{r^2}\,\hat{\phi}\mathcal{M}_m\hat{\phi}^*\right){\rm d}r. \end{align}

Now,

(4.35)\begin{align} \int_a^{\infty}\frac{1}{r^2}\,\chi_0'\hat{\phi}\mathcal{M}_m\hat{\phi}^*\,{\rm d}r=\int_a^{\infty} \left((1+\chi_0)\hat{\phi}(r(\hat{\phi}^*)')'+r\chi_0'\hat{\phi}(\hat{\phi}^*)'-\frac{m^2(1+\chi_0)|\hat{\phi}|^2}{r}\right){\rm d}r, \end{align}

and integration by parts gives

(4.36)\begin{equation} \int_a^{\infty} \frac{1}{r^2}\,\chi_0'\hat{\phi}\mathcal{M}_m\,\hat{\phi}^*\,{\rm d}r={-}m|\hat{\phi}(a)|^2-\int_a^{\infty} r(1+\chi_0)\left(|\hat{\phi}'|^2+\frac{m^2|\hat{\phi}|^2}{r^2}\right){\rm d}r. \end{equation}

Equation (4.36) and the self-adjoint property of $\mathcal {L}_m$ allow (4.34) to be written as

(4.37)\begin{align} & s^2\,{{Re}}\int_a^{\infty}\left(|(\mathcal{M}_m\hat{\phi})'|^2 +\frac{m^2|\mathcal{M}_m\hat{\phi}|^2}{r^2}\right)r\,{\rm d}r +s\int_a^{\infty}\left|\mathcal{L}_m(\mathcal{M}_m\hat{\phi})\right|^2 r\,{\rm d}r \nonumber\\ &\quad +{B}m^2\left(m^3|\hat{\phi}(a)|^2+\int_a^{\infty}\left[\frac{-\chi_0'\,|\mathcal{M}_m\hat{\phi}|^2}{r^4}+m^2r(1+\chi_0)\left(|\hat{\phi}'|^2+\frac{m^2|\hat{\phi}|^2}{r^2}\right)\right]{\rm d}r\right)=0. \end{align}

It follows that if $\chi _0'<0$ everywhere, then two-dimensional modes are stable. Yet if

(4.38)\begin{equation} \int_a^{\infty}\frac{\chi_0'\,|\mathcal{M}_m\hat{\phi}|^2}{r^4}\,r\,{\rm d}r>m^3|\hat{\phi}(a)|^2+m^2\int_a^{\infty} r(1+\chi_0)\left(|\hat{\phi}'|^2+\frac{m^2|{\hat{\phi}}|^2}{r^2}\right){\rm d}r \end{equation}

holds for an eigenfunction $\hat {\phi }$, then two-dimensional modes are unstable.

Furthermore, by considering the functional

(4.39) \begin{equation} I(y)=\frac{-\int_a^{\infty}\left(\mathcal{L}_m(\mathcal{M}_my)\right)^2 r\,{\rm d}r+\sqrt{(\int_a^{\infty}\left(\mathcal{L}_m(\mathcal{M}_my)\right)^2 r\,{\rm d}r)^2-W_1}}{2\,{{Re}}\int_a^{\infty}\left(\left(\mathcal{M}_my'\right)^2+\dfrac{m^2}{r^2}\left(\mathcal{M}_my\right)^2\right)r\,{\rm d}r}, \end{equation}

where

(4.40) \begin{align} W_1 & = 4m^2B\,{{Re}}\left(\int_a^{\infty}\left((\mathcal{M}_my')^2+\frac{m^2}{r^2}\left(\mathcal{M}_my\right)^2\right)r\,{\rm d}r\right)\left(m^3|\hat{\phi}(a)|^2 \vphantom{\left[\frac{-\chi_0'\,|\mathcal{M}_m\hat{\phi}|^2}{r^4}+m^2r(1+\chi_0)\left(|\hat{\phi}'|^2+\frac{m^2|\hat{\phi}|^2}{r^2}\right)\right]}\right.\nonumber\\ &\quad \left.+\int_a^{\infty}\left[\frac{-\chi_0'\,|\mathcal{M}_m\hat{\phi}|^2}{r^4}+m^2r(1+\chi_0)\left(|\hat{\phi}'|^2+\frac{m^2|\hat{\phi}|^2}{r^2}\right)\right]{\rm d}r\right), \end{align}

for all real functions $y$ satisfying the boundary conditions of $\hat {\phi }$, we prove that the stationary points of $I(y)$ correspond to the real eigenvalues of (4.31) by an argument analogous to that in § 4.1. Again, if (4.29) is true, then by picking a $y=\hat {y}(r)$, where $\hat {y}$ satisfies (4.30), with oscillatory behaviour for $r$ in the interval $[r_1, r_2]$, then $I(\hat {y})>0$ and (4.38) holds as $(\mathcal {M}_m\hat {y})^2\gg \hat {y}'\gg \hat {y}$. Consequently, if and only if $\chi _0'>0$ anywhere in the fluid, every mode where $k=0$ is unstable.

4.3. Highly viscous and inviscid limits

Although we have proven that if $\chi _0'>0$ anywhere, all axisymmetric or two-dimensional disturbances are unstable, we are unable to prove, for arbitrary Reynolds number, that if $\chi _0'<0$ everywhere, then all modes are stable; we have proved only that axisymmetric or two-dimensional disturbances are stable. A global energy argument could be applied here by posing an argument analogous to the Rayleigh's stability argument for centrifugal instability. One considers the change in energy when two parcels of ferrofluid at different radii are interchanged while conserving $\chi (r)$, where the resulting condition for stability is that $\boldsymbol {\nabla } \chi$ must decrease continuously radially. Yet this argument would not account for viscous forces or three-dimensional disturbances.

By considering the inviscid limit of (4.11)–(4.17) and following an argument analogous to that in § 4.1, we prove that for all $m,k$, if and only if $\chi _0'>0$ anywhere in the fluid, the system is unstable, and it can be proved that ${s}^2$ and $\hat {\phi }$ are real. In the highly viscous limit of (4.11)–(4.17), the eigenvalues are also proven to be real for axisymmetric disturbances and modes $k=0$, and these modes are shown to be unstable if and only if $\chi _0'>0$ anywhere in the fluid.

4.4. Stabilisation with an axial field

We now show that by adding an axial field

(4.41)\begin{equation} \boldsymbol{H}_0=\left(0,\frac{1}{r},Z\right), \end{equation}

we can suppress unstable disturbances. It follows that

(4.42)\begin{equation} H_0=\left(\frac{1}{r^2}+Z^2\right)^{1/2}, \quad H_1=\frac{{\rm i}m+{\rm i}kZr^2}{r^2H_0}\,\hat{\phi}(r)\,\zeta, \end{equation}

and (4.7) in component form is

(4.43a)$$\begin{gather} (s\,{{Re}}-\mathcal{L})\hat{\omega}_r+\frac{2{\rm i}m}{r^2}\,\hat{\omega}_{\theta} = 0, \end{gather}$$

(4.43b)$$\begin{gather}(s\,{{Re}}-\mathcal{L})\hat{\omega}_{\theta}-\frac{2{\rm i}m}{r^2}\,\hat{\omega}_{r} = {B}\left(\frac{k(m+kZr^2)\chi_0'\hat{\phi}}{r^2}+{\rm i}kH_0H_0'\hat{\chi}\right), \end{gather}$$

(4.43c)$$\begin{gather}(s\,{{Re}}-\mathcal{D})\hat{\omega}_z = {B}\left(-\frac{m(m+kZr^2)\chi_0'\hat{\phi}}{r^3}-\frac{{\rm i}mH_0H_0'\hat{\chi}}{r}\right). \end{gather}$$

4.4.1. Axisymmetric disturbances

For solely axisymmetric disturbances, we have

(4.44)\begin{equation} H_1=\frac{{\rm i}kZ}{H_0}\,\hat{\phi}(r)\exp({{\rm i}(m\theta+kz)+st}), \end{equation}

and (4.43) becomes

(4.45a)$$\begin{gather} (s\,{{Re}}-\mathcal{L}_0)\hat{\omega}_r = 0, \end{gather}$$

(4.45b)$$\begin{gather}(s\,{{Re}}-\mathcal{L}_0)\hat{\omega}_{\theta} = {B}(k^2Z\chi_0'\hat{\phi}+{\rm i}kH_0H_0'\hat{\chi}), \end{gather}$$

(4.45c)$$\begin{gather}(s\,{{Re}}-\mathcal{D}_0)\hat{\omega}_z = 0. \end{gather}$$

Equation (2.4) gives

(4.46)\begin{equation} \hat{\chi}=\frac{{\rm i}((1+\chi_0){\nabla}^2\hat{\phi}+\chi_0'\hat{\phi}')}{kZ}, \end{equation}

and it follows from (4.15), for axisymmetric disturbances, that

(4.47)\begin{equation} {\hat{\psi}}=\frac{-{\rm i}s\hat{\chi}}{k\chi_0'}, \end{equation}

where $\hat {\psi }$ was defined in § 4.1.

Equations (4.45)–(4.47) give an eigenvalue equation for $\hat {\phi }$:

(4.48)\begin{equation} (s^2\,{{Re}}\,\mathcal{L}_0-s\mathcal{L}_0^2){\mathcal{P}_0\hat{\phi}}=k^2{B}\chi_0'(H_0H_0'\mathcal{P}_0\hat{\phi}-k^2Z^2\hat{\phi}), \end{equation}

where

(4.49)\begin{equation} \mathcal{P}_0\hat{\phi}=\frac{(1+\chi_0)}{\chi_0'}\left(\frac{1}{r}\,(r\hat{\phi}')'-k^2\hat{\phi}\right)+\hat{\phi}'. \end{equation}

When $m=0$,

(4.50)\begin{equation} \hat{\phi}'=\frac{kI_1\hat{\phi}}{(1+\chi_0)I_0} \end{equation}

at $r=a$, and $\boldsymbol {\nabla }\hat {\phi } \rightarrow 0$ as $r\rightarrow \infty$. Due to the boundary conditions for $\hat {\chi }$, $\mathcal {P}_0\hat {\phi }=0$ at $r=a$ and as $r\rightarrow \infty$.

Multiply (4.48) by $r\mathcal {P}_0\hat {\phi }^*$, and integrate over the domain to give

(4.51) \begin{align} \int_a^{\infty}(s^2\,{{Re}}\,\mathcal{P}_0\hat{\phi}^*\mathcal{L}_0{\mathcal{P}_0\hat{\phi}}-s\mathcal{P}_0\hat{\phi}^*\mathcal{L}_0^2{\mathcal{P}_0\hat{\phi}})r\,{\rm d}r=k^2{B}\int_a^{\infty}r\chi_0'(H_0H_0'|\mathcal{P}_0\hat{\phi}|^2-k^2Z^2\hat{\phi}\mathcal{P}_0\hat{\phi}^*){\rm d}r. \end{align}

Due to the self-adjoint property of $\mathcal {L}_0$ and the boundary conditions for $\hat {\phi }$, we can write (4.51) as

(4.52) \begin{align} & \int_a^{\infty}\left(s^2\,{{Re}}\left(|(\mathcal{P}_0\hat{\phi})'|^2 +\left(k^2+\frac{1}{r^2}\right)|{\mathcal{P}_0\hat{\phi}}|^2\right)+s\,|\mathcal{L}_0{\mathcal{P}_0\hat{\phi}}|^2\right)r\,{\rm d}r\nonumber\\ &\quad +k^2{B}\int_a^{\infty}(H_0H_0'\chi_0'|\mathcal{P}_0\hat{\phi}|^2 +k^2Z^2(1+\chi_0)(|\hat{\phi}'|^2+k^2|\hat{\phi}|^2))r\,{\rm d}r=0. \end{align}

Thus for sufficiently large $kZ$, $\mathrm {Re}(s)<0$ provided that $B\neq 0$, therefore axisymmetric modes are stable. Moreover, in the inviscid regime, this result holds for all disturbances $m,k$.