2.1 Equilibrium state and notation

The cylindrical cell under consideration in this paper is sketched in figure 1. We align the axis of symmetry with the $z$ -axis of a cylindrical coordinate system ( $r,\unicode[STIX]{x1D703},z$ ) with unit vectors ( $\boldsymbol{e}_{r},\boldsymbol{e}_{\unicode[STIX]{x1D703}},\boldsymbol{e}_{z}$ ). At equilibrium, the cryolite and aluminium layers (subscripts $i=1,2$ ) occupy the volumes

(2.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\mathcal{V}}_{1}:(r,\unicode[STIX]{x1D703},z)\in [0,R]\times [0,2\unicode[STIX]{x03C0})\times [0,H_{1}],\\ {\mathcal{V}}_{2}:(r,\unicode[STIX]{x1D703},z)\in [0,R]\times [0,2\unicode[STIX]{x03C0})\times [-H_{2},0],\end{array}\right\}\end{eqnarray}$$

where $R$ is the radius of the cell and $H_{i}$ are the heights of the two layers. We denote $\unicode[STIX]{x1D6FF}V_{i}$ the boundary of each fluid domain and split each boundary $\unicode[STIX]{x1D6FF}{\mathcal{V}}_{i}$ as follows; $\unicode[STIX]{x1D6FF}{\mathcal{V}}_{i}=\unicode[STIX]{x1D6F4}_{i}\cup {\mathcal{S}}$ , with $\unicode[STIX]{x1D6F4}_{i}$ the rigid boundary and ${\mathcal{S}}$ the interface at rest. We denote as $\boldsymbol{n}_{i}$ the outward normal on $\unicode[STIX]{x1D6FF}{\mathcal{V}}_{i}.$

The electrical conductivity, density, kinematic viscosity of both fluids are denoted by $\unicode[STIX]{x1D70E}_{i},\unicode[STIX]{x1D70C}_{i},\unicode[STIX]{x1D708}_{i}$ , respectively. The surface tension of the interface between the two fluids is denoted by $\unicode[STIX]{x1D6FE}_{1|2}$ . We suppose that the interface touches the lateral wall with a contact angle of $90^{\circ }$ and that it slides freely.

If the cell is at equilibrium, both fluids are at rest. The fluids are subject to gravity $\boldsymbol{g}=-g\boldsymbol{e}_{z}$ . A homogenous electrical current density $\boldsymbol{J}=J\boldsymbol{e}_{z}$ runs through the cell (in reduction cells, $J<0$ always). This current generates an azimuthal magnetic field $B_{\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D707}_{0}Jr/2$ and we suppose that a uniform external magnetic field $\boldsymbol{B}^{e}=B_{z}\boldsymbol{e}_{z}$ is present. Horizontal components of the external magnetic fields are not considered here because one can show that they have no impact on the instability, see appendix C. The total magnetic field is denoted $\boldsymbol{B}=B_{\unicode[STIX]{x1D703}}\boldsymbol{e}_{\unicode[STIX]{x1D703}}+B_{z}\boldsymbol{e}_{z}$ . The above state is a magnetohydrodynamic equilibrium configuration if the hydrodynamic pressure is $P_{i}=P_{\ast }-\unicode[STIX]{x1D70C}_{i}gz-\unicode[STIX]{x1D707}_{0}J^{2}r^{2}/4$ and the electrical potential is $\unicode[STIX]{x1D6F7}_{i}=\unicode[STIX]{x1D6F7}_{\ast }-Jz/\unicode[STIX]{x1D70E}_{i}$ with $P_{\ast },\unicode[STIX]{x1D6F7}_{\ast }$ arbitrary and constant offsets.

Although it appears as a detail, the hypothesis on the contact angle is in fact rather crucial to our stability theory. A more realistic contact angle, possibly different from $90^{\circ }$ , may result in a very different base state (such as figure 6 of Flueck et al. (Reference Flueck, Hofer, Picasso, Rappaz and Steiner2009)) with a deformed interface, inhomogeneous current density and possibly even capillary induced electro-vortex flows. The stability analysis of this kind of base state is evidently much more difficult, and this explains why we have limited our model to the simplified case of a freely sliding contact angle of $90^{\circ }$ . This configuration admits a simple base state with a flat interface, homogenous current and no flow. This implies that we only deal with some possible effects of capillarity but certainly not all.

2.2 The instability mechanism that we model

Figure 2.

Illustration of a metal pad roll instability mechanism in an idealized cylindrical reduction cell. (a) Equilibrium cell with electrolysis current $\boldsymbol{J}$ and external vertical magnetic field $\boldsymbol{B}^{e}=B_{z}\boldsymbol{e}_{z}$ . (b) Snapshot of a small amplitude rotating wave with rotation speed $\unicode[STIX]{x1D714}$ and accompanying flow $\boldsymbol{u}$ in top and bottom layers. (c) Due to the wave, current density redistributes to preferentially pass through the shallower parts of the cryolite. (d) The current excess $\boldsymbol{j}$ is horizontal in the aluminium and interacts there with $\boldsymbol{B}^{e}$ to create a Lorentz force $\boldsymbol{j}\times \boldsymbol{B}^{e}$ that pushes in the direction of the flow $\boldsymbol{u}$ , amplifying the rotating wave as a result.

In figure 2, we illustrate the instability mechanism that is captured by our model. Overall, it is very similar to Sele’s mechanism (Sele Reference Sele1977), but the cylindrical geometry allows an essential, simplifying difference in presentation. Rather than finding the amplification as a consequence of a coupling between transverse standing waves, we can, in a cylindrical geometry, more simply show that rotating waves amplify themselves.

Let us start with figure 2(a), that shows the equilibrium configuration: cryolite in the top layer, aluminium in the bottom layer and the reduction current flowing from top to bottom ( $J<0$ ). We choose the external vertical magnetic field $\boldsymbol{B}^{e}$ to be upward. In (b), we perturb this equilibrium with a small rotating wave, rotating clockwise when seen from above, as suggested by the figure. Along with this wave comes the instantaneous fluid flow $\boldsymbol{u}$ , pointing as sketched, either towards (top) or away from (bottom) the observer. Through time, this entire pattern will rotate, and, in the absence of electromagnetic forces, slowly decay in magnitude because of viscous dissipation. In (c,d), we explain how the Lorentz force associated with this wave necessarily amplifies this wave. Since cryolite is a much worse conductor than aluminium, the current density $\boldsymbol{J}$ will deviate due to a current excess $\boldsymbol{j}$ , in order to preferentially pass through the shallower part of the cryolite. This is suggested by the thin and thick arrows of (c). (d) Shows how the current excess $\boldsymbol{j}$ typically loops in reduction cells, being mainly horizontal in the aluminium and mainly vertical in the cryolite. In the aluminium layer, this horizontal $\boldsymbol{j}$ interacts with the external vertical field $\boldsymbol{B}^{e}$ to create a Lorentz force $\boldsymbol{j}\times \boldsymbol{B}^{e}$ which pushes exactly in the direction of the instantaneous flow $\boldsymbol{u}$ . Precisely the same picture can be made at any rotated interface position and this means that positive power is permanently being injected into the rotating wave by the Lorentz force: $\boldsymbol{u}\boldsymbol{\cdot }(\boldsymbol{j}\times \boldsymbol{B}^{e})>0$ . Initially small rotating waves can therefore be electromagnetically amplified. Notice that the direction of rotation of the wave is a crucial element in the destabilization process: inverting $\unicode[STIX]{x1D714}$ , inverts the flow $\boldsymbol{u}$ and with $\boldsymbol{u}\boldsymbol{\cdot }(\boldsymbol{j}\times \boldsymbol{B}^{e})<0$ , the wave rotating in the opposite sense will be damped rather than magnified. Only waves rotating in one direction will be amplified and this direction obviously depends on the sign of $JB_{z}$ . In our theory, this amplification is captured in the growth rate denoted  $\unicode[STIX]{x1D706}_{\mathit{v}}$ .

With the same kind of physical reasoning, one can understand why the azimuthal magnetic field, $B_{\unicode[STIX]{x1D703}}\boldsymbol{e}_{\unicode[STIX]{x1D703}}$ , cannot amplify rotating waves (nor standing waves). The Lorentz force $\boldsymbol{j}\times (B_{\unicode[STIX]{x1D703}}\boldsymbol{e}_{\unicode[STIX]{x1D703}})$ is radial in the cryolite and vertical in the aluminium layer, i.e. always perpendicular to the instantaneous flow $\boldsymbol{u}$ . The azimuthal field can therefore not cause an amplification nor a damping. In our theory we find that the azimuthal field causes a shift in the frequency of the wave (denoted by $\unicode[STIX]{x1D6FF}_{\mathit{h}}$ ), in agreement with the theoretical predictions of Sneyd (Reference Sneyd1985) and Sneyd & Wang (Reference Sneyd and Wang1994).

2.3 Linearized problem in the quasi-static limit

Now that we know how the equilibrium state can become unstable, we want to perform the linear stability analysis that takes into account viscous dissipation, magnetic dissipation and surface tension, and without making the shallow layer approximation. We start by linearizing the governing equations about the equilibrium state. Denoting $B=\max (\unicode[STIX]{x1D707}_{0}JR,B_{z})$ the typical scale for the magnetic field strength, we suppose that the Lundquist numbers $Lu_{i}$ in both layers $i=1,2$ remain very low,

(2.2) $$\begin{eqnarray}Lu_{i}=\unicode[STIX]{x1D70E}_{i}\unicode[STIX]{x1D707}_{0}BR/\sqrt{\unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x1D707}_{0}}\ll 1.\end{eqnarray}$$

This hypothesis allows us to use the quasi-static approximation of MHD, in which the electrical field derives from the electrical potential. Denoting by $(\boldsymbol{u}_{i},p_{i},\boldsymbol{b}_{i},\boldsymbol{j}_{i},\unicode[STIX]{x1D711}_{i})$ perturbations for the velocity, the pressure, the magnetic induction, the current density and the electrical potential, we find the linearized perturbation equations of the magnetohydrodynamic stability problem to be

(2.3a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x2202}_{t}\boldsymbol{u}_{i}+\unicode[STIX]{x1D735}p_{i} & = & \displaystyle \boldsymbol{J}\times \boldsymbol{b}_{i}+\boldsymbol{j}_{i}\times \boldsymbol{B}+\unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x1D708}_{i}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{i},\end{eqnarray}$$

(2.3b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{i} & = & \displaystyle 0,\end{eqnarray}$$

(2.3c ) $$\begin{eqnarray}\displaystyle \boldsymbol{j}_{i} & = & \displaystyle \unicode[STIX]{x1D70E}_{i}(-\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}_{i}+\boldsymbol{u}_{i}\times \boldsymbol{B}),\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{j}_{i}=0,\end{eqnarray}$$

(2.3d ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\times \boldsymbol{b}_{i} & = & \displaystyle \unicode[STIX]{x1D707}_{0}\boldsymbol{j}_{i},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{b}_{i}=0.\end{eqnarray}$$

$\unicode[STIX]{x1D6F4}_{i}$

(2.4) $$\begin{eqnarray}\boldsymbol{u}_{i}\boldsymbol{\cdot }\boldsymbol{n}_{i}=0|_{\unicode[STIX]{x1D6F4}_{i}},\end{eqnarray}$$

to which we must add $\boldsymbol{u}_{\boldsymbol{i}}\times \boldsymbol{n}_{i}=0|_{\unicode[STIX]{x1D6F4}_{i}}$ when the viscous effects are accounted for. On the interface ${\mathcal{S}}$ , which we locate at $z=\unicode[STIX]{x1D702}(r,\unicode[STIX]{x1D703},t)$ , the essential inviscid boundary conditions are

(2.5a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D702} & = & \displaystyle u_{1,z}|_{z=0},\end{eqnarray}$$

(2.5b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D702} & = & \displaystyle u_{2,z}|_{z=0},\end{eqnarray}$$

(2.5c ) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x1D6FE}_{1|2}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D702}+(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g\unicode[STIX]{x1D702} & = & \displaystyle p_{2}|_{z=0}-p_{1}|_{z=0}.\end{eqnarray}$$

$\unicode[STIX]{x1D6F4}_{i}$

(2.6a ) $$\begin{eqnarray}\displaystyle j_{1,r}|_{r=R} & = & \displaystyle 0,\end{eqnarray}$$

(2.6b ) $$\begin{eqnarray}\displaystyle j_{2,r}|_{r=R} & = & \displaystyle 0,\end{eqnarray}$$

(2.6c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D711}_{1}|_{z=H_{1}} & = & \displaystyle 0,\end{eqnarray}$$

(2.6d ) $$\begin{eqnarray}\displaystyle j_{2,z}|_{z=-H_{2}} & = & \displaystyle 0.\end{eqnarray}$$

${\mathcal{S}}$

$\boldsymbol{n}\boldsymbol{\cdot }(\boldsymbol{J}+\boldsymbol{j}_{i})$

$\unicode[STIX]{x1D6F7}_{i}+\unicode[STIX]{x1D711}_{i}$

$z=0$

(2.7a ) $$\begin{eqnarray}\displaystyle 0 & = & \displaystyle j_{1,z}|_{z=0}-j_{2,z}|_{z=0},\end{eqnarray}$$

(2.7b ) $$\begin{eqnarray}\displaystyle J(\unicode[STIX]{x1D70E}_{2}^{-1}-\unicode[STIX]{x1D70E}_{1}^{-1})\unicode[STIX]{x1D702} & = & \displaystyle \unicode[STIX]{x1D711}_{2}|_{z=0}-\unicode[STIX]{x1D711}_{1}|_{z=0}.\end{eqnarray}$$

$\unicode[STIX]{x1D711}_{i}$

$\boldsymbol{j}_{i}$

Finally, for the magnetic induction driven by the electrical currents, we suppose the following boundary condition on the cylindrical boundary:

(2.8) $$\begin{eqnarray}\boldsymbol{e}_{r}\times \boldsymbol{b}_{i}|_{r=R}=0.\end{eqnarray}$$

Magnetic field behaviour at other boundaries is not specified, because it has no influence in the theoretical model that we develop. This boundary condition is obviously artificial but it avoids us having to explicitly calculate the magnetic field perturbation, which is never an easy task. The magnetic field perturbation $\boldsymbol{b}_{i}$ is anyway of secondary importance: it only appears through the term $\boldsymbol{J}\times \boldsymbol{b}_{i}$ that is not destabilizing according to our model. Many previous metal pad roll theories have directly ignored this term and according to Sun et al. (Reference Sun, Zikanov and Ziegler2004), this seems justified.

The linearized problem is now completely defined and so we can look for linear eigenmodes. Due to the stationarity of the equilibrium state, we look for fundamental solutions with a simple exponential time dependence. More precisely, if $f$ represents any of the relevant fields or components, we propose the ansatz

(2.9) $$\begin{eqnarray}f=\hat{f}\text{e}^{st},\end{eqnarray}$$

and try to find the admissible (complex) eigenfunctions $\hat{f}$ and (complex) eigenvalues $s\in \mathbb{C}$ . We call $\text{Re}(s)$ the growth rate and $\text{Im}(s)$ the frequency. If $\text{Re}(s)>0$ , the eigenmode is linearly unstable and if $\text{Im}(s)\neq 0$ the solution is wave like. Using a numerical approach it is possible to solve this linear eigenvalue problem. In the following sections we rather use a perturbative approach to find a suitable approximation to $s$ , valid in some parameter range.

2.4 Free inviscid gravito-capillary waves

Our perturbative approach is based on the free inviscid gravito-capillary waves that form a family of eigensolutions in the absence of an electrolysis current ( $J=0$ ), vertical magnetic field ( $B_{z}=0$ ) and without viscosity ( $\unicode[STIX]{x1D708}=0$ ). In that case, the following hydrodynamic field

(2.10) $$\begin{eqnarray}[\boldsymbol{u}_{i},p_{i},\unicode[STIX]{x1D702}]=[\widehat{\boldsymbol{u}}_{i},\widehat{p}_{i},\widehat{\unicode[STIX]{x1D702}}]\text{e}^{\text{i}\unicode[STIX]{x1D714}t},\end{eqnarray}$$

is a solution to the linearized problem defined above. The hydrodynamic fields $\widehat{\boldsymbol{u}}_{i},\widehat{p}_{i},\widehat{\unicode[STIX]{x1D702}}$ correspond to the classical potential wave solution (see Lamb Reference Lamb1945),

(2.11) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}\widehat{\boldsymbol{u}}_{i}\\ \widehat{p}_{i}\\ \widehat{\unicode[STIX]{x1D702}}\end{array}\right]=\left[\begin{array}{@{}c@{}}\unicode[STIX]{x1D735}\widehat{\unicode[STIX]{x1D719}}_{i}\\ -\unicode[STIX]{x1D70C}_{i}(\text{i}\unicode[STIX]{x1D714})\widehat{\unicode[STIX]{x1D719}}_{i}\\ \unicode[STIX]{x2202}_{z}\widehat{\unicode[STIX]{x1D719}}_{i}/(\text{i}\unicode[STIX]{x1D714})\end{array}\right],\end{eqnarray}$$

with

(2.12) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}\widehat{\unicode[STIX]{x1D719}}_{1}\\ \widehat{\unicode[STIX]{x1D719}}_{2}\\ \widehat{\unicode[STIX]{x1D702}}\end{array}\right]=A\,R\left[\begin{array}{@{}r@{}}(\unicode[STIX]{x1D714}/k)\cosh (k(z-H_{1}))/\sinh (kH_{1})\\ -(\unicode[STIX]{x1D714}/k)\cosh (k(z+H_{2}))/\sinh (kH_{2})\\ \text{i}\end{array}\right]\text{J}_{m}(kr)\text{e}^{\text{i}m\unicode[STIX]{x1D703}}.\end{eqnarray}$$

The non-dimensional amplitude $A$ is arbitrary here. The functions $\text{J}_{m}$ are Bessel functions, $m\in \mathbb{Z}$ is the azimuthal wavenumber and $k$ is a radial wavenumber that takes the discrete values

(2.13a,b ) $$\begin{eqnarray}k=\unicode[STIX]{x1D705}_{mn}/R,\quad \text{J}_{m}^{\prime }(\unicode[STIX]{x1D705}_{mn})=0,\end{eqnarray}$$

to ensure that the boundary condition $u_{r}|_{r=R}=0$ is satisfied. We note by $\unicode[STIX]{x1D705}_{mn}$ the $n$ th zero of the derivative of the Bessel function $\text{J}_{m}$ ( $\unicode[STIX]{x1D705}_{11}\approx 1.84$ , $\unicode[STIX]{x1D705}_{21}\approx 3.05$ , $\unicode[STIX]{x1D705}_{31}\approx 4.20$ ). The natural frequency $\unicode[STIX]{x1D714}$ of the (free) waves is

(2.14) $$\begin{eqnarray}\unicode[STIX]{x1D714}=\pm \sqrt{\frac{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})gk+\unicode[STIX]{x1D6FE}_{1|2}k^{3}}{{\displaystyle \frac{\unicode[STIX]{x1D70C}_{2}}{\tanh (kH_{2})}}+{\displaystyle \frac{\unicode[STIX]{x1D70C}_{1}}{\tanh (kH_{1})}}}},\end{eqnarray}$$

and may be either positive or negative for a fixed $m$ . The solution presented here (see (2.10)) has a simple harmonic form, but the real flow is always real valued. Taking the real part of (2.10), one finds rotating waves with positive rotation speed if $m\unicode[STIX]{x1D714}<0$ (anti-clockwise when seen from above) and negative rotation speed if $m\unicode[STIX]{x1D714}>0$ (clockwise when seen from above). Two oppositely rotating waves can also be superposed to give standing waves.

2.5 Sufficient conditions to apply perturbation methods

Perturbation methods may only be used when the Lorentz force and the viscous forces appear as small perturbations in the momentum balance. Symbolically, this means that

(2.15) $$\begin{eqnarray}[\boldsymbol{J}\times \boldsymbol{b}_{i}],[\boldsymbol{j}_{i}\times \boldsymbol{B}],[\unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x1D708}_{i}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{i}]\ll [\unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x2202}_{t}\boldsymbol{u}_{i}],[\unicode[STIX]{x1D735}p_{i}],\end{eqnarray}$$

in terms of orders of magnitude (denoted using square brackets). This limits the parameter space that can be explored and we want to estimate these limits here.

For simplicity, we assume $J,B,\unicode[STIX]{x1D714}\geqslant 0$ to avoid the use of absolute values. We use $[\boldsymbol{x}]=R$ as the space scale rather than $k^{-1}$ since this is more explicit (and anyway similar to $k^{-1}$ for long-wavelength waves). We denote $[\boldsymbol{u}_{i}]=U$ the velocity scale, which is arbitrary in this linear approach. For hydrodynamic waves, we can estimate

(2.16a,b ) $$\begin{eqnarray}[\unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x2202}_{t}\boldsymbol{u}_{i}]\sim \unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x1D714}U,\quad [p_{i}]\sim \unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x1D714}RU.\end{eqnarray}$$

The viscous term has the following order of magnitude:

(2.17) $$\begin{eqnarray}[\unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x1D708}_{i}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{i}]\sim \unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x1D708}_{i}U/R^{2}.\end{eqnarray}$$

We introduce $B$ as the order of magnitude of the imposed magnetic field

(2.18) $$\begin{eqnarray}[\boldsymbol{B}]\sim B=\left\{\begin{array}{@{}rcl@{}}\unicode[STIX]{x1D707}_{0}JR & \longrightarrow & \text{azimuthal},\\ B_{z} & \longrightarrow & \text{vertical}.\end{array}\right.\end{eqnarray}$$

We can use two independent estimates related to either the azimuthal or to the vertical magnetic field. From Ohm’s law (2.3c ), we then estimate the magnitude of the induction term as

(2.19) $$\begin{eqnarray}[\unicode[STIX]{x1D70E}_{i}\boldsymbol{u}_{i}\times \boldsymbol{B}]\sim \unicode[STIX]{x1D70E}_{i}UB.\end{eqnarray}$$

Still in Ohm’s law (2.3c ), the electrical potential term may have the same order of magnitude as $\boldsymbol{u}_{i}\times \boldsymbol{B}$ , but the interface deformations caused by the waves can also lead to a different magnitude which can be estimated from the boundary conditions (2.5a ), (2.5b ) and (2.7b ). This yields the following two different estimates

(2.20) $$\begin{eqnarray}[-\unicode[STIX]{x1D70E}_{i}\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}_{i}]\sim \left\{\begin{array}{@{}rcl@{}}JU/(\unicode[STIX]{x1D714}R) & \longrightarrow & \text{interface deformations},\\ \unicode[STIX]{x1D70E}_{i}UB & \longrightarrow & \text{induction by the flow}.\end{array}\right.\end{eqnarray}$$

Using Ampère’s law (2.3d ), we can calculate the associated magnetic field perturbation magnitudes

(2.21) $$\begin{eqnarray}[\boldsymbol{b}_{i}]\sim \left\{\begin{array}{@{}l@{}}\unicode[STIX]{x1D707}_{0}JU/\unicode[STIX]{x1D714},\\ \unicode[STIX]{x1D707}_{0}\unicode[STIX]{x1D70E}_{i}URB.\end{array}\right.\end{eqnarray}$$

All the orders of magnitude of the different fields have been estimated so we can now reformulate the inequalities (2.15). The trivial but sufficient condition to model viscous effects perturbatively is

(2.22) $$\begin{eqnarray}[\unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x1D708}_{i}\unicode[STIX]{x1D6FB}^{2}u_{i}]\ll [\unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x2202}_{t}\boldsymbol{u}_{i}]\quad \Leftrightarrow \quad Re_{i}=\frac{\unicode[STIX]{x1D714}R^{2}}{\unicode[STIX]{x1D708}_{i}}\gg 1,\end{eqnarray}$$

which we recognize as a constraint on Reynolds numbers based on the properties of the wave. A set of sufficient conditions that allows us to model the effects of the Lorentz force using perturbation methods is then

(2.23a ) $$\begin{eqnarray}\displaystyle [\boldsymbol{J}\times \boldsymbol{b}_{i}]\ll [\unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x2202}_{t}\boldsymbol{u}_{i}] & \quad \Leftrightarrow \quad & \displaystyle \left\{\begin{array}{@{}l@{}}{\displaystyle \frac{\unicode[STIX]{x1D707}_{0}J^{2}}{\unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x1D714}^{2}}}\ll 1,\\ {\displaystyle \frac{\unicode[STIX]{x1D707}_{0}\unicode[STIX]{x1D70E}_{i}JRB}{\unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x1D714}}}\ll 1,\end{array}\right.\end{eqnarray}$$

(2.23b ) $$\begin{eqnarray}\displaystyle [\boldsymbol{j}_{i}\times \boldsymbol{B}]\ll [\unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x2202}_{t}\boldsymbol{u}_{i}] & \quad \Leftrightarrow \quad & \displaystyle \left\{\begin{array}{@{}l@{}}{\displaystyle \frac{JB}{\unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x1D714}^{2}R}}\ll 1,\\ {\displaystyle \frac{\unicode[STIX]{x1D70E}_{i}B^{2}}{\unicode[STIX]{x1D70C}_{i}\unicode[STIX]{x1D714}}}\ll 1.\end{array}\right.\end{eqnarray}$$

2.6 Perturbative calculation of growth rates and frequency shifts

In this technical section, we calculate the theoretical expression for the complex eigenvalue $s$ that characterizes the exponential behaviour ( ${\sim}\text{exp}(st)$ ) of the linear eigenmode. Since the perturbation methods that we are going to use are rather involved, we first give a brief overview of the methodology. The theoretical expression for $s$ will be split as

(2.24) $$\begin{eqnarray}s=\text{i}\unicode[STIX]{x1D714}+\underbrace{(\unicode[STIX]{x1D706}+\text{i}\unicode[STIX]{x1D6FF})}_{\unicode[STIX]{x1D6FC}}.\end{eqnarray}$$

Here $\unicode[STIX]{x1D714}$ is the inviscid gravity wave frequency and $\unicode[STIX]{x1D6FC}$ is a small complex shift induced by the Lorentz and the viscous forces. We define $\unicode[STIX]{x1D706}=\text{Re}(\unicode[STIX]{x1D6FC})$ , the growth rate of the gravity wave, and $\unicode[STIX]{x1D6FF}=\text{Im}(\unicode[STIX]{x1D6FC})$ , the frequency shift. The growth rate is decomposed into several independent terms as follows:

(2.25) $$\begin{eqnarray}\unicode[STIX]{x1D706}=\underbrace{\unicode[STIX]{x1D706}_{\mathit{v}}}_{destab.}+\underbrace{\unicode[STIX]{x1D706}_{\mathit{vv}}+\unicode[STIX]{x1D706}_{\mathit{hh}}}_{quasi\text{-}static}+\underbrace{\unicode[STIX]{x1D706}_{\mathit{visc}}}_{visc.damping}.\end{eqnarray}$$

The first term $\unicode[STIX]{x1D706}_{\mathit{v}}\lessgtr 0$ is the potentially destabilizing term and directly relates to the instability mechanism sketched in figure 2, where the external vertical magnetic field $B_{z}$ (hence the suffix $~_{\mathit{v}}$ ) plays an important role in destabilizing a rotating wave. The terms $\unicode[STIX]{x1D706}_{\mathit{vv}},\unicode[STIX]{x1D706}_{\mathit{hh}}$ are negative in general. These corrections are quadratic in the vertical and horizontal magnetic field magnitudes (hence the suffixes $~_{\mathit{vv}}$ and $~_{\mathit{hh}}$ ) and are often ignored in metal pad roll theory. They relate to the inductive terms $\unicode[STIX]{x1D70E}_{i}\boldsymbol{u}_{i}\times \boldsymbol{B}$ that exist in the quasi-static version of Ohm’s law and that in general cause a magnetic damping. The term $\unicode[STIX]{x1D706}_{\mathit{visc}}<0$ captures the viscous damping of the gravity waves. We similarly decompose the frequency shift as follows:

(2.26) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D6FF}_{\mathit{h}}+\unicode[STIX]{x1D6FF}_{\mathit{vh}}+\unicode[STIX]{x1D6FF}_{\mathit{visc}}.\end{eqnarray}$$

The term $\unicode[STIX]{x1D6FF}_{\mathit{h}}$ is the result of the Lorentz force interaction of $\boldsymbol{j}_{i}$ , caused by the interface deformation, with the background horizontal (azimuthal) magnetic field $B_{\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D707}_{0}Jr/2$ (hence the suffix $~_{\mathit{h}}$ ). The term $\unicode[STIX]{x1D6FF}_{\mathit{vh}}$ is a small quasi-static correction related to both $B_{\unicode[STIX]{x1D703}}$ and $B_{z}$ , and $\unicode[STIX]{x1D6FF}_{\mathit{visc}}$ is a frequency shift caused by viscous effects.

Above, we have provided sufficient conditions to apply a perturbation method. Necessary conditions can be derived a posteriori by verifying that

(2.27) $$\begin{eqnarray}\left|\frac{\unicode[STIX]{x1D706}_{\mathit{v}}}{\unicode[STIX]{x1D714}}\right|,\quad \left|\frac{\unicode[STIX]{x1D706}_{\mathit{vv}}}{\unicode[STIX]{x1D714}}\right|,\quad \left|\frac{\unicode[STIX]{x1D706}_{\mathit{hh}}}{\unicode[STIX]{x1D714}}\right|,\quad \left|\frac{\unicode[STIX]{x1D706}_{\mathit{visc}}}{\unicode[STIX]{x1D714}}\right|,\quad \left|\frac{\unicode[STIX]{x1D6FF}_{\mathit{h}}}{\unicode[STIX]{x1D714}}\right|,\quad \left|\frac{\unicode[STIX]{x1D6FF}_{\mathit{vh}}}{\unicode[STIX]{x1D714}}\right|,\quad \left|\frac{\unicode[STIX]{x1D6FF}_{\mathit{visc}}}{\unicode[STIX]{x1D714}}\right|\ll 1.\end{eqnarray}$$

When these conditions are not met, the theoretical formulas we deduce in the following sections may not be accurate or may even be wrong. This observation must be kept in mind.

2.6.1 Metal pad roll instability in the magneto-static and inviscid limit: $\unicode[STIX]{x1D706}_{\mathit{v}}$ and $\unicode[STIX]{x1D6FF}_{\mathit{h}}$

In this section, we focus on the leading terms in the growth rate and frequency shift formulas, $\unicode[STIX]{x1D706}_{\mathit{v}}$ and $\unicode[STIX]{x1D6FF}_{\mathit{h}}$ . This can be done in the inviscid limit $\unicode[STIX]{x1D708}_{i}=0$ , and, assuming that the magnetostatic approximation of Ohm’s law

(2.28) $$\begin{eqnarray}\boldsymbol{j}_{i}\approx -\unicode[STIX]{x1D70E}_{i}\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}_{i}\end{eqnarray}$$

can be used. This a common supplementary assumption in most metal pad roll theory. Using $[\unicode[STIX]{x1D70E}_{i}\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}_{i}]=JU/(\unicode[STIX]{x1D714}R)$ and $[\unicode[STIX]{x1D70E}_{i}\boldsymbol{u}_{i}\times \boldsymbol{B}]=\unicode[STIX]{x1D70E}_{i}UB$ , from our order of magnitude analysis we find a sufficient condition

(2.29) $$\begin{eqnarray}[\unicode[STIX]{x1D70E}_{i}\boldsymbol{u}_{i}\times \boldsymbol{B}]\ll [-\unicode[STIX]{x1D70E}_{i}\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}_{i}]\quad \Leftrightarrow \quad \unicode[STIX]{x1D70E}_{i}\unicode[STIX]{x1D714}RB/J\ll 1,\end{eqnarray}$$

which we suppose valid in this first section. As a result of the inviscid and magneto-static hypothesis, both viscous and Joule dissipation will be absent in this entire section.

We propose the following ansatz to the inviscid and magneto-static limit of the linearized problem defined in § 2.3:

(2.30) $$\begin{eqnarray}[\boldsymbol{u}_{i},p_{i},\unicode[STIX]{x1D702},\boldsymbol{j}_{i},\unicode[STIX]{x1D711}_{i},\boldsymbol{b}_{i}]=\left(\left[\widehat{\boldsymbol{u}}_{i},\widehat{p}_{i},\widehat{\unicode[STIX]{x1D702}},\widehat{\boldsymbol{j}}_{i},\widehat{\unicode[STIX]{x1D711}}_{i},\widehat{\boldsymbol{b}}_{i}\right]+\left[\widetilde{\boldsymbol{u}}_{i},\widetilde{p}_{i},\widetilde{\unicode[STIX]{x1D702}},\widetilde{\boldsymbol{j}}_{i},\widetilde{\unicode[STIX]{x1D711}}_{i},\widetilde{\boldsymbol{b}}_{i}\right]\right)\text{e}^{\text{i}\unicode[STIX]{x1D714}t}\text{e}^{\unicode[STIX]{x1D6FC}t}.\end{eqnarray}$$

In this notation, the tilded variables are small perturbations with respect to the hatted variables, and the complex shift $\unicode[STIX]{x1D6FC}$ is small with respect to $\unicode[STIX]{x1D714}$ . Here $\unicode[STIX]{x1D714}$ is the frequency of the gravity waves. We inject this ansatz into the governing equations and we treat the Lorentz force term perturbatively. At leading order, we find the following balance equations:

(2.31) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D70C}_{i}(i\unicode[STIX]{x1D714})\widehat{\boldsymbol{u}}_{i}+\unicode[STIX]{x1D735}\widehat{p}_{i}=0,\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\widehat{\boldsymbol{u}}_{i}=0,\\ \widehat{\boldsymbol{j}}_{i}=-\unicode[STIX]{x1D70E}_{i}\unicode[STIX]{x1D735}\widehat{\unicode[STIX]{x1D711}}_{i},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\widehat{\boldsymbol{j}}_{i}=0,\\ \unicode[STIX]{x1D735}\times \widehat{\boldsymbol{b}}_{i}=\unicode[STIX]{x1D707}_{0}\widehat{\boldsymbol{j}}_{i},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\widehat{\boldsymbol{b}}_{i}=0.\end{array}\right\}\end{eqnarray}$$

Notice that the Lorentz force is absent in the momentum equation. The boundary conditions for the hatted fields can be copied from ((2.4) $\rightarrow$ (2.8)) after replacing $\unicode[STIX]{x2202}_{t}\rightarrow \text{i}\unicode[STIX]{x1D714}$ . We have to solve this problem in order to specify all the first-order quantities. The hydrodynamic fields $\widehat{\boldsymbol{u}}_{i},\widehat{p}_{i},\widehat{\unicode[STIX]{x1D702}}$ and the frequency $\unicode[STIX]{x1D714}$ have already been given above. The electrical potential $\widehat{\unicode[STIX]{x1D711}}_{i}$ is harmonic ( $\unicode[STIX]{x1D6FB}^{2}\widehat{\unicode[STIX]{x1D711}}_{i}=0$ ) and is related to the surface elevation $\widehat{\unicode[STIX]{x1D702}}$ by the jump condition (2.7b ). After some computations we find that

(2.32) $$\begin{eqnarray}\displaystyle \left[\begin{array}{@{}c@{}}\widehat{\unicode[STIX]{x1D711}}_{1}\\ \widehat{\unicode[STIX]{x1D711}}_{2}\end{array}\right] & = & \displaystyle JRA\frac{\unicode[STIX]{x1D70E}_{1}^{-1}-\unicode[STIX]{x1D70E}_{2}^{-1}}{\unicode[STIX]{x1D70E}_{1}^{-1}\tanh (kH_{1})+(\unicode[STIX]{x1D70E}_{2}\tanh (kH_{2}))^{-1}}\nonumber\\ \displaystyle & & \displaystyle \times \left[\begin{array}{@{}c@{}}-\text{i}\,\sinh (k(z-H_{1}))/(\unicode[STIX]{x1D70E}_{1}\cosh (kH_{1}))\\ -\text{i}\,\cosh (k(z+H_{2}))/(\unicode[STIX]{x1D70E}_{2}\sinh (kH_{2}))\end{array}\right]\text{J}_{m}(kr)\text{e}^{\text{i}m\unicode[STIX]{x1D703}}.\end{eqnarray}$$

The insulating side wall condition is automatically satisfied since $\text{J}_{m}^{\prime }(kR)=0$ by definition (see (2.13)). The electrical top and bottom boundary conditions are satisfied as well. The electrical currents, $\widehat{\boldsymbol{j}}_{i}=-\unicode[STIX]{x1D70E}_{i}\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}_{i}$ are easily deduced from the potentials. The magnetic field perturbations $\widehat{\boldsymbol{b}}_{i}$ can be calculated by using for example the Biot–Savart law. This computation is very technical, but it is not necessary here: $\widehat{\boldsymbol{b}}_{i}$ is not explicitly required to evaluate the complex shift $\unicode[STIX]{x1D6FC}$ , provided that the boundary condition (2.8) applies.

It is instructive to notice the systematic real or purely imaginary character of the different leading-order fields after dividing by $A\exp (\text{i}m\unicode[STIX]{x1D703})$ ,

(2.33a,b ) $$\begin{eqnarray}\frac{\{\widehat{\unicode[STIX]{x1D719}}_{i},\widehat{u}_{i,r},\widehat{u}_{i,z},\widehat{j}_{i,\unicode[STIX]{x1D703}},\widehat{b}_{i,r},\widehat{b}_{i,z}\}}{A\exp (\text{i}m\unicode[STIX]{x1D703})}\in \mathbb{R},\quad \frac{\{\widehat{\unicode[STIX]{x1D702}}_{i},\widehat{u}_{i,\unicode[STIX]{x1D703}},\widehat{\unicode[STIX]{x1D711}}_{i},\widehat{j}_{i,r},\widehat{j}_{i,z},\widehat{b}_{i,\unicode[STIX]{x1D703}}\}}{A\exp (\text{i}m\unicode[STIX]{x1D703})}\in \text{i}\mathbb{R}.\end{eqnarray}$$

All the leading-order fields that are needed for the calculation of the shift $\unicode[STIX]{x1D6FC}$ are now specified. We are now ready to write the next order in the perturbation problem. From the balance of momentum and mass, we have

(2.34a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{i}\widehat{\boldsymbol{u}}_{i}+\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{i}\widetilde{\boldsymbol{u}}_{i}+\unicode[STIX]{x1D735}\widetilde{p}_{i} & = & \displaystyle \widehat{\boldsymbol{j}}_{i}\times \boldsymbol{B}+\boldsymbol{J}\times \widehat{\boldsymbol{b}}_{i},\end{eqnarray}$$

(2.34b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\widetilde{\boldsymbol{u}}_{i} & = & \displaystyle 0.\end{eqnarray}$$

Note that the (small) Lorentz force is now present on the right-hand side of the momentum balance equation, and it can be calculated by using the first-order fields $\widehat{\boldsymbol{j}}_{i}$ and $\widehat{\boldsymbol{b}}_{i}$ . We are not considering the corrections $\widetilde{\boldsymbol{j}}_{i},\widetilde{\unicode[STIX]{x1D711}}_{i},\widetilde{\boldsymbol{b}}_{i}$ at this stage as they do not affect the eigenvalue shift $\unicode[STIX]{x1D6FC}$ . The hydrodynamic boundary conditions for the tilded variables can be copied from (2.4) and (2.5) except for the kinematic boundary conditions which now reads

(2.35) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}\widehat{\unicode[STIX]{x1D702}}+\text{i}\unicode[STIX]{x1D714}\widetilde{\unicode[STIX]{x1D702}}=\widetilde{u}_{i,z}|_{z=0}.\end{eqnarray}$$

Next to $\text{i}\unicode[STIX]{x1D714}\widetilde{\unicode[STIX]{x1D702}}$ , we find the term $\unicode[STIX]{x1D6FC}\widehat{\unicode[STIX]{x1D702}}$ that is proportional to the small shift $\unicode[STIX]{x1D6FC}$ . To find an equation for $\unicode[STIX]{x1D6FC}$ , we express the solvability condition (also known as the Fredholm alternative). This condition is obtained by summing over the two fluid regions, $i=1,2$ , the linear combination $\int _{V_{i}}[\widehat{\boldsymbol{u}}_{i}^{\ast }\boldsymbol{\cdot }(2.34\mathit{a})+\widehat{p}_{i}^{\ast }(2.34\mathit{b})]\,\text{d}V$ . We obtain

(2.36) $$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{i=1,2}\int _{V_{i}}[\widehat{\boldsymbol{u}}_{i}^{\ast }\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{i}\widehat{\boldsymbol{u}}_{i}+\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{i}\widetilde{\boldsymbol{u}}_{i}+\unicode[STIX]{x1D735}\widetilde{p}_{i})+\widehat{p}_{i}^{\ast }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\widetilde{\boldsymbol{u}}_{i}]\,\text{d}V\nonumber\\ \displaystyle & & \displaystyle \quad =\mathop{\sum }_{i=1,2}\int _{V_{i}}\widehat{\boldsymbol{u}}_{i}^{\ast }\boldsymbol{\cdot }(\widehat{\boldsymbol{j}}_{i}\times \boldsymbol{B}+\boldsymbol{J}\times \widehat{\boldsymbol{b}}_{i})\,\text{d}V.\end{eqnarray}$$

Through quite a number of simplifications, involving integration by parts, use of leading-order balances and the boundary conditions for the hatted and tilded variables, we end up with a formula for the complex frequency shift

(2.37) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}=\underbrace{\frac{{\mathcal{P}}_{\mathit{v}}}{2{\mathcal{K}}}}_{\unicode[STIX]{x1D706}_{\mathit{v}}}+\underbrace{\frac{{\mathcal{P}}_{\mathit{h}}}{2{\mathcal{K}}}}_{\text{i}\unicode[STIX]{x1D6FF}_{\mathit{v}}},\end{eqnarray}$$

where

(2.38a ) $$\begin{eqnarray}{\mathcal{K}}=[(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g+\unicode[STIX]{x1D6FE}_{1|2}k^{2}]\int _{{\mathcal{S}}}|\widehat{\unicode[STIX]{x1D702}}|^{2}\,\text{d}S\end{eqnarray}$$

is proportional to the kinetic energy and

(2.38b ) $$\begin{eqnarray}\displaystyle {\mathcal{P}}_{\mathit{v}} & = & \displaystyle \mathop{\sum }_{i=1,2}\int _{V_{i}}\widehat{\boldsymbol{u}}_{i}^{\ast }\boldsymbol{\cdot }(\widehat{\boldsymbol{j}}_{i}\times B_{z}\boldsymbol{e}_{z})\,\text{d}V,\end{eqnarray}$$

(2.38c ) $$\begin{eqnarray}\displaystyle {\mathcal{P}}_{\mathit{h}} & = & \displaystyle \mathop{\sum }_{i=1,2}\int _{V_{i}}\left\{\widehat{\boldsymbol{u}}_{i}^{\ast }\boldsymbol{\cdot }\left[\widehat{\boldsymbol{j}}_{i}\times (\unicode[STIX]{x1D707}_{0}Jr/2)\boldsymbol{e}_{\unicode[STIX]{x1D703}}\right]+\unicode[STIX]{x1D707}_{0}J\widehat{\unicode[STIX]{x1D719}}_{i}^{\ast }\,\widehat{j}_{i,z}\right\}\,\text{d}V,\end{eqnarray}$$

to the power injected by the Lorentz force. The condition expressed by this solvability condition is in fact very similar to what the mechanical energy equation would yield at leading order. The mathematical operations behind these simplifications are detailed in appendix A. Partial integration and the artificial magnetic boundary condition (2.8) were necessary to rewrite the power density $\widehat{\boldsymbol{u}}_{i}^{\ast }\boldsymbol{\cdot }(\boldsymbol{J}\times \widehat{\boldsymbol{b}}_{i})$ as the second term $\unicode[STIX]{x1D707}_{0}J\widehat{\unicode[STIX]{x1D719}}_{i}^{\ast }\,\widehat{j}_{i,z}$ in the integrandum of  ${\mathcal{P}}_{\mathit{h}}$ .

Let us give a first interpretation of these formulas. Using the real or purely imaginary nature of the different field components (2.33), we immediately deduce that ${\mathcal{K}},{\mathcal{P}}_{\mathit{v}}$ are real and ${\mathcal{P}}_{\mathit{h}}$ imaginary, explaining why we can split $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D706}_{\mathit{v}}+\text{i}\unicode[STIX]{x1D6FF}_{\mathit{v}}$ . Physically this implies that at this order of the expansion, only the externally imposed, vertical magnetic field component $B_{z}$ can destabilize the gravity waves and the corresponding growth rate is $\unicode[STIX]{x1D706}_{\mathit{v}}$ . The self-generated, azimuthal magnetic field $B_{\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D707}_{0}Jr/2$ is not destabilizing but shifts the frequency of the wave by $\unicode[STIX]{x1D6FF}_{\mathit{h}}$ . This result is coherent with the analysis of Sneyd (Reference Sneyd1985) and Sneyd & Wang (Reference Sneyd and Wang1994) who also found that the self-generated field is not destabilizing. We further note that the artificial magnetic boundary condition (2.8) has no impact on $\unicode[STIX]{x1D706}_{\mathit{v}}$ but it may affect  $\unicode[STIX]{x1D6FF}_{\mathit{h}}$ .

A few pages of calculations are necessary to evaluate all the integrals in ${\mathcal{K}}$ , ${\mathcal{P}}_{\mathit{v}}$ and ${\mathcal{P}}_{\mathit{h}}$ analytically, and this allows us to arrive at the following explicit formulas for the inviscid and magneto-static growth rate and frequency shift:

(2.39a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}_{\mathit{v}} & = & \displaystyle \frac{\unicode[STIX]{x1D714}}{2}\frac{JB_{z}}{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g+\unicode[STIX]{x1D6FE}_{1|2}k^{2}}\frac{m}{(kR)^{2}-m^{2}}\left(\tanh (kH_{1})+\frac{kH_{2}}{\sinh ^{2}(kH_{2})}+\frac{1}{\tanh (kH_{2})}\right)\unicode[STIX]{x1D6EC},\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(2.39b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}_{\mathit{h}} & = & \displaystyle \frac{\unicode[STIX]{x1D714}}{4}\frac{\unicode[STIX]{x1D707}_{0}J^{2}k^{-1}}{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g+\unicode[STIX]{x1D6FE}_{1|2}k^{2}}\frac{(kR)^{2}-2m^{2}}{(kR)^{2}-m^{2}}\left(\frac{kH_{1}}{\sinh (kH_{1})\cosh (kH_{1})}\right)\unicode[STIX]{x1D6EC},\end{eqnarray}$$

where

(2.39c ) $$\begin{eqnarray}\unicode[STIX]{x1D6EC}=\frac{\unicode[STIX]{x1D70E}_{1}^{-1}-\unicode[STIX]{x1D70E}_{2}^{-1}}{\unicode[STIX]{x1D70E}_{1}^{-1}\tanh (kH_{1})+(\unicode[STIX]{x1D70E}_{2}\tanh (kH_{2}))^{-1}}.\end{eqnarray}$$

In both the growth rate $\unicode[STIX]{x1D706}_{\mathit{v}}$ and frequency shift $\unicode[STIX]{x1D6FF}_{\mathit{h}}$ , we see dimensionless factors that balance the typical strength of the Lorentz force ( $JB_{z}$ in $\unicode[STIX]{x1D706}_{\mathit{v}}$ , $\unicode[STIX]{x1D707}_{0}J^{2}k^{-1}$ in $\unicode[STIX]{x1D6FF}_{\mathit{h}}$ ) with respect to gravitational and capillary restoring forces ( $(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g+\unicode[STIX]{x1D6FE}_{1|2}k^{2}$ ). We see that the growth rate $\unicode[STIX]{x1D706}_{\mathit{v}}\sim m$ , so axisymmetric waves with $m=0$ can never be destabilized. The number $\unicode[STIX]{x1D6EC}$ is a conductivity-dependent factor that is not often seen in the metal pad literature. It clearly shows that a jump in electrical conductivity $\unicode[STIX]{x1D70E}_{1}\neq \unicode[STIX]{x1D70E}_{2}$ is essential for instability. Note that we have verified the exactness of the formulas (2.39) by comparing them with numerical evaluations of the integrals of ${\mathcal{P}}_{\mathit{v}}$ , ${\mathcal{P}}_{\mathit{h}}$ , ${\mathcal{K}}$ computed using a simple quadrature rule on a two-dimensional $r-z$ grid.

We can derive an inviscid instability criterion for the case $\unicode[STIX]{x1D70E}_{1}\ll \unicode[STIX]{x1D70E}_{2}$ commonly considered in the literature and realistic for Hall–Héroult cells. Since $(kR)^{2}-m^{2}>0$ always (see appendix A), inviscid instability, that is $\unicode[STIX]{x1D706}_{\mathit{v}}>0$ , requires

(2.40) $$\begin{eqnarray}\text{Sgn}(JB_{z}m\unicode[STIX]{x1D714})>0\quad \Rightarrow \quad \text{wave is MPR-unstable}.\end{eqnarray}$$

In reduction cells, we further have $J<0$ , as the reduction current flows from top to bottom, and if $B_{z}>0$ , this criterion tells us that only rotating waves with $m\unicode[STIX]{x1D714}<0$ can be destabilized. This corresponds to waves rotating in the positive direction, or anti-clockwise when seen from above, and this coincides with what can be expected from the Sele mechanism (Sele Reference Sele1977) or from our presentation of the instability mechanism. Note also that there is no threshold in this dissipationless limit. Waves that are unstable always exist, which is coherent with Davidson & Lindsay (Reference Davidson and Lindsay1998) and Lukyanov et al. (Reference Lukyanov, El and Molokov2001) who find $\unicode[STIX]{x1D6FD}_{c}=0$ in cylindrical cells and without dissipation.

The general formula for the growth rate takes simpler asymptotic forms in the limit of large conductivity jumps and for both deep fluid layers and shallow fluid layers. In the deep limit $kH_{i}\gg 1$ and with $\unicode[STIX]{x1D70E}_{1}/\unicode[STIX]{x1D70E}_{2}\ll 1$ we find

(2.41a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}_{\mathit{v},\mathit{deep}} & {\approx} & \displaystyle \frac{JB_{z}}{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g+\unicode[STIX]{x1D6FE}_{1|2}k^{2}}\frac{m\unicode[STIX]{x1D714}_{\mathit{deep}}}{(kR)^{2}-m^{2}},\end{eqnarray}$$

(2.41b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}_{h,\mathit{deep}} & {\approx} & \displaystyle 0,\end{eqnarray}$$

where

(2.42) $$\begin{eqnarray}\unicode[STIX]{x1D714}_{\mathit{deep}}\approx \sqrt{\frac{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})gk+\unicode[STIX]{x1D6FE}_{1|2}k^{3}}{\unicode[STIX]{x1D70C}_{1}+\unicode[STIX]{x1D70C}_{2}}}.\end{eqnarray}$$

For shallow layers $kH_{i}\ll 1$ and with $\unicode[STIX]{x1D70E}_{1}/\unicode[STIX]{x1D70E}_{2}\ll (kH_{1})(kH_{2})$ we have

(2.43a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}_{\mathit{v},\mathit{shallow}} & {\approx} & \displaystyle \frac{JB_{z}}{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g+\unicode[STIX]{x1D6FE}_{1|2}k^{2}}\frac{m\unicode[STIX]{x1D714}_{\mathit{shallow}}}{(kR)^{2}-m^{2}}\frac{1}{(kH_{1})(kH_{2})},\end{eqnarray}$$

(2.43b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}_{h,\mathit{shallow}} & {\approx} & \displaystyle \frac{\unicode[STIX]{x1D714}_{\mathit{shallow}}}{4}\frac{\unicode[STIX]{x1D707}_{0}J^{2}k^{-1}}{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g+\unicode[STIX]{x1D6FE}_{1|2}k^{2}}\frac{(kR)^{2}-2m^{2}}{(kR)^{2}-m^{2}}\frac{1}{kH_{1}},\end{eqnarray}$$

where

(2.44) $$\begin{eqnarray}\unicode[STIX]{x1D714}_{\mathit{shallow}}\approx \sqrt{\frac{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})gk^{2}+\unicode[STIX]{x1D6FE}_{1|2}k^{4}}{\unicode[STIX]{x1D70C}_{1}H_{1}^{-1}+\unicode[STIX]{x1D70C}_{2}H_{2}^{-1}}}.\end{eqnarray}$$

The shallow approximation is commonly adopted in the metal pad roll literature. For long-wavelength waves, the surface tension is also usually ignored. In this limit, the growth rate can be approximated as follows:

(2.45) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{\mathit{v},\mathit{shallow}}\approx \frac{JB_{z}}{\sqrt{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g(\unicode[STIX]{x1D70C}_{1}H_{1}^{-1}+\unicode[STIX]{x1D70C}_{2}H_{2}^{-1})}}\frac{R}{H_{1}H_{2}}\frac{m}{\unicode[STIX]{x1D705}_{mn}(\unicode[STIX]{x1D705}_{mn}-m^{2})}.\end{eqnarray}$$

In this formula, we have separated a wave-dependent factor that only depends on $\unicode[STIX]{x1D705}_{mn}$ and $m$ , from a factor that groups all other physical and geometrical parameters. This allows us to emphasize that the growth rate increases proportionally to the lateral extent $R$ of the shallow cell if current density $J$ is held fixed. Alternatively, we can also write

(2.46a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{\mathit{v},\mathit{shallow}}\approx \frac{m\unicode[STIX]{x1D714}_{\mathit{shallow}}}{\unicode[STIX]{x1D705}_{mn}^{2}-m^{2}}\frac{\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D705}_{mn}^{2}},\quad \unicode[STIX]{x1D6FD}=\frac{IB_{z}}{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})gH_{1}H_{2}},\end{eqnarray}$$

with $I=J\unicode[STIX]{x03C0}R^{2}$ the total current running through the cell. The shallow-limit growth rate is proportional to the non-dimensional group $\unicode[STIX]{x1D6FD}$ that appeared in many previous studies on the metal pad roll instability and that we discussed in the Introduction. In appendix B, we develop further the small gap limit and show that the growth rate formula (2.46) is compatible with the existing dissipationless shallow layer models. In the double limit of asymptotically small $kH_{i}\ll 1$ and $\unicode[STIX]{x1D6FD}\ll 1$ and without dissipation, our formula and the models of Davidson & Lindsay (Reference Davidson and Lindsay1998) and Lukyanov et al. (Reference Lukyanov, El and Molokov2001) yield precisely the same growth rates.

In the supplementary material to this article (available online athttps://doi.org/10.1017/jfm.2019.642), we explore the effect on the metal pad roll instability of having symmetric boundary conditions at the top and the bottom of the cell (fixed potential or fixed normal current), as investigated in Munger & Vincent (Reference Munger and Vincent2006b ). In the shallow limit, the mixed electrical boundary conditions we use in the present paper (fixed potential at the top and fixed normal current at the bottom) always yield growth rates that are larger than when the boundary conditions are symmetric. In deep cells, the growth rates naturally depend very little on the nature of the top and bottom boundary conditions.

2.6.2 Quasi-static corrections: $\unicode[STIX]{x1D706}_{\mathit{vv}}$ , $\unicode[STIX]{x1D706}_{\mathit{hh}}$ and $\unicode[STIX]{x1D6FF}_{\mathit{vh}}$

When the magnetic field strength $B=\unicode[STIX]{x1D707}_{0}JR~\text{or}~B_{z}$ is large, the simplifying assumption $\unicode[STIX]{x1D70E}_{i}\unicode[STIX]{x1D714}RB/J\ll 1$ is no longer satisfied. In this case, the quasi-static approximation of MHD must include the retroaction of the fluid flow on the electrical current density through induction. The perturbations on the electrical current density must take the following form:

(2.47) $$\begin{eqnarray}\boldsymbol{j}_{i}=\unicode[STIX]{x1D70E}_{i}(-\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}_{i}+\boldsymbol{u}_{i}\times \boldsymbol{B}),\end{eqnarray}$$

instead of the static approximation $\boldsymbol{j}_{i}\approx -\unicode[STIX]{x1D70E}_{i}\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}_{i}$ used in the previous section, at least in the aluminium layer where $\unicode[STIX]{x1D70E}_{2}\gg \unicode[STIX]{x1D70E}_{1}$ . The growth rate, $\unicode[STIX]{x1D706}_{\mathit{v}}$ , and frequency shift, $\unicode[STIX]{x1D6FF}_{\mathit{h}}$ , are then augmented with the quasi-static corrections $\unicode[STIX]{x1D706}_{\mathit{vv}}\sim B_{z}^{2}$ , $\unicode[STIX]{x1D706}_{\mathit{hh}}\sim B_{\unicode[STIX]{x1D703}}^{2}$ and $\unicode[STIX]{x1D6FF}_{\mathit{vh}}\sim B_{\unicode[STIX]{x1D703}}B_{z}$ quadratic in magnetic field magnitude.

To find the quasi-static corrections, we modify the perturbation ansatz (2.30) accordingly,

(2.48) $$\begin{eqnarray}\displaystyle [\boldsymbol{u}_{i},p_{i},\unicode[STIX]{x1D702},\boldsymbol{j}_{i},\unicode[STIX]{x1D711}_{i}] & = & \displaystyle \left(\left[\widehat{\boldsymbol{u}}_{i},\widehat{p}_{i},\widehat{\unicode[STIX]{x1D702}},\widehat{\boldsymbol{j}}_{i}+\widehat{\boldsymbol{{\mathcal{J}}}}_{i}^{\mathit{v}}+\widehat{\boldsymbol{{\mathcal{J}}}}_{i}^{\mathit{h}},\widehat{\unicode[STIX]{x1D711}}_{i}+\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{v}}+\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{h}},\widehat{\boldsymbol{b}}_{i}+\widehat{\boldsymbol{{\mathcal{B}}}}_{i}^{\mathit{v}}+\widehat{\boldsymbol{{\mathcal{B}}}}_{i}^{\mathit{h}}\right]\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\left[\widetilde{\boldsymbol{u}}_{i},\widetilde{p}_{i},\widetilde{\unicode[STIX]{x1D702}},\widetilde{\boldsymbol{j}}_{i},\widetilde{\unicode[STIX]{x1D711}}_{i},\widetilde{\boldsymbol{b}}_{i}\right]\right)\text{e}^{\text{i}\unicode[STIX]{x1D714}t}\text{e}^{\unicode[STIX]{x1D6FC}t}.\end{eqnarray}$$

The fields $\widehat{\boldsymbol{u}}_{i},\widehat{p}_{i},\widehat{\unicode[STIX]{x1D702}},\widehat{\unicode[STIX]{x1D711}}_{i},\widehat{\boldsymbol{j}}_{i},\widehat{\boldsymbol{b}}_{i}$ are identical to those defined in § 2.6.1 and capture all the magnetohydrodynamic effects related to surface elevation. The quasi-static corrections for the electrical potential and the current densities are related to either the vertical (superscripts $^{\mathit{v}}$ ) or the horizontal magnetic field (superscript $^{\mathit{h}}$ ). That is to say, we define the current density corrections as follows:

(2.49a ) $$\begin{eqnarray}\displaystyle \widehat{\boldsymbol{{\mathcal{J}}}}_{i}^{\mathit{v}} & = & \displaystyle \unicode[STIX]{x1D70E}_{i}(-\unicode[STIX]{x1D735}\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{v}}+\widehat{\boldsymbol{u}}_{i}\times B_{z}\boldsymbol{e}_{z}),\end{eqnarray}$$

(2.49b ) $$\begin{eqnarray}\displaystyle \widehat{\boldsymbol{{\mathcal{J}}}}_{i}^{\mathit{h}} & = & \displaystyle \unicode[STIX]{x1D70E}_{i}(-\unicode[STIX]{x1D735}\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{h}}+\widehat{\boldsymbol{u}}_{i}\times (\unicode[STIX]{x1D707}_{0}Jr/2)\boldsymbol{e}_{\unicode[STIX]{x1D703}}),\end{eqnarray}$$

with the constraints: $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\widehat{\boldsymbol{{\mathcal{J}}}}_{i}^{v}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\widehat{\boldsymbol{{\mathcal{J}}}}_{i}^{h}=0$ . This in turn implies that electrical potentials satisfy the following Laplace and Poisson problems:

(2.50a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{v}}=0,\quad \unicode[STIX]{x1D6FB}^{2}\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{h}}=-\unicode[STIX]{x1D707}_{0}J\widehat{u}_{z,i},\end{eqnarray}$$

with the following transmission and inhomogeneous boundary conditions:

(2.51) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\widehat{\unicode[STIX]{x1D6F9}}_{1}^{\mathit{v}}|_{z=H_{1}}=0,\quad \unicode[STIX]{x2202}_{z}\widehat{\unicode[STIX]{x1D6F9}}_{2}^{\mathit{v}}|_{z=-H_{2}}=0,\quad \unicode[STIX]{x2202}_{r}\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{v}}|_{r=R}=B_{z}\widehat{u}_{i,\unicode[STIX]{x1D703}}|_{r=R},\\ \widehat{\unicode[STIX]{x1D6F9}}_{1}^{\mathit{v}}|_{z=0}=\widehat{\unicode[STIX]{x1D6F9}}_{2}^{\mathit{v}}|_{z=0},\quad \unicode[STIX]{x1D70E}_{1}\unicode[STIX]{x2202}_{z}\widehat{\unicode[STIX]{x1D6F9}}_{1}^{\mathit{v}}|_{z=0}=\unicode[STIX]{x1D70E}_{2}\unicode[STIX]{x2202}_{z}\widehat{\unicode[STIX]{x1D6F9}}_{2}^{\mathit{v}}|_{z=0},\end{array}\right\}\end{eqnarray}$$

and

(2.52) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\widehat{\unicode[STIX]{x1D6F9}}_{1}^{\mathit{h}}|_{z=H_{1}}=0,\quad \unicode[STIX]{x2202}_{z}\widehat{\unicode[STIX]{x1D6F9}}_{2}^{\mathit{h}}|_{z=-H_{2}}=(\unicode[STIX]{x1D707}_{0}Jr/2)\widehat{u}_{2,r}|_{z=-H_{2}},\\ \unicode[STIX]{x2202}_{r}\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{h}}|_{r=R}=-(\unicode[STIX]{x1D707}_{0}JR/2)\widehat{u}_{i,z}|_{r=R},\quad \widehat{\unicode[STIX]{x1D6F9}}_{1}^{\mathit{h}}|_{z=0}=\widehat{\unicode[STIX]{x1D6F9}}_{2}^{\mathit{h}}|_{z=0},\\ \unicode[STIX]{x1D70E}_{1}(-\unicode[STIX]{x2202}_{z}\widehat{\unicode[STIX]{x1D6F9}}_{1}^{\mathit{h}}+(\unicode[STIX]{x1D707}_{0}Jr/2)\widehat{u}_{1,r})_{z=0}=\unicode[STIX]{x1D70E}_{2}(-\unicode[STIX]{x2202}_{z}\widehat{\unicode[STIX]{x1D6F9}}_{2}^{\mathit{h}}+(\unicode[STIX]{x1D707}_{0}Jr/2)\widehat{u}_{2,r})_{z=0}.\end{array}\right\}\end{eqnarray}$$

The problems defined here have no simple analytical solution, and below we discuss how we handle this in practice. Without any numerical calculation, we can deduce

(2.53) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\displaystyle \frac{\{\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{h}},\widehat{{\mathcal{J}}}_{i,r}^{\mathit{h}},\widehat{{\mathcal{J}}}_{i,z}^{\mathit{h}},\widehat{{\mathcal{J}}}_{i,\unicode[STIX]{x1D703}}^{\mathit{v}},\widehat{{\mathcal{B}}}_{i,r}^{\mathit{v}},\widehat{{\mathcal{B}}}_{i,z}^{\mathit{v}},\widehat{{\mathcal{B}}}_{i,\unicode[STIX]{x1D703}}^{\mathit{h}}\}}{A\exp (\text{i}m\unicode[STIX]{x1D703})}}\in \mathbb{R},\\ {\displaystyle \frac{\{\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{v}},\widehat{{\mathcal{J}}}_{i,r}^{\mathit{v}},\widehat{{\mathcal{J}}}_{i,z}^{\mathit{v}},\widehat{{\mathcal{J}}}_{i,\unicode[STIX]{x1D703}}^{\mathit{h}},\widehat{{\mathcal{B}}}_{i,r}^{\mathit{h}},\widehat{{\mathcal{B}}}_{i,z}^{\mathit{h}},\widehat{{\mathcal{B}}}_{i,\unicode[STIX]{x1D703}}^{\mathit{v}}\}}{A\exp (\text{i}m\unicode[STIX]{x1D703})}}\in \text{i}\mathbb{R}.\end{array}\right\}\end{eqnarray}$$

Let us briefly suppose that the quasi-static field corrections are known and pursue the calculations to find the impact on $\unicode[STIX]{x1D6FC}$ , the complex frequency shift. In this augmented, quasi-static version of the problem, we find a modified solvability condition,

(2.54) $$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{i=1,2}\int _{V_{i}}[\widehat{\boldsymbol{u}}_{i}^{\ast }\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{i}\widehat{\boldsymbol{u}}_{i}+\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{i}\widetilde{\boldsymbol{u}}_{i}+\unicode[STIX]{x1D735}\widetilde{p}_{i})+\widehat{p}_{i}^{\ast }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\widetilde{\boldsymbol{u}}_{i}]\,\text{d}V\nonumber\\ \displaystyle & & \displaystyle \quad =\mathop{\sum }_{i=1,2}\int _{V_{i}}\widehat{\boldsymbol{u}}_{i}^{\ast }\boldsymbol{\cdot }\left[(\widehat{\boldsymbol{j}}_{i}+\widehat{\boldsymbol{{\mathcal{J}}}}_{i}^{\mathit{v}}+\widehat{\boldsymbol{{\mathcal{J}}}}_{i}^{\mathit{h}})\times \boldsymbol{B}+\boldsymbol{J}\times (\widehat{\boldsymbol{b}}_{i}+\widehat{\boldsymbol{{\mathcal{B}}}}_{i}^{\mathit{v}}+\widehat{\boldsymbol{{\mathcal{B}}}}_{i}^{\mathit{h}})\right]\,\text{d}V.\end{eqnarray}$$

Similar simplifications as before lead to a new formula for the complex frequency shift

(2.55) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}=\underbrace{\frac{{\mathcal{P}}_{\mathit{v}}}{2{\mathcal{K}}}}_{\unicode[STIX]{x1D706}_{\mathit{v}}}+\underbrace{\frac{{\mathcal{Q}}_{\mathit{vv}}}{2{\mathcal{K}}}}_{\unicode[STIX]{x1D706}_{\mathit{vv}}}+\underbrace{\frac{{\mathcal{Q}}_{\mathit{hh}}}{2{\mathcal{K}}}}_{\unicode[STIX]{x1D706}_{\mathit{hh}}}+\underbrace{\frac{{\mathcal{P}}_{\mathit{h}}}{2{\mathcal{K}}}}_{\text{i}\unicode[STIX]{x1D6FF}_{\mathit{v}}}+\underbrace{\frac{{\mathcal{Q}}_{\mathit{vh}}}{2{\mathcal{K}}}}_{\text{i}\unicode[STIX]{x1D6FF}_{\mathit{vh}}},\end{eqnarray}$$

with ${\mathcal{K}},{\mathcal{P}}_{v},{\mathcal{P}}_{h}$ as before and

(2.56a ) $$\begin{eqnarray}\displaystyle {\mathcal{Q}}_{\mathit{vv}} & = & \displaystyle \mathop{\sum }_{i=1,2}\int _{V_{i}}\widehat{\boldsymbol{u}}_{i}^{\ast }\boldsymbol{\cdot }\left(\widehat{\boldsymbol{{\mathcal{J}}}}_{i}^{\mathit{v}}\times B_{z}\boldsymbol{e}_{z}\right)\text{d}V,\end{eqnarray}$$

(2.56b ) $$\begin{eqnarray}\displaystyle {\mathcal{Q}}_{\mathit{hh}} & = & \displaystyle \mathop{\sum }_{i=1,2}\int _{V_{i}}\left\{\widehat{\boldsymbol{u}}_{i}^{\ast }\boldsymbol{\cdot }\left[\widehat{\boldsymbol{{\mathcal{J}}}}_{i}^{\mathit{h}}\times (\unicode[STIX]{x1D707}_{0}Jr/2)\boldsymbol{e}_{\unicode[STIX]{x1D703}}\right]+\unicode[STIX]{x1D707}_{0}J\widehat{\unicode[STIX]{x1D719}}_{i}^{\ast }\,\widehat{{\mathcal{J}}}_{i,z}^{\mathit{h}}\right\}\,\text{d}V,\end{eqnarray}$$

(2.56c ) $$\begin{eqnarray}\displaystyle {\mathcal{Q}}_{\mathit{vh}} & = & \displaystyle \mathop{\sum }_{i=1,2}\int _{V_{i}}\left\{\widehat{\boldsymbol{u}}_{i}^{\ast }\boldsymbol{\cdot }\left[\widehat{\boldsymbol{{\mathcal{J}}}}_{i}^{\mathit{h}}\times B_{z}\boldsymbol{e}_{z}+\widehat{\boldsymbol{{\mathcal{J}}}}_{i}^{\mathit{v}}\times (\unicode[STIX]{x1D707}_{0}Jr/2)\boldsymbol{e}_{\unicode[STIX]{x1D703}}\right]+\unicode[STIX]{x1D707}_{0}J\widehat{\unicode[STIX]{x1D719}}_{i}^{\ast }\,\widehat{{\mathcal{J}}}_{i,z}^{\mathit{v}}\right\}\,\text{d}V.\quad\end{eqnarray}$$

The artificial magnetic boundary condition (2.8) was used to write the formulas for ${\mathcal{Q}}_{\mathit{hh}}$ and ${\mathcal{Q}}_{\mathit{vh}}$ , but does not affect ${\mathcal{Q}}_{\mathit{vv}}$ . Using the real and imaginary character of the fields (2.53), we readily deduce that ${\mathcal{Q}}_{\mathit{vv}},{\mathcal{Q}}_{\mathit{hh}}$ are real whereas ${\mathcal{Q}}_{\mathit{vh}}$ is imaginary, thus yielding the quasi-static corrections $\unicode[STIX]{x1D706}_{\mathit{vv}},\unicode[STIX]{x1D706}_{\mathit{hh}}$ to the growth rate $\unicode[STIX]{x1D706}_{\mathit{v}}$ and the quasi-static correction $\unicode[STIX]{x1D6FF}_{\mathit{vh}}$ to the frequency shift $\unicode[STIX]{x1D6FF}_{h}$ .

In practice, we still need to solve for the fields $\widehat{\unicode[STIX]{x1D6F9}}_{i}^{v}$ and $\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{h}}$ to be able to evaluate these integrals and here we have chosen two different approaches. The first approach is adapted to arbitrary values of $\unicode[STIX]{x1D70E}_{i}$ : we just calculate $\widehat{\unicode[STIX]{x1D6F9}}_{i}^{v}$ and $\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{h}}$ numerically solving the mentioned Laplace and Poisson problems. We use a Fourier expansion in the azimuthal direction; that is, we set $\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{v}}=f_{i}^{\mathit{v}}(r,z)\exp (\text{i}m\unicode[STIX]{x1D703})$ , $\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{h}}=f_{i}^{\mathit{h}}(r,z)\exp (\text{i}m\unicode[STIX]{x1D703})$ , and we solve the two-dimensional problems for $f_{i}^{\mathit{v}}$ and $f_{i}^{\mathit{h}}$ using a standard second-order finite-difference method on a uniform grid with $2\times (M+1)^{2}$ points,

(2.57) $$\begin{eqnarray}\left\{\begin{array}{@{}l@{}}\text{region 1:}\quad (i\unicode[STIX]{x1D6FF}r,j\unicode[STIX]{x1D6FF}z_{1}),\quad \\ \text{region 2:}\quad (i\unicode[STIX]{x1D6FF}r,-j\unicode[STIX]{x1D6FF}z_{2}),\quad \end{array}\right.\quad \forall i,j\in \{0,1,\ldots ,M\}.\end{eqnarray}$$

Here $\unicode[STIX]{x1D6FF}r=R/M$ , $\unicode[STIX]{x1D6FF}z_{1}=H_{1}/M$ , $\unicode[STIX]{x1D6FF}z_{2}=H_{2}/M$ . Radial and vertical derivatives appearing in the boundary conditions are discretized using second-order uncentred finite-difference formulas. At the axis, we enforce the regularity condition $\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{v}}|_{r=0}=\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{h}}|_{r=0}=0$ for $m\neq 0$ . ${\mathcal{Q}}_{\mathit{vv}}$ , ${\mathcal{Q}}_{\mathit{hh}}$ , ${\mathcal{Q}}_{\mathit{vh}}$ are then calculated numerically using the field profiles provided by the finite-difference code, using a simple trapezoidal quadrature.

The need to use a numerical approach to find the quasi-static corrections is slightly unsatisfying. Therefore we also propose a second approach applicable to the realistic limit where $\unicode[STIX]{x1D70E}_{1}\ll \unicode[STIX]{x1D70E}_{2}$ . In that case, the induction in the top cryolite layer is very weak, and we can use the approximation

(2.58) $$\begin{eqnarray}\widehat{\boldsymbol{{\mathcal{J}}}}_{1}^{\mathit{v}}\approx \widehat{\boldsymbol{{\mathcal{J}}}}_{1}^{\mathit{h}}\approx \mathbf{0}.\end{eqnarray}$$

We then need to solve (2.50) in the bottom aluminium layer only, with the boundary conditions stated in (2.51), (2.52), and the transmission conditions replaced by

(2.59a,b ) $$\begin{eqnarray}\unicode[STIX]{x2202}_{z}\widehat{\unicode[STIX]{x1D6F9}}_{2}^{\mathit{v}}|_{z=0}\approx 0,\quad \unicode[STIX]{x2202}_{z}\widehat{\unicode[STIX]{x1D6F9}}_{2}^{\mathit{h}}|_{z=0}\approx (\unicode[STIX]{x1D707}_{0}Jr/2)\widehat{u}_{2,r}|_{z=0}.\end{eqnarray}$$

Analytical solutions can be obtained in terms of a series of harmonic functions. The calculation for $\widehat{\unicode[STIX]{x1D6F9}}_{2}^{\mathit{h}}$ is quite tedious and is not presented here for brevity. The computation for $\widehat{\unicode[STIX]{x1D6F9}}_{2}^{\mathit{v}}$ is easier and we can propose

(2.60) $$\begin{eqnarray}\widehat{\unicode[STIX]{x1D6F9}}_{2}^{\mathit{v}}=\left[C_{0}\left(\frac{r}{R}\right)^{m}+\mathop{\sum }_{j=1}^{+\infty }C_{j}\text{I}_{m}(j\unicode[STIX]{x03C0}r/H_{2})\cos (j\unicode[STIX]{x03C0}z/H_{2})\right]\text{e}^{\text{i}m\unicode[STIX]{x1D703}},\end{eqnarray}$$

as a candidate harmonic solution. Here $\text{I}_{m}$ is a modified Bessel function. The coefficients $C_{0},C_{1},C_{2},\ldots$ are calculated by projecting the inhomogeneous boundary condition $\unicode[STIX]{x2202}_{r}\widehat{\unicode[STIX]{x1D6F9}}_{2}^{\mathit{v}}|_{r=R}=\widehat{u}_{\unicode[STIX]{x1D703},2}|_{r=R}B_{z}$ onto the set of the basis functions $\{\cos (j\unicode[STIX]{x03C0}z/H_{2}),j\in \mathbb{N}\}$ using the orthogonality relation

(2.61) $$\begin{eqnarray}\int _{-H_{2}}^{0}\cos \frac{j\unicode[STIX]{x03C0}z}{H_{2}}\cos \frac{j^{\prime }\unicode[STIX]{x03C0}z}{H_{2}}=\unicode[STIX]{x1D6FF}_{j,j^{\prime }}\left\{\begin{array}{@{}ll@{}}H_{2},\quad & j=0\\ H_{2}/2,\quad & j\neq 0.\end{array}\right.\end{eqnarray}$$

This yields

(2.62) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}C_{0}=-\text{i}{\displaystyle \frac{\unicode[STIX]{x1D714}ARB_{z}}{k^{2}H_{2}}}\text{J}_{m}(kR),\\ C_{j}=-\text{i}{\displaystyle \frac{\unicode[STIX]{x1D714}ARB_{z}}{kj\unicode[STIX]{x03C0}}}{\displaystyle \frac{4m\text{J}_{m}(kR)}{(kR)(I_{m+1}(j\unicode[STIX]{x03C0}R/H_{2})+I_{m-1}(j\unicode[STIX]{x03C0}R/H_{2}))}}{\displaystyle \frac{k^{2}}{k^{2}+(j\unicode[STIX]{x03C0}/H_{2})^{2}}}\sinh (kH_{2}).\end{array}\right\}\end{eqnarray}$$

With this information we know $\widehat{\unicode[STIX]{x1D6F9}}_{2}^{\mathit{v}}$ and so we can calculate the integral ${\mathcal{Q}}_{\mathit{vv}}$ ignoring contributions from the top layer. This integral can be rewritten as

(2.63) $$\begin{eqnarray}{\mathcal{Q}}_{\mathit{vv}}\approx -\unicode[STIX]{x1D70E}_{2}\int _{V_{2}}\Vert \widehat{\boldsymbol{u}}_{2}^{\ast }\times B_{z}\boldsymbol{e}_{z}\Vert ^{2}\,\text{d}V+\int _{V_{2}}\widehat{\boldsymbol{u}}_{2}^{\ast }\boldsymbol{\cdot }\left(-\unicode[STIX]{x1D70E}_{2}\unicode[STIX]{x1D735}\widehat{\unicode[STIX]{x1D6F9}}_{2}^{\mathit{v}}\times B_{z}\boldsymbol{e}_{z}\right)\text{d}V.\end{eqnarray}$$

In containers with perfectly conducting side walls, we would have $\widehat{\unicode[STIX]{x1D6F9}}_{2}^{\mathit{v}}=0$ and only the first term would then contribute. Here the side walls are insulating and therefore $\widehat{\unicode[STIX]{x1D6F9}}_{2}^{\mathit{v}}\neq 0$ . This second term can drastically reduce the magnetic damping as explained by Sreenivasan, Davidson & Etay (Reference Sreenivasan, Davidson and Etay2005) and should certainly not be ignored. Analytical evaluation of all the integrals finally yields the following formula

(2.64) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}_{\mathit{vv}} & = & \displaystyle {\mathcal{Q}}_{\mathit{vv}}/2{\mathcal{K}}\nonumber\\ \displaystyle & {\approx} & \displaystyle \frac{\unicode[STIX]{x1D70E}_{2}B_{z}^{2}}{{\displaystyle \frac{\unicode[STIX]{x1D70C}_{1}}{\tanh (kH_{1})}}+{\displaystyle \frac{\unicode[STIX]{x1D70C}_{2}}{\tanh (kH_{2})}}}\left\{-\frac{1}{4}\left({\displaystyle \frac{1}{\tanh (kH_{2})}}+{\displaystyle \frac{kH_{2}}{\sinh ^{2}(kH_{2})}}\right)+\left(\frac{1}{k^{2}R^{2}-m^{2}}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \times \left.\left[\frac{m}{kH_{2}}+4m^{2}\mathop{\sum }_{j=1}^{+\infty }\left(\frac{1}{j\unicode[STIX]{x03C0}}\right)\!\frac{k^{3}}{R(k^{2}+(j\unicode[STIX]{x03C0}/H_{2})^{2})^{2}}\frac{\text{I}_{m}(j\unicode[STIX]{x03C0}R/H_{2})}{I_{m+1}(j\unicode[STIX]{x03C0}R/H_{2})+I_{m-1}(j\unicode[STIX]{x03C0}R/H_{2})}\!\right]\!\right\}.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

In the infinite sum over $j$ , we observe a fast exponential decay of the terms and the shallower the cell, the faster the decay. All the computations reported in the paper are done by truncating the sum over $j$ to the first 50 terms. Below we successfully compare this formula to the numerical finite-difference approach. In the shallow limit, $kH_{i}\ll 1$ , a single term is dominant and yields

(2.65) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{\mathit{vv},\mathit{shallow}}\approx \frac{\unicode[STIX]{x1D70E}_{2}B_{z}^{2}}{\unicode[STIX]{x1D70C}_{1}\left({\displaystyle \frac{H_{2}}{H_{1}}}\right)+\unicode[STIX]{x1D70C}_{2}}\left(-\frac{1}{2}+\frac{m}{k^{2}R^{2}-m^{2}}\right)\end{eqnarray}$$

for the magnetic damping caused by the vertical magnetic field. No significant simplification occurs in (2.64) in the deep limit $kH_{i}\gg 1$ because the cosine basis used in the expansion is not adapted to the infinitely deep case.

2.6.3 Viscous corrections: damping $\unicode[STIX]{x1D706}_{\mathit{visc}}$ , frequency shift $\unicode[STIX]{x1D6FF}_{\mathit{visc}}$

We finally include the viscous effects. To leading order we can ignore all magnetohydrodynamic effects. Viscosity will act on the waves just as in a purely hydrodynamic set-up. The earliest calculations of viscous damping rates of waves in cylinders filled with one fluid go back to Case & Parkinson (Reference Case and Parkinson1957), but no explicit formula is available for the case of two fluid layers. In appendix D we detail our calculation of the viscous effects on gravity waves. First, we find the leading-order spatial structure of the viscous boundary layers (Stokes layers). Near the solid walls, these boundary layers resemble those we can find in a cylinder filled with one liquid. Different is that, near the interface between both liquids, a viscous boundary layer is required to ensure the continuity of the flow field and that of the tangential viscous stress. By extending the method of Case & Parkinson (Reference Case and Parkinson1957) to the two-layer case, we then calculate the viscous damping rates from a perturbative expansion on the mechanical energy balance. A second, alternative approach, passes by a solvability condition. It confirms the formula for the viscous damping rate and also provides access to the frequency shift caused by viscosity. Our formulas are only valid for a freely sliding contact line making an angle of $90^{\circ }$ with the solid wall and we ignore contact line hysteresis (see Viola & Gallaire Reference Viola and Gallaire2018) and suppose an unpolluted interface. We find

(2.66a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{\mathit{visc}}=\unicode[STIX]{x1D706}_{\mathit{visc}}^{s}+\unicode[STIX]{x1D706}_{\mathit{visc}}^{i},\quad \unicode[STIX]{x1D6FF}_{\mathit{visc}}=\text{Sgn}(\unicode[STIX]{x1D714})\unicode[STIX]{x1D706}_{\mathit{visc}},\end{eqnarray}$$

with

(2.66c ) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{\mathit{visc}}^{s}=-\frac{1}{2\sqrt{2}}\mathop{\sum }_{i=1,2}\frac{\sqrt{{\displaystyle \frac{\unicode[STIX]{x1D708}_{i}|\unicode[STIX]{x1D714}|}{R^{2}}}}\unicode[STIX]{x1D70C}_{i}\left[{\displaystyle \frac{k(R-H_{i})}{\sinh ^{2}(kH_{i})}}+\left({\displaystyle \frac{k^{2}R^{2}+m^{2}}{k^{2}R^{2}-m^{2}}}\right){\displaystyle \frac{1}{\tanh (kH_{i})}}\right]}{\left[{\displaystyle \frac{\unicode[STIX]{x1D70C}_{1}}{\tanh (kH_{1})}}+{\displaystyle \frac{\unicode[STIX]{x1D70C}_{2}}{\tanh (kH_{2})}}\right]}\end{eqnarray}$$

and

(2.66d ) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{\mathit{visc}}^{i}=-\frac{\sqrt{|\unicode[STIX]{x1D714}|}k}{2\sqrt{2}}\frac{\left[{\displaystyle \frac{1}{\tanh (kH_{1})}}+{\displaystyle \frac{1}{\tanh (kH_{2})}}\right]^{2}}{\left[{\displaystyle \frac{\unicode[STIX]{x1D70C}_{1}}{\tanh (kH_{1})}}+{\displaystyle \frac{\unicode[STIX]{x1D70C}_{2}}{\tanh (kH_{2})}}\right]\left({\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{1}\sqrt{\unicode[STIX]{x1D708}_{1}}}}+{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{2}\sqrt{\unicode[STIX]{x1D708}_{2}}}}\right)}.\end{eqnarray}$$

The term $\unicode[STIX]{x1D706}_{\mathit{visc}}^{s}$ relates to viscous dissipation caused in boundary layers near solid walls, whereas the term $\unicode[STIX]{x1D706}_{\mathit{visc}}^{i}$ is caused by the boundary layers at the interface. Both terms can have comparable magnitudes in liquids with similar densities and kinematic viscosities. A much weaker higher-order damping term (often referred to as the interior damping term) is here ignored, but is calculated in appendix D. In the limit of one fluid layer, $\unicode[STIX]{x1D70C}_{1}\rightarrow 0$ , we recover the known damping rates (Case & Parkinson Reference Case and Parkinson1957; Ibrahim Reference Ibrahim2005). In the shallow limit $kH_{i}\ll 1$ , we can derive the asymptotic form

(2.67a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}_{\mathit{visc},\mathit{shallow}}^{s} & {\approx} & \displaystyle -\frac{1}{2\sqrt{2}}\frac{\unicode[STIX]{x1D70C}_{1}H_{1}^{-2}\sqrt{\unicode[STIX]{x1D708}_{1}}+\unicode[STIX]{x1D70C}_{2}H_{2}^{-2}\sqrt{\unicode[STIX]{x1D708}_{2}}}{\unicode[STIX]{x1D70C}_{1}H_{1}^{-1}+\unicode[STIX]{x1D70C}_{2}H_{2}^{-1}}\sqrt{|\unicode[STIX]{x1D714}_{\mathit{shallow}}|},\end{eqnarray}$$

(2.67b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}_{\mathit{visc},\mathit{shallow}}^{i} & {\approx} & \displaystyle -\frac{1}{2\sqrt{2}}\frac{(H_{1}^{-1}+H_{2}^{-1})^{2}}{(\unicode[STIX]{x1D70C}_{1}H_{1}^{-1}+\unicode[STIX]{x1D70C}_{2}H_{2}^{-1})\left({\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{1}\sqrt{\unicode[STIX]{x1D708}_{1}}}}+{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{2}\sqrt{\unicode[STIX]{x1D708}_{2}}}}\right)}\sqrt{|\unicode[STIX]{x1D714}_{\mathit{shallow}}|},\end{eqnarray}$$

and in the deep limit, $kH_{i}\gg 1$ , we can use

(2.68a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}_{\mathit{visc},\mathit{deep}}^{s} & {\approx} & \displaystyle -\frac{1}{2\sqrt{2}}\left({\displaystyle \frac{k^{2}R^{2}+m^{2}}{k^{2}R^{2}-m^{2}}}\right)\frac{1}{R}\frac{\unicode[STIX]{x1D70C}_{1}\sqrt{\unicode[STIX]{x1D708}_{1}}+\unicode[STIX]{x1D70C}_{2}\sqrt{\unicode[STIX]{x1D708}_{2}}}{\unicode[STIX]{x1D70C}_{1}+\unicode[STIX]{x1D70C}_{2}}\sqrt{|\unicode[STIX]{x1D714}_{\mathit{deep }}|},\end{eqnarray}$$

(2.68b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}_{\mathit{visc},\mathit{deep}}^{i} & {\approx} & \displaystyle -\sqrt{2}k\left[(\unicode[STIX]{x1D70C}_{1}+\unicode[STIX]{x1D70C}_{2})\left({\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{1}\sqrt{\unicode[STIX]{x1D708}_{1}}}}+{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{2}\sqrt{\unicode[STIX]{x1D708}_{2}}}}\right)\right]^{-1}\sqrt{|\unicode[STIX]{x1D714}_{\mathit{deep }}|}\end{eqnarray}$$

as an approximation.

4.1 Theoretical predictions

Due to its small radius, limited $J$ and relatively high $B_{z}$ , the horizontal magnetic field is negligible: $B_{\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D707}_{0}Jr/2\ll B_{z}$ . As a result, we obtain a very good approximation for the growth rate and frequency shift with

(4.3a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D706}\approx \unicode[STIX]{x1D706}_{\mathit{v}}+\unicode[STIX]{x1D706}_{\mathit{vv}}+\unicode[STIX]{x1D706}_{\mathit{visc}},\quad \unicode[STIX]{x1D6FF}\approx \unicode[STIX]{x1D6FF}_{\mathit{visc}}.\end{eqnarray}$$

Table 1 provides numerical values obtained from our theory (in SI units) for the wavenumbers $k$ and the frequencies $\unicode[STIX]{x1D714}$ of several waves $(m,n)$ together with $\unicode[STIX]{x1D706}_{\mathit{v}}/JB_{z}$ , $\unicode[STIX]{x1D706}_{\mathit{vv}}/B_{z}^{2}$ , $\unicode[STIX]{x1D706}_{\mathit{visc}}(=\unicode[STIX]{x1D6FF}_{visc})$ and the ratio $\unicode[STIX]{x1D706}_{\mathit{visc}}^{i}/\unicode[STIX]{x1D706}_{\mathit{visc}}$ that characterizes how much of the damping is caused by dissipation near the interface. Some quantities in this table depend on surface tension, the superscript $(a)$ corresponds to the case without surface tension, the superscript $(b)$ to the case with surface tension. By inspecting the columns $\unicode[STIX]{x1D706}_{\mathit{v}}/JB_{z}$ we see that the large-scale sloshing mode $(m,n)=(1,1)$ is not the only wave that is unstable without dissipation. In absence of surface tension, the higher $m$ modes are even more unstable, but with surface tension, the mode $(1,1)$ is the most unstable. Waves with smaller radial structures, i.e. $n\geqslant 2$ will always be less destabilized. The magnetic damping $\unicode[STIX]{x1D706}_{\mathit{vv}}/B_{z}^{2}$ is independent of the surface tension and increases with $m$ and $n$ . Viscous damping $\unicode[STIX]{x1D706}_{\mathit{visc}}$ increases with $m$ and $n$ and surface tension always increases the viscous dissipation, mainly because the frequency of the wave increases. In the last column of table 1, we show the percentage of damping that is caused by viscous dissipation near the interface. For the large-scale waves with $n=1$ , dissipation near the interface and near the solid walls are grossly equivalent, $\unicode[STIX]{x1D706}_{\mathit{visc}}^{i}\approx \unicode[STIX]{x1D706}_{\mathit{visc}}^{s}$ , but for small-scale waves with $n=2$ , dissipation in interfacial layers clearly dominates.

In figure 5, we show curves of marginal instability in the $J$ – $B_{z}$ plane calculated using the data of table 1. The first mode that becomes unstable is the large-scale sloshing wave, $(m,n)=(1,1)$ , followed by the other waves $(m,1)$ with $m\geqslant 2$ . In absence of surface tension, the marginal stability curves for the different waves are very close to each other: many different waves quickly destabilize, one after the other. When the surface tension is active, the $(1,1)$ mode will be dominant in the shown part of parameter space. The horizontal dashed line in figure 5(a), shows the path that is followed in figure 9(a), where we show the growth rate $\unicode[STIX]{x1D706}$ of the (most unstable) $(1,1)$ -mode as a function of $J$ with fixed $B_{z}=16~\text{mT}$ . The theoretical line $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D706}_{\mathit{v}}+\unicode[STIX]{x1D706}_{\mathit{vv}}+\unicode[STIX]{x1D706}_{\mathit{visc}}$ , grows linearly in $J$ because $\unicode[STIX]{x1D706}_{\mathit{v}}\sim |JB_{z}|$ . This figure will be further discussed below.

Figure 5.

Small cell. Marginal stability curves in the $J$ – $B_{z}$ plane for metal pad roll instability of waves $(m,n)$ without surface tension (a) or with surface tension (b). In the regions shaded in grey the cell is stable. Surface tension has stabilizing effects. The fundamental wave $(m,n)=(1,1)$ always becomes unstable first.

In figure 6(a), we show the currents $\widehat{\boldsymbol{j}}_{i}$ , $\widehat{\boldsymbol{{\mathcal{J}}}}_{i}^{\mathit{v}}$ that are associated with the fundamental wave $(1,1)$ for an equilibrium electrolysis current $\boldsymbol{J}$ flowing from the top to the bottom. The deformed interface is represented by the black line. The current lines $\widehat{\boldsymbol{j}}_{i}$ , due to the interface motion, loop in a symmetrical way in both fluid layers. This is very different from what we usually see in shallow cells. As expected, the quasi-static correction $\widehat{\boldsymbol{{\mathcal{J}}}}_{i}^{\mathit{v}}$ is mainly confined inside the aluminium, where induction by the flow takes place. In figure 6(b), we show the normalized power densities that appear in the integrands of ${\mathcal{P}}_{\mathit{v}}$ , ${\mathcal{Q}}_{\mathit{vv}}$ . This allows us to visually locate where the instability is powered and where Joule damping takes place. We see that both the cryolite and the aluminium participate in a symmetric way in the destabilization through ${\mathcal{P}}_{\mathit{v}}$ . The density of ${\mathcal{Q}}_{\mathit{vv}}$ is negative in the bottom aluminium layer, where the induced current $\widehat{\boldsymbol{{\mathcal{J}}}}_{i}^{\mathit{v}}$ is ohmically dissipated.

Figure 6.

The fundamental wave $(1,1)$ in the small cell (suggested by the thick black line) causes current perturbations. In (a) we show streamlines of the current perturbation $\widehat{\boldsymbol{j}}_{i}$ related to the interface motion (magneto-static part) and $\widehat{\boldsymbol{{\mathcal{J}}}}_{i}^{\mathit{v}}$ related to induction by the flow in the bottom layer. In (b), we show normalized plots of the power densities that occur in the integrands of ${\mathcal{P}}_{\mathit{v}}$ and ${\mathcal{Q}}_{\mathit{vv}}$ .

Table 2.

Comparisons between the theoretical formulas in (2.39) and the deep and shallow limits for inviscid frequencies $\unicode[STIX]{x1D714}$ (in $\text{s}^{-1}$ ), relative growth rates $\unicode[STIX]{x1D706}_{\mathit{v}}/JB_{z}$ (in $\text{s}~\text{m}^{2}~\text{kg}^{-1}$ ), relative damping term $\unicode[STIX]{x1D706}_{\mathit{vv}}/B_{z}^{2}$ (in $\text{s}^{-1}~\text{T}^{-2}$ ) and viscous damping rates $\unicode[STIX]{x1D706}_{\mathit{visc}}$ (in $\text{s}^{-1}$ ). Surface tension is ignored. The deep limit is an excellent approximation in the small cell.

In table 2, we show that the deep-limit formulas apply very well to this cell. We focus our attention on the two waves $(m,n)=(1,1),(2,1)$ . Shallow-limit formulas (2.43) are also shown, only to illustrate that they are not adapted here: no existing shallow theory would be able to accurately model the metal pad roll instability in this small cell.

Figure 7.

Variation of the relative growth rate $\unicode[STIX]{x1D706}_{\mathit{v}}/JB_{z}$ (a) and magnetic damping term $\unicode[STIX]{x1D706}_{\mathit{vv}}/B_{z}^{2}$ (b) when $\unicode[STIX]{x1D70E}_{1}$ is artificially increased. From $\unicode[STIX]{x1D70E}_{1}=250~\text{S}~\text{m}^{-1}$ to $\unicode[STIX]{x1D70E}_{1}=\unicode[STIX]{x1D70E}_{2}/100$ there is only a weak modification of the growth rate and magnetic damping. At $\unicode[STIX]{x1D70E}_{1}=\unicode[STIX]{x1D70E}_{2}$ , the instability vanishes: a strong difference in conductivity is essential for the metal pad roll instability.

Table 3.

The relative magnetic damping rates $\unicode[STIX]{x1D706}_{\mathit{vv}}/B_{z}^{2}$ (in $\text{s}^{-1}~\text{T}^{-2}$ ) are calculated numerically for different waves $(m,n)$ using a finite-difference solver with variable spatial resolution $M=100$ , 200, 400, 800 and compared to the theoretical estimate (2.64) obtained by assuming $\unicode[STIX]{x1D70E}_{1}\ll \unicode[STIX]{x1D70E}_{2}$ . The numerical values converge towards the theoretical estimate.

In table 3 we compare the theoretical estimation of $\unicode[STIX]{x1D706}_{\mathit{vv}}/B_{z}^{2}$ given by formula (2.64) to numerically calculated values, obtained using numerical quadrature and the field profiles $\widehat{\unicode[STIX]{x1D6F9}}_{i}^{\mathit{v}}$ that result from the finite-difference code. Increasing the spatial resolution $M=100$ , 200, 400, 800 in the code, we observe that the values of the damping rates converge towards the theoretical prediction of (2.64). This validates our finite difference solver, and also suggests that our theoretical formula (2.64) for the magnetic damping rate is correct.

In figure 7, we inspect how the rates $\unicode[STIX]{x1D706}_{\mathit{v}}/JB_{z}$ and $\unicode[STIX]{x1D706}_{\mathit{vv}}/B_{z}^{2}$ of the fundamental wave $(m,n)=(1,1)$ may vary with the electrical conductivity of the top layer. In reality, we obviously cannot choose the conductivity of the cryolite, but we want to know how important it is to use the exact value for $\unicode[STIX]{x1D70E}_{1}$ . We investigate the range $\unicode[STIX]{x1D70E}_{1}\in [250~\text{S}~\text{m}^{-1},\unicode[STIX]{x1D70E}_{2}]$ . We use the numerical approach to calculate $\unicode[STIX]{x1D706}_{\mathit{vv}}/B_{z}^{2}$ , since the theoretical approximation (2.64) is not valid when $\unicode[STIX]{x1D70E}_{1}\not \ll \unicode[STIX]{x1D70E}_{2}$ . For $\unicode[STIX]{x1D70E}_{1}$ in the range $[250~\text{S}~\text{m}^{-1},\unicode[STIX]{x1D70E}_{2}/100]$ , we observe little change in both $\unicode[STIX]{x1D706}_{\mathit{v}}/JB_{z}$ and $\unicode[STIX]{x1D706}_{\mathit{vv}}/B_{z}^{2}$ . Increasing further $\unicode[STIX]{x1D70E}_{1}$ does lower the growth rate $\unicode[STIX]{x1D706}_{\mathit{v}}$ significantly and when $\unicode[STIX]{x1D70E}_{1}=\unicode[STIX]{x1D70E}_{2}$ , we ultimately get $\unicode[STIX]{x1D706}_{\mathit{v}}=0$ . Magnetic damping only gets stronger with increasing $\unicode[STIX]{x1D70E}_{1}$ , since ohmic dissipation may then also occur in the top layer. We conclude here that the low conductivity of the electrolyte $\unicode[STIX]{x1D70E}_{1}$ can be significantly increased without quantitively modifying the instability, but that a true jump in conductivity is crucial for the metal pad roll instability.

4.2 Numerically solved problem

The numerical methods used in SFEMaNS and OpenFOAM were previously described in Guermond et al. (Reference Guermond, Laguerre, Léorat and Nore2007, Reference Guermond, Laguerre, Léorat and Nore2009), Weber et al. (Reference Weber, Galindo, Stefani, Weier and Wondrak2013, Reference Weber, Galindo, Stefani and Weier2014, Reference Weber, Beckstein, Galindo, Herreman, Nore, Stefani and Weier2017a ,Reference Weber, Beckstein, Herreman, Horstmann, Nore, Stefani and Weier b , Reference Weber, Beckstein, Galindo, Starace and Weier2018), Nore et al. (Reference Nore, Quiroz, Cappanera and Guermond2016), Ashour et al. (Reference Ashour, Kelley, Salas, Starace, Weber and Weier2018), Cappanera et al. (Reference Cappanera, Guermond, Herreman and Nore2018) and Horstmann et al. (Reference Horstmann, Weber and Weier2018). SFEMaNS (for spectral finite element code for Maxwell and Navier–Stokes equations) is a hybrid finite element–spectral solver that can handle axisymmetric fluid domains, OpenFOAM is a public library that uses a finite volume method. Both codes are used here to solve the nonlinear momentum equation and mass balance

(4.4) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D70C}(\unicode[STIX]{x2202}_{t}\boldsymbol{U}+(\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{U})=-\unicode[STIX]{x1D735}P+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}\unicode[STIX]{x1D708}(\unicode[STIX]{x1D735}\boldsymbol{U}+\unicode[STIX]{x1D735}\boldsymbol{U}^{\text{T}}))+\boldsymbol{J}\times \boldsymbol{B},\\ \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{U}=0,\end{array}\right\}\end{eqnarray}$$

in a cylindrical fluid domain. The fields $\boldsymbol{U}$ , $P$ , $\boldsymbol{B}$ , $\boldsymbol{J}=\unicode[STIX]{x1D707}_{0}^{-1}\unicode[STIX]{x1D735}\times \boldsymbol{B}$ correspond here to the total flow, pressure, magnetic induction and current density. To deal with the multi-phase character, the material parameters $\unicode[STIX]{x1D70C},\unicode[STIX]{x1D708}$ are considered as space- and time-dependent functions that may be reconstructed using a scalar field $f\in [0,1]$ that sharply varies from 0 to 1 on the interface. This function allows us to reconstruct the material parameters of the different liquid phases, e.g. $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{1}f+\unicode[STIX]{x1D70C}_{2}(1-f)$ for the density in the two-phase system. By materially advecting the indicator function $f$ , that is, solving

(4.5) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}f+(\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735})f\approx 0,\end{eqnarray}$$

we keep track of the different phases as time evolves. This is done in SFEMaNS by using a level set method. OpenFOAM uses a volume of fluid method. Surface tension may be dealt with in both codes, but SFEMaNS nor OpenFOAM use particular numerical models for the contact line motion. As a result no realistic contact angle other than $90^{\circ }$ can be used. Therefore, we limit all the direct numerical simulations to the case without surface tension: $\unicode[STIX]{x1D6FE}_{1|2}=0$ . On the side wall and top and bottom surfaces, we impose the no-slip condition.

(4.6) $$\begin{eqnarray}\boldsymbol{U}=\mathbf{0}|_{\unicode[STIX]{x1D6F4}_{1}\cup \unicode[STIX]{x1D6F4}_{2}}.\end{eqnarray}$$

This no-slip condition is imposed in a weak way and, as a result, the ‘contact line’ (region where $f$ changes from 0 to 1 near $r=R$ ) can still freely slide on the cylindrical wall. The magnetic problem is solved differently in both codes. SFEMaNS solves the full induction equation and Gauss’ law

(4.7) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x2202}_{t}\boldsymbol{B}=\unicode[STIX]{x1D735}\times (\boldsymbol{U}\times \boldsymbol{B})-\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D70E}^{-1}\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D707}_{0}^{-1}\boldsymbol{B})),\\ 0=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B},\end{array}\right\}\end{eqnarray}$$

together with the boundary conditions

(4.8a ) $$\begin{eqnarray}\displaystyle B_{z}|_{r=R} & = & \displaystyle B_{z},\quad B_{\unicode[STIX]{x1D703}}|_{r=R}=\frac{\unicode[STIX]{x1D707}_{0}JR}{2},\end{eqnarray}$$

(4.8b ) $$\begin{eqnarray}\displaystyle \text{J}_{r}|_{z=H_{1}} & = & \displaystyle 0,\quad \text{J}_{\unicode[STIX]{x1D703}}|_{z=H_{1}}=0,\end{eqnarray}$$

(4.8c ) $$\begin{eqnarray}\displaystyle B_{r}|_{z=-H_{2}} & = & \displaystyle 0,\quad B_{\unicode[STIX]{x1D703}}|_{z=-H_{2}}=\frac{\unicode[STIX]{x1D707}_{0}Jr}{2},\end{eqnarray}$$

$\unicode[STIX]{x1D70E}^{-1}$

$f$

$\boldsymbol{E}\approx -\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}$

(4.9) $$\begin{eqnarray}\boldsymbol{J}\approx \unicode[STIX]{x1D70E}\left(-\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}+\boldsymbol{U}\times \left(\frac{\unicode[STIX]{x1D707}_{0}Jr}{2}\boldsymbol{e}_{\unicode[STIX]{x1D703}}+B_{z}\boldsymbol{e}_{z}\right)\right).\end{eqnarray}$$

Combining this with $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{J}=0$ , we find a Poisson problem for the total electrical potential $\unicode[STIX]{x1D6F7}$ , that can be solved with the boundary conditions

(4.10a-c ) $$\begin{eqnarray}\text{J}_{r}=0|_{r=R},\quad \unicode[STIX]{x1D6F7}|_{z=H_{1}}=0,\quad \text{J}_{z}|_{z=-H_{2}}=J.\end{eqnarray}$$

Once the current density is known, the total magnetic field $\boldsymbol{B}$ can be calculated using the Biot–Savart integral. In this $\boldsymbol{B}$ , the base state magnetic field $(\unicode[STIX]{x1D707}_{0}Jr/2)\boldsymbol{e}_{\unicode[STIX]{x1D703}}+B_{z}\boldsymbol{e}_{z}$ always dominates in the quasi-static limit, but the small deviations in the magnetic field still need to be calculated in order to have a consistent approximation of the Lorentz force $\boldsymbol{J}\times \boldsymbol{B}$ (this correction was crucial for the Tayler instability (Weber et al. Reference Weber, Galindo, Stefani, Weier and Wondrak2013, Reference Weber, Galindo, Stefani and Weier2014) but is likely less important here).

In both codes, we have varied the grid size to assess the numerical convergence. We introduce the non-dimensional measure $h=\unicode[STIX]{x1D6FF}x/R$ to characterize the grid size. In SFEMaNS the mesh size $\unicode[STIX]{x1D6FF}x$ is non-uniform in the plane $(r,z)$ and the mesh cells are refined near the interface. Therefore we provide intervals $h=(h_{min}\rightarrow h_{max})$ to characterize the SFEMaNS grid. All the simulations reported in the paper have been done with 30 real Fourier modes in the azimuthal direction. In SFEMaNS, we relax the conductivity jump by taking $\unicode[STIX]{x1D70E}_{1}=\unicode[STIX]{x1D70E}_{2}/100$ . This stabilizes the numerical calculations and, as explained above, should not modify the linear instability mechanism. In OpenFOAM only purely orthogonal cells are used with $h=0.05$ or smaller. The cells are flattened towards the interface with an aspect ratio of 5 at most. Interface compression ensures that the interface is typically smeared over 2–3 cells. Information on the initialization of both codes is provided below.

Table 4.

Using both SFEMaNS and OpenFOAM, we simulate the decay of large-scale gravity waves $(m,n)=(1,1),(2,1)$ with or without an imposed magnetic field $B_{z}=0$ , 10, 25, 50 mT and using different grids with typical mesh sizes $h=\unicode[STIX]{x1D6FF}x/R$ . We measure the frequency $\unicode[STIX]{x1D714}+\unicode[STIX]{x1D6FF}_{\mathit{visc}}$ ( $\text{s}^{-1}$ ), the viscous damping rate $\unicode[STIX]{x1D706}_{\mathit{visc}}$ ( $\text{s}^{-1}$ ) and the magnetic damping rates $\unicode[STIX]{x1D706}_{\mathit{vv}}$ ( $\text{s}^{-1}$ ) of the waves and we compare with the theoretical estimates provided by the perturbative analysis.

4.3 Numerical results and comparison with theory

4.3.2 Viscous and magnetically damped gravity waves

In a second series of tests, we measure the damping rates of gravity waves in the presence of non-zero external magnetic fields $B_{z}\neq 0$ , but without an electrolysis current ( $J=0$ ). We perform simulations for increasing magnetic field $B_{z}=10$ , 25, 50 mT. Denoting by $\unicode[STIX]{x1D706}$ the growth rate of the waves in the numerical simulations using a mesh with cell size $h$ and $\unicode[STIX]{x1D706}_{\mathit{visc}}$ the viscous growth rate measured for this mesh, we then calculate the numerical values $\unicode[STIX]{x1D706}_{\mathit{vv}}=\unicode[STIX]{x1D706}-\unicode[STIX]{x1D706}_{\mathit{visc}}$ in order to filter the code- and grid-dependent viscous damping that was considered in the previous section. Table 4 gathers the numerical measures for $\unicode[STIX]{x1D706}_{\mathit{vv}}$ . These measures are compared with the theoretical estimates for waves $(1,1)$ and $(2,1)$ . The agreement with the theoretical value is very good on sufficiently fine grids and excellent for the small values of $B_{z}=10$ , 25 mT. For the largest value $B_{z}=50~\text{mT}$ there is a small mismatch, which we suspect is due to the fact that we leave the asymptotic regime of validity of the theoretical model. We notice that the magnetic damping rates $-\unicode[STIX]{x1D706}_{\mathit{vv}}$ are slightly lower in OpenFOAM than in the simulations done with SFEMaNS.

Figure 8.

Visualization of the interface between the cryolite and the aluminium in the numerical simulations. (a,b) Eight snapshots over one period $T$ in the saturated regime for $J=-7795~\text{A}~\text{m}^{-2}$ , $B_{z}=16~\text{mT}$ . (c) Far from the stability threshold, the interface violently destabilizes and develops large deformations, $J=-20\,000~\text{A}~\text{m}^{-2}$ , $B_{z}=16~\text{mT}$ .

Figure 9.

(a) Comparison of theoretical and numerical growth rates for the large-scale sloshing wave $(m,n)=(1,1)$ as a function of $|J|$ . (b) Saturation amplitudes of rotating wave as a function of $|J|$ . Small cell configuration with $B_{z}=16~\text{mT}$ and ignoring surface tension.

4.3.3 Metal pad roll instability

In a third series of tests, we perform simulations of the metal pad roll instability with both $J$ and $B_{z}$ active. We fix the magnetic field strength to $B_{z}=16~\text{mT}$ which is used in Steiner (Reference Steiner2009). We increase the current density in the interval $J\in -[7795,25\,984]~\text{A}~\text{m}^{-2}$ which corresponds to electrolysis currents in the range $I=[30,100]~\text{A}$ . The current runs downwards through the cell. From the theory, we estimate that the metal pad roll instability should occur when $|J|$ is above the threshold $\text{J}_{c}=4430~\text{A}~\text{m}^{-2}$ .

Figure 8 shows snapshots from simulations at $|J|=7795~\text{A}~\text{m}^{-2}>\text{J}_{c}$ . The metal pad roll instability is clearly visible. We observe a saturated wave similar in shape to the fundamental wave $(m,n)=(1,1)$ , and this wave rotates in the anti-clockwise direction $(m\unicode[STIX]{x1D714}<0)$ when seen from above. This behaviour corresponds to what we expect from the theory. In a more systematic series of computations we have measured the growth rate of the metal pad roll instability as a function of $J$ . The results are shown in figure 9(a). We observe an excellent agreement with theory. The fact that both SFEMaNS and OpenFOAM find straight lines with the same slope is an indication that both codes capture the destabilizing mechanism at the origin of $\unicode[STIX]{x1D706}_{\mathit{v}}$ very well. The growth rates from SFEMaNS are slightly beneath those from OpenFOAM, which is not unexpected owing to the small difference in viscous and magnetic damping rates we have observed before.

The numerical simulations further provide insights in the nonlinear regime, but since this regime is not modelled by the proposed theory we will not study it in detail here. In figure 9(b), we show a diagram that gathers the maximal amplitude $\unicode[STIX]{x1D702}_{max}$ of the saturated rotating gravity wave as a function of current density. This curve suggests a supercritical bifurcation. We propose a weakly nonlinear fit $\unicode[STIX]{x1D702}_{max}\approx 1.5\times 10^{-4}\sqrt{|J|-\text{J}_{c}}$ in MKS units (based on the OpenFOAM data). In non-dimensional form and using the theoretical formula for the growth rate $\unicode[STIX]{x1D706}$ , we rewrite this fit as follows:

(4.11) $$\begin{eqnarray}\text{If }\unicode[STIX]{x1D706}>0:\quad \frac{\unicode[STIX]{x1D702}_{max}}{R}\approx 2\sqrt{\frac{\unicode[STIX]{x1D706}}{|\unicode[STIX]{x1D714}|}}.\end{eqnarray}$$

It is likely that this non-dimensional relation is more generally applicable to estimate the weakly nonlinear saturation amplitude of the fundamental wave $(1,1)$ in deep reduction cells, but we certainly do not expect it to be valid for shallow cells.

The weakly nonlinear regime exists only in a small parameter range. Above $|J|\sim 15\,000~\text{A}~\text{m}^{-2}$ , violent secondary instabilities are observed with both numerical solvers. The interface undergoes strong deformations which ultimately lead to a blow up of the numerical calculations. To illustrate this phenomenon, we show in figure 8(c) the interface obtained with $J=-20\,000~\text{A}~\text{m}^{-2}$ using SFEMaNS, right before the blow-up. The precise origin of this secondary instability is unknown, but in this strongly nonlinear regime we do find the ‘explosive’ combination of

(4.12a,b ) $$\begin{eqnarray}\Vert \boldsymbol{J}_{NL}-J\boldsymbol{e}_{z}\Vert \sim 1,\quad \frac{|JB_{z}|}{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g}\geqslant 1.\end{eqnarray}$$

The total electrical current density ( $\boldsymbol{J}_{NL}$ ) may significantly differ from the equilibrium state and it is quite likely that the resulting Lorentz force can then no longer be equilibrated by pressure forces alone. If this local Lorentz force may, in addition, strongly overpower gravitational restoring forces, there is little reason for the cell to remain in a stably stratified state. Note that similar explosive behaviour was also observed in a liquid metal battery geometry (Weber et al. Reference Weber, Beckstein, Galindo, Herreman, Nore, Stefani and Weier2017a ,Reference Weber, Beckstein, Herreman, Horstmann, Nore, Stefani and Weier b ).

5.1 Large, industrial-scale cells

Industrial reduction cells are very large in lateral extent and each fluid layer is extremely shallow. As such, they are very different from the small reduction cell studied in the previous section (see the to-scale sketch of figure 10). We inspect how our theory behaves in these more realistic conditions and consider large cells, with dimensions

(5.1a-c ) $$\begin{eqnarray}R\in [0.5,5]~\text{m},\quad H_{1}=0.05~\text{m},\quad H_{2}=0.30~\text{m}.\end{eqnarray}$$

All the material parameters are taken the same as in (4.2). We ignore surface tension since the focus will be on long-wavelength waves.

Figure 10.

To-scale side view of the two reduction cell configurations studied in the article. (a) Small experimental cell with radius $R=0.035~\text{m}$ used for the validation of the theory. (b) Large industrial-scale cell (in this plot $R=1~\text{m}$ ).

Table 5.

Theoretical values for all the relevant quantities characterizing the stability and frequency of various waves $(m,n)$ in the large cell with radius $R=1~\text{m}$ . Interfacial tension is ignored, $\unicode[STIX]{x1D6FE}_{1|2}=0$ . Units are in the MKS system.

In table 5 we give the numerical values for all the relevant theoretical quantities determining the stability and frequency shifts for different waves $(m,n)$ in a cell of radius $R=1~\text{m}$ . In the column for $\unicode[STIX]{x1D706}_{\mathit{v}}/JB_{z}$ , we see that the fundamental wave $(1,1)$ is always the most unstable. Since the viscous dissipation is also the lowest for this wave, we expect this mode to be dominant in an unstable cell. The relative frequency shift $\unicode[STIX]{x1D6FF}_{\mathit{h}}/J^{2}$ can be either positive or negative and it decays as $m$ or $n$ increases. Although the numbers are small, this effect is no longer negligible for current densities of the order $J\sim 10^{4}~\text{A}~\text{m}^{-2}$ . A noteworthy detail is the positive value of $\unicode[STIX]{x1D706}_{\mathit{hh}}/J^{2}=8.06\times 10^{-14}~\text{s}^{-1}~\text{m}^{2}~\text{A}^{-2}$ for the wave $(1,1)$ : for large enough $J$ , this wave may be unstable even without an external magnetic field $B_{z}=0$ . However, it seems that this possibility is anecdotal: electrical currents would need to be so large that $\unicode[STIX]{x1D6FF}_{\mathit{h}}/\unicode[STIX]{x1D714}\gg 1$ so that our perturbation method would no longer be trustworthy. In table 6, we show that the shallow-limit formulas (see equations (2.43), (2.44), (2.65), (2.67)) apply very well to this shallow cell. Deep-limit formulas are not shown here but are not adapted.

Figure 11.

Growth rate $\unicode[STIX]{x1D706}$ and shifted frequency $\unicode[STIX]{x1D714}+\unicode[STIX]{x1D6FF}$ for various gravity waves $(m,n)$ in large cells with dimensions $R\in [0.5,2.5]~\text{m}$ , $H_{1}=0.05~\text{m}$ , $H_{2}=0.35~\text{m}$ . We take realistic values for the current density $J=1.5\times 10^{4}~\text{A}~\text{m}^{-2}$ and the vertical magnetic field $B_{z}=1~\text{mT}$ . The dashed lines in (b) correspond to the unperturbed gravity wave frequency $\unicode[STIX]{x1D714}$ ; the solid lines correspond to $\unicode[STIX]{x1D714}+\unicode[STIX]{x1D6FF}$ .

Table 6.

Comparison of the exact theoretical formulas with shallow limits for inviscid frequencies $\unicode[STIX]{x1D714}$ (in $\text{s}^{-1}$ ) and relative growth rates $\unicode[STIX]{x1D706}_{\mathit{v}}/JB_{z}$ (in $\text{s}^{-1}~\text{A}^{-1}~\text{m}^{2}~\text{T}^{-1}$ ), relative magnetic damping term $\unicode[STIX]{x1D706}_{\mathit{vv}}/B_{z}^{2}$ (in $\text{s}^{-1}~\text{T}^{-2}$ ), viscous damping term $\unicode[STIX]{x1D706}_{\mathit{visc}}$ (in $\text{s}^{-1}$ ) and relative frequency shift $\unicode[STIX]{x1D6FF}_{\mathit{h}}/J^{2}$ (in $\text{s}^{-1}~\text{A}^{-2}~\text{m}^{4}$ ). In the large cell under study, the shallow-limit formulas provide a very accurate description.

We now fix $J=1.5\times 10^{4}~\text{A}~\text{m}^{-2}$ and $B_{z}=1~\text{mT}$ , which we consider to be realistic values. Due to the low strength of the magnetic field, we can just neglect all the quasi-static MHD effects in the large cell as is done usually in the metal pad roll theory for shallow cells. The following approximation applies well:

(5.2a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D706}\approx \unicode[STIX]{x1D706}_{\mathit{v}}+\unicode[STIX]{x1D706}_{\mathit{visc}},\quad \unicode[STIX]{x1D6FF}\approx \unicode[STIX]{x1D6FF}_{\mathit{h}}+\unicode[STIX]{x1D6FF}_{\mathit{visc}}.\end{eqnarray}$$

In figure 11(a) we show how the theoretical growth rate $\unicode[STIX]{x1D706}$ increases almost linearly with $R$ in agreement with what was indicated by the shallow-limit formula (2.45). The large-scale sloshing mode $(1,1)$ is the most unstable mode over the entire $R$ -span and should be the selected mode in this cell. In figure 11(b), the dashed lines show the inviscid frequency $\unicode[STIX]{x1D714}$ and need to be compared with the solid lines that show the shifted frequency, $\unicode[STIX]{x1D714}+\unicode[STIX]{x1D6FF}$ . The frequency shift is clearly visible and mostly due to the contribution $\unicode[STIX]{x1D6FF}_{\mathit{h}}$ , i.e. it is caused by the Lorentz force interaction of the electrical currents $\widehat{\boldsymbol{j}}_{i}$ with the horizontal magnetic field $\unicode[STIX]{x1D707}_{0}Jr/2\,\boldsymbol{e}_{\unicode[STIX]{x1D703}}$ .

Figure 12.

Large cell geometry $R=1~\text{m}$ . (a) Streamlines of the current perturbation $\widehat{\boldsymbol{j}}_{i}$ for a downward base state electrolysis current density $\boldsymbol{J}$ . The interface deformation is suggested by the thicker black line. (b) Normalized power density for the fundamental wave $(1,1)$ .

In figure 12, we show the current perturbation $\widehat{\boldsymbol{j}}_{i}$ together with the power density in the integrand of ${\mathcal{P}}_{\mathit{v}}$ . Comparison with figure 6 allows us to highlight some crucial differences between deep and shallow cells. In the shallow cell, $\widehat{\boldsymbol{j}}_{i}$ is mainly vertical in the electrolyte and horizontal in the aluminium. The very low ${\mathcal{P}}_{\mathit{v}}$ -power density in the electrolyte layer indicates that no significant mechanical power transfer occurs in the thin electrolyte layer. This is all well known for shallow cells but noticeably different from what we had in the small (deep) cell of figure 6 (symmetrical $\widehat{\boldsymbol{j}}_{i}$ , mechanical power transfer in both the cryolite and the aluminium).

A common statement in metal pad roll theory is that cylindrical cells are always unstable, often expressed in terms of the critical $\unicode[STIX]{x1D6FD}$ number as $\unicode[STIX]{x1D6FD}_{c}=0$ ( $\unicode[STIX]{x1D6FD}$ defined in (1.1)). This statement obviously only holds when all dissipation is ignored, which is not the case in our theory. Focussing on large shallow cells, ignoring surface tension and magnetic damping (small), we can rearrange the equation $\unicode[STIX]{x1D706}_{\mathit{v},shallow}+\unicode[STIX]{x1D706}_{\mathit{visc},shallow}\approx 0$ for the onset of instability into an equation for the critical $\unicode[STIX]{x1D6FD}$ number in the presence of viscous damping,

(5.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}_{c,mn} & {\approx} & \displaystyle \frac{\unicode[STIX]{x03C0}}{2\sqrt{2}}\frac{\unicode[STIX]{x1D705}_{mn}^{3/2}(\unicode[STIX]{x1D705}_{mn}^{2}-m^{2})}{m}\nonumber\\ \displaystyle & & \displaystyle \times \,\frac{H_{1}^{-2}\left(2+{\displaystyle \frac{\unicode[STIX]{x1D70C}_{1}\sqrt{\unicode[STIX]{x1D708}_{1}}}{\unicode[STIX]{x1D70C}_{2}\sqrt{\unicode[STIX]{x1D708}_{2}}}}\right)+2(H_{1}H_{2})^{-1}+H_{2}^{-2}\left(2+{\displaystyle \frac{\unicode[STIX]{x1D70C}_{2}\sqrt{\unicode[STIX]{x1D708}_{2}}}{\unicode[STIX]{x1D70C}_{1}\sqrt{\unicode[STIX]{x1D708}_{1}}}}\right)}{((\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g)^{1/4}(\unicode[STIX]{x1D70C}_{1}H_{1}^{-1}+\unicode[STIX]{x1D70C}_{2}H_{2}^{-1})^{3/4}\left({\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{1}\sqrt{\unicode[STIX]{x1D708}_{1}}}}+{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{2}\sqrt{\unicode[STIX]{x1D708}_{2}}}}\right)}R^{1/2}.\end{eqnarray}$$

This notation is to be interpreted as follows: when $\unicode[STIX]{x1D6FD}>\unicode[STIX]{x1D6FD}_{c,mn}$ the wave $(m,n)$ is destabilized. It is evident from this formula that the critical $\unicode[STIX]{x1D6FD}$ numbers increase with radius $R$ , which is further visualized in figure 13. This figure also shows that we can define the critical $\unicode[STIX]{x1D6FD}$ number as $\unicode[STIX]{x1D6FD}_{c}=\min _{mn}(\unicode[STIX]{x1D6FD}_{c,mn})=\unicode[STIX]{x1D6FD}_{c,11}$ since the fundamental wave $(1,1)$ always is the first to be destabilized. With the present geometrical choices and material parameters, we calculate $\unicode[STIX]{x1D6FD}_{c}\in [0.33,1.06]$ for $R\in [0.5,5]~\text{m}$ . The values found here for $\unicode[STIX]{x1D6FD}_{c}$ are not so small and so it is best that we check if these high values of $\unicode[STIX]{x1D6FD}$ still remain accessible to our perturbation theory. Using (2.46) with $m=1$ , $\unicode[STIX]{x1D705}_{11}\approx 1.84$ , we find $\unicode[STIX]{x1D706}_{\mathit{v}}/\unicode[STIX]{x1D714}=0.0394\unicode[STIX]{x1D6FD}$ . Even with $\unicode[STIX]{x1D6FD}\sim 1$ , the necessary condition $\unicode[STIX]{x1D706}_{\mathit{v}}/\unicode[STIX]{x1D714}\ll 1$ remains fulfilled so we are confident about our prediction of  $\unicode[STIX]{x1D6FD}_{c}$ .

Figure 13.

In the presence of viscous damping, the critical $\unicode[STIX]{x1D6FD}$ -number for the onset of metal pad roll instability is non-zero. Using realistic material properties and fluid layer heights, we show here how critical $\unicode[STIX]{x1D6FD}$ numbers $\unicode[STIX]{x1D6FD}_{c,mn}$ for the destabilization of different waves $(m,n)$ , vary with radius $R$ .

The formula for $\unicode[STIX]{x1D6FD}_{c}$ and figure 13 give the impression that larger cells are more stable, but this is without taking into account the fact that the total current $I$ appearing in $\unicode[STIX]{x1D6FD}$ can still depend on $R$ . In electrolysis cells, just as in liquid metal batteries, it is likely that an upper bound exists for the current density $J$ which is of the order $10^{4}~\text{A}~\text{m}^{-2}$ . Replacing $I=J\unicode[STIX]{x03C0}R^{2}$ in the definition (1.1) of $\unicode[STIX]{x1D6FD}$ , we can then rewrite the condition $\unicode[STIX]{x1D6FD}>\unicode[STIX]{x1D6FD}_{c}$ for the onset of instability as a more practical condition on the radius of a cylindrical cell: $R>R_{c}$ . Using the material properties and the typical layer heights of our modelled reduction cell, $J=1.5\times 10^{4}~\text{A}~\text{m}^{-2}$ and $B_{z}=1~\text{mT}$ , we find $R_{c}=0.37~\text{m}$ . This is below the dimensions of actual reduction cells, but it is not so small either: viscosity would be rather efficient in keeping moderate-scale reduction cells stable.

5.2 Mg–Sb liquid metal battery

The references Weber et al. (Reference Weber, Beckstein, Galindo, Herreman, Nore, Stefani and Weier2017a ,Reference Weber, Beckstein, Herreman, Horstmann, Nore, Stefani and Weier b ) numerically investigate the metal pad roll instability in a $\text{Mg}|\text{NaCl}-\text{KCl}-\text{MgCl}_{2}|\text{Sb}$ liquid metal battery (LMB) of cylindrical shape. These batteries are structurally similar to reduction cells, except that they are composed of three layers of conducting fluid (top, middle, bottom $=$  indices 1, 2, 3) instead of two. The cell studied in the above references has dimensions

(5.4a,b ) $$\begin{eqnarray}R=0.05~\text{m},\quad (H_{1},H_{2},H_{3})=(0.045,0.01,0.045)~\text{m}.\end{eqnarray}$$

The material properties are

(5.5) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}(\unicode[STIX]{x1D70C}_{1},\unicode[STIX]{x1D70C}_{2},\unicode[STIX]{x1D70C}_{3})=(1577,1715,6270)~\text{kg}~\text{m}^{-3},\\ (\unicode[STIX]{x1D70E}_{1},\unicode[STIX]{x1D70E}_{2},\unicode[STIX]{x1D70E}_{3})=(3.6\times 10^{6},80,8.7\times 10^{5})~\text{S}~\text{m}^{-1},\\ (\unicode[STIX]{x1D708}_{1},\unicode[STIX]{x1D708}_{2},\unicode[STIX]{x1D708}_{3})=(6.7,6.8,2.0)\times 10^{-7}~\text{m}^{2}~\text{s}^{-1},\end{array}\right\}\end{eqnarray}$$

with the surface tensions

(5.6) $$\begin{eqnarray}(\unicode[STIX]{x1D6FE}_{1|2},\unicode[STIX]{x1D6FE}_{2|3})=(0.19,0.095)~\text{N}~\text{m}^{-1}.\end{eqnarray}$$

The top (Mg) and bottom (Sb) layers are very good conductors, but the middle layer, composed of molten salt $(\text{NaCl}-\text{KCl}-\text{MgCl}_{2})$ , is a badly conducting electrolyte. The numerical simulations performed in Weber et al. (Reference Weber, Beckstein, Galindo, Herreman, Nore, Stefani and Weier2017a ) are done with $B_{z}=10~\text{mT}$ and varying electrical currents $I=J\unicode[STIX]{x03C0}R^{2}$ . The metal pad roll instability is observed and growth rates are reported therein. Due to the large density ratio between the second and third layers, the instability mainly deforms the upper interface, separating the first and the second layers. This phenomenon is visible in figure 14(a), showing a snapshot of a saturated rotating wave. The lower interface only slightly deforms and the whole bottom layer does not really participate in the dynamics. Actually, closer inspection reveals that the frequency of the wave matches very closely that of the two-layer system composed of the first and the second layer. This observation has motivated us to try to apply our two-layer reduction cell theory to this three-layer LMB.

The following ‘technical’ operations are necessary in order to apply the theory. We entirely ignore region 3 and all that may occur there. At the interface of regions 2 and 3, the electrical boundary condition, as seen from the electrolyte is very similar to the one we had on the top cap of our cylindrical Hall–Héroult cell: the electrolyte meets there a much better conductor which we henceforth assume to be perfect. This means that if we exchange regions 1 and 2 and invert the sense of gravity i.e. $g\rightarrow -g$ , the theory and the formulas developed in the previous sections apply. The only difference is that the sense of rotation of the unstable wave has to be inverted. Finally, all the horizontal field effects are negligible here, just as in the small reduction cell model.

Figure 14.

(a) Snapshot of a numerical simulation of the metal pad roll instability in a Mg–Sb liquid metal battery. $J=10^{4}~\text{A}~\text{m}^{-2}$ , $B_{z}=10~\text{mT}$ . (b) Growth rate of the metal pad roll instability in the numerical simulations and according to our two-layer theory.

In figure 14, we compare the theoretical growth rates with simulations done with OpenFOAM for various values of the current. The dashed line corresponds to our theoretical estimate $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D706}_{\mathit{v}}+\unicode[STIX]{x1D706}_{\mathit{vv}}+\unicode[STIX]{x1D706}_{\mathit{visc}}$ . The data points are those from figure 9 in Weber et al. (Reference Weber, Beckstein, Galindo, Herreman, Nore, Stefani and Weier2017a ). The agreement is surprisingly good. The fact that we find the right slope with our two-layer theory implies that $\unicode[STIX]{x1D706}_{\mathit{v}}$ is correctly calculated. The threshold also seems approximatively right, which suggests that the combined magnetic and viscous damping rates $\unicode[STIX]{x1D706}_{\mathit{vv}}+\unicode[STIX]{x1D706}_{\mathit{visc}}$ are also more or less correct.

Since our model does quite a decent job in predicting the onset of instability, we want to use it to make some estimates for a much larger, more realistic battery system. We imagine a cell with shallow layers, keeping, for simplicity, the same fluid heights $H_{i}$ and material properties $\unicode[STIX]{x1D70C}_{i},\unicode[STIX]{x1D70E}_{i},\unicode[STIX]{x1D708}_{i}$ as in the cell of Weber et al. (Reference Weber, Beckstein, Galindo, Herreman, Nore, Stefani and Weier2017a ,Reference Weber, Beckstein, Herreman, Horstmann, Nore, Stefani and Weier b ) but greatly increasing its radius $R$ . We limit the current density to the realistic value of (i) $J=3~\text{k}\text{A}\text{m}^{-2}$ for Mg–Sb cells or to a more optimistic one (ii) $J=10~\text{k}\text{A}\text{m}^{-2}$ (both scenarios were also considered in Herreman et al. (Reference Herreman, Nore, Cappanera and Guermond2015)). We also suppose a weak external magnetic field $B_{z}\approx 1~\text{mT}$ that is more realistic. Using the formula (5.3) (invert $1\leftrightarrow 2$ and change $g\rightarrow -g$ ) we find the critical $\unicode[STIX]{x1D6FD}$ number

(5.7) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}_{c}\approx 2.81R^{1/2}\quad (R\text{ in units of m}),\end{eqnarray}$$

which yields

(5.8a,b ) $$\begin{eqnarray}\text{(i) }R_{c}=0.32~\text{m},\quad \text{(ii) }R_{c}=0.14~\text{m}\end{eqnarray}$$

as critical radii in the two current density scenarios. These values suggest that metal pad roll instability can indeed appear in liquid metal batteries with relatively moderate sizes, approximatively in the same range as the critical radii for the Tayler instability (Herreman et al. Reference Herreman, Nore, Cappanera and Guermond2015). A more detailed nonlinear study is required to be able to find the critical sizes at which the interface motion becomes too large and potentially is detrimental for the liquid metal battery’s integrity.