1. Introduction
Stably stratified liquid layers carrying strong vertical electric current are the core element of Hall–Héroult aluminium reduction cells (Evans & Ziegler Reference Evans and Ziegler2007) and the recently invented liquid metal batteries (Kim et al. Reference Kim2013). It is well known that an aluminium reduction cell may become unstable when the depth of electrolyte layer is decreased below a certain critical threshold which depends on the current strength and the cell design. The metal pad instability, as it is generally known (Urata Reference Urata1985), is caused by a rotating interfacial wave (Sele Reference Sele1977) and manifests itself as a sloshing of liquids in the cell (Davidson Reference Davidson2000). If the sloshing becomes too strong, it can disrupt the operation of the cell by creating a short circuit between the molten aluminium and the carbon anode, which is immersed in the electrolyte at the top. There are concerns that metal pad instability could affect also liquid metal batteries once they reach a sufficiently large size (Horstmann, Weber & Weier Reference Horstmann, Weber and Weier2018; Kelley & Weier Reference Kelley and Weier2018; Molokov Reference Molokov2018; Tucs, Bojarevics & Pericleous Reference Tucs, Bojarevics and Pericleous2018a,Reference Tucs, Bojarevics and Pericleousb). In contrast to the Tayler (internal pinch) instability, which is driven by the interaction of electric current with its own magnetic field and could also affect large-size liquid metal batteries (Stefani et al. Reference Stefani, Weier, Gundrum and Gerbeth2011; Herreman et al. Reference Herreman, Nore, Cappanera and Guermond2015; Weber et al. Reference Weber, Galindo, Priede, Stefani and Weier2015; Priede Reference Priede2017), the metal pad instability is caused by the presence of an external magnetic field, in particular, its vertical component, in the cell.
The mechanism behind the rotating interfacial wave is clear. Firstly, as the electric current flows through the poorly conducting electrolyte layer following the path of least resistance, a slight initial tilt of the interface redistributes current density from the depressed to the elevated parts of the interface. Secondly, in order to leave the system uniformly through the relatively poorly conducting carbon cathode, the current spreads out horizontally through the aluminium layer. As the horizontal electric current flows across the vertical magnetic field, a transverse electromagnetic force arises which tries to tilt the interface cross-wise to its initial slope so making it rotate.
This basic mechanism is captured by the experimental model of Pedchenko et al. (Reference Pedchenko, Molokov, Priede, Lukyanov and Thomas2009). Although the same mechanism is reproduced by the electromechanical pendulum model of Davidson & Lindsay (Reference Davidson and Lindsay1998), which exhibits a roll-type instability when the frequencies of two perpendicular oscillation modes occur sufficiency close, this model misses the fact that the metal pad instability is essentially a boundary effect (Lukyanov, El & Molokov Reference Lukyanov, El and Molokov2001). Namely, in the basic instability model, the electromagnetic force is mathematically reducible to a boundary condition which describes electromagnetically modified reflection of interfacial gravity waves by the sidewalls.
The simplest theoretical model of the metal pad instability is due to Lukyanov et al. (Reference Lukyanov, El and Molokov2001) who considered a semi-infinite domain (half-plane) bounded by a single lateral wall. They studied the effect of electromagnetic force on the reflection of interfacial gravity waves from the sidewalls, but missed the metal pad instability itself. This instability was first identified in the semi-infinite model by Morris & Davidson (Reference Morris and Davidson2003) and then revisited by Molokov, El & Lukyanov (Reference Molokov, El and Lukyanov2011). In this model, the interface can be destabilised by an arbitrary weak electromagnetic effect which gives rise to a nearly transverse wave travelling along the wall. In the infinite channel bounded by two parallel walls, which was considered first by Morris & Davidson (Reference Morris and Davidson2003), a sufficiently strong electromagnetic effect is required to destabilise the interface. In this model, the instability emerges as a longitudinal standing wave. These two basic models are reviewed and analysed in more detail in Appendix A.
The threshold of the metal pad instability in realistic finite-size cells is expected to depend on their geometry. For example, the square (Bojarevics & Romerio Reference Bojarevics and Romerio1994) and cylindrical (Davidson & Lindsay Reference Davidson and Lindsay1998; Lukyanov et al. Reference Lukyanov, El and Molokov2001) cells are known to be inherently unstable; like the semi-infinite cell, they can be destabilised by a relatively weak electromagnetic effect. The low stability threshold of the square cell is attributed to the occurrence of several pairs of gravity wave modes with the same frequency (Urata, Mori & Ikeuchi Reference Urata, Mori and Ikeuchi1976). The same applies also to cylindrical cells, which are considered as an option for the liquid metal batteries (Herreman et al. Reference Herreman, Nore, Guermond, Cappanera, Weber and Horstmann2019; Horstmann, Wylega & Weier Reference Horstmann, Wylega and Weier2019). According to the electromechanical model of Davidson & Lindsay (Reference Davidson and Lindsay1998), such cells are expected to be inherently unstable. However, their numerical results as well as those of Sneyd & Wang (Reference Sneyd and Wang1994) indicate that this is not in general the case as some rectangular cells with degenerate gravity wave spectra appear to have significantly higher stability thresholds than the square cell. Unfortunately, Sneyd & Wang (Reference Sneyd and Wang1994) did not notice this key result in their low-resolution numerical solution, which displays a very irregular variation of the instability threshold with the cell aspect ratio. Like Bojarevics & Romerio (Reference Bojarevics and Romerio1994) they claim all aspect ratios where two natural gravity wave frequencies coincide to be critical. Although the same result is predicted by the electromechanical pendulum model of Davidson & Lindsay (Reference Davidson and Lindsay1998), their numerical results obtained using the shallow-water model reveal that the gravity wave modes with the closest frequencies are not always the most unstable ones in rectangular cells. They explain this by the observation that only one in approximately five mode pairs are actually coupled and thus could become unstable.
Despite the relatively long history of this problem, it is not yet clear how the aspect ratio of rectangular cells affects their stability. This is a question of fundamental as well as practical importance which is addressed in the present study by carrying out a comprehensive theoretical and numerical analysis of the basic aluminium reduction cell model.
We show that the degeneracy is necessary but not in general sufficient for the instability because not all degenerate gravity wave modes are electromagnetically coupled. The critical aspect ratios predicted by the basic metal pad instability model are found to be limited to 
$\alpha =\sqrt {m/n},$ where 
$m$ and 
$n$ are any two odd numbers. These unstable aspect ratios form an absolutely discontinuous dense set of points which intersperse aspect ratios with finite (non-zero) stability thresholds. Although the fractality vanishes when viscous friction is taken into account, the instability threshold is smoothed out gradually and preserves its principal structure up to relatively high values of viscous friction.
The paper is organised as follows. In the next section, we formulate the problem and introduce a fully nonlinear two-layer shallow-water approximation which involves a novel integro-differential equation for the electric potential. Section 3 presents linear stability analysis of rectangular cells which is carried out analytically using an eigenvalue perturbation method as well as numerically using the Galerkin and the Chebyshev collocation methods. The main results are summarised and discussed in § 4. To make the paper self-contained, a brief review of analytical solutions for the semi-infinite and channel models is included in Appendix A.
2. Formulation of problem
2.1. Basic model and governing equations
Consider a shallow rectangular cell of size 
$L_{x}\times L_{y}$ bounded horizontally by insulating sidewalls and vertically by two solid electrodes which are separated by a constant distance 
$H\ll L$, as shown in figure 1. The electrodes carry electric current of a constant density 
$\boldsymbol {j}_{0}=-j_{0}\boldsymbol {e}_{z}$ which passes through a layer of liquid metal (aluminium) of constant depth 
$\bar {h}_{+}$ and density 
$\rho _{+}$ superposed by a layer of electrolyte of depth 
$\bar {h}_{-}=H-\bar {h}_{+}$ and density 
$\rho _{-}<\rho _{+}$. The electrical conductivity of the liquid metal 
$\sigma _{+}$ is supposed to be much higher than that of the electrodes 
$\sigma _{0}$ which, in turn, is much higher than that of the electrolyte: 
$\sigma _{-}\ll \sigma _{0}\ll \sigma _{+}.$ This approximation, which is used already by Sneyd & Wang (Reference Sneyd and Wang1994), is applicable to aluminium reduction cells because 
$\sigma _{-}/\sigma _{0}\sim \sigma _{0}/\sigma _{+}\sim 10^{-2}$ (Gerbeau, Le Bris & Lelièvre Reference Gerbeau, Le Bris and Lelièvre2006; Molokov et al. Reference Molokov, El and Lukyanov2011).
Sketch of the problem.
The two-layer system is subject to a downward gravity force with the free fall acceleration 
$g$ and external magnetic field 
$\boldsymbol {B}_{0}=B_{0}\boldsymbol {e}_{z}$ which is assumed to be vertical and uniform. Horizontal components of the magnetic field, which are known to have a stabilizing effect as long as they are predominantly generated by the current passing through the system (Sneyd Reference Sneyd1985; Herreman et al. Reference Herreman, Nore, Guermond, Cappanera, Weber and Horstmann2019), are as usual assumed to be negligible for the metal pad instability. This is because the effect of the horizontal magnetic field being determined by its horizontal gradient (Sneyd Reference Sneyd1985) is 
$\sim H/L=\varepsilon \ll 1$ relative to that of the vertical magnetic field which is determined by the magnitude of this field (Sneyd & Wang Reference Sneyd and Wang1994). The system admits a stably stratified quiescent equilibrium state with the layers separated by a plane horizontal interface. There is no electromagnetic force in this base state because the electric current is perfectly aligned with the magnetic field.
Note that a magnetostatic equilibrium in this system is in general possible if the magnetic field is vertically invariant. This, in turn, means that the vertical magnetic field, which is supposed to be generated by external currents, has to be also horizontally invariant (Sneyd & Wang Reference Sneyd and Wang1994). This is not the case in the model introduced by Urata et al. (Reference Urata, Mori and Ikeuchi1976) and also adopted by Bojarevics & Romerio (Reference Bojarevics and Romerio1994) who consider a horizontally inhomogeneous vertical component of the magnetic field but still assume a static base state.
Despite the highly idealised nature, this simple model captures the salient features of the Sele instability mechanism which is commonly believed to be relevant to aluminium reduction cells as well as to large-scale liquid metal batteries. For a comprehensive discussion of the relevance of the Sele instability mechanism (as well as a more detailed analysis of the assumptions underlying this model), we refer the reader to Gerbeau et al. (Reference Gerbeau, Le Bris and Lelièvre2006, § 6.1.4.3) who note the following: ‘Indeed, the model behind the Sele mechanism has been successfully used to improve the stability of some industrial cells… In addition, the Sele mechanism has been confirmed by complete 3D MHD [magnetohydrodynamic] simulations’.
Let us now consider a perturbed interface with the height 
$z=\zeta (\boldsymbol {x},t)$ which equals the depth of the bottom layer, 
$h_{+}(\boldsymbol {x},t)=\zeta (\boldsymbol {x},t)$, while the depth of the top layer is 
$h_{-}=H-h_{+}$. Following the shallow-water (long-wave) approximation, we assume that the perturbation has a characteristic longitudinal length scale 
$L\gg H$. Then the mass conservation implies that the vertical velocity 
$w$ is of relatively small magnitude 
$\sim H/L=\varepsilon \ll 1$. Respectively, the associated fluid flow is predominantly horizontal with the velocity 
$\boldsymbol {u}=u\boldsymbol {e}_{x}+v\boldsymbol {e}_{y}$. Such a flow has a negligible effect on the vertical pressure distribution which is thus purely hydrostatic:
\begin{equation} p_{{\pm}}(\boldsymbol{x},z,t)={\varPi}(\boldsymbol{x},t)+\rho_{{\pm}}g(\zeta-z), \end{equation}
where 
${\varPi }(\boldsymbol {x},t)=p_{\pm }(\boldsymbol {x},z,t)|_{z=\zeta }$ is the pressure distribution along the interface and 
$\boldsymbol {x}$ is the position vector in the 
$(x,y)$ plane. Substituting this pressure distribution into the horizontal fluid flow equation, we obtain
\begin{equation} \rho(\partial_{t}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u})={-}\boldsymbol{\nabla}({\varPi}+\rho g\zeta)+\boldsymbol{j}\times\boldsymbol{B}, \end{equation}
where 
$\boldsymbol {j}$ is the current density and B is the magnetic field; the 
$\pm$ indices are henceforth omitted for the sake of brevity.
The associated perturbation of the electric current density in each layer is governed by Ohm's law for a moving medium
\begin{equation}\boldsymbol{j}=\sigma(\boldsymbol{E}+\boldsymbol{u}\times\boldsymbol{B}), \end{equation}
where 
$\boldsymbol {E}$ is the perturbation of the electric field in the laboratory frame of reference. The oscillation period of interfacial waves, 
$\tau _{0}\sim L/c$, where 
$c$ is the wave speed (2.31), is assumed to be much longer than the characteristic magnetic diffusion time 
$\tau _{m}\sim \mu _{0}\sigma H^{2}$, where 
$\mu _{0}=4{\rm \pi} \times 10^{-7}\ \textrm{H} \ {\rm m}^{-1}$ is the vacuum permeability. This is the case in aluminium reduction cells because 
$\tau _{m}/\tau _{0}\sim \mu _{0}\sigma cH^{2}/L\sim \varepsilon Rm\ll 1,$ where not only 
$\varepsilon =H/L$ but also the magnetic Reynolds number 
$Rm=\mu _{0}\sigma cH$ is usually small. It means that the induction effect due to the temporal variation of the magnetic field perturbation is negligible. Hence, the Faraday law of induction reduces to 
${\boldsymbol {\nabla }\times \boldsymbol {E}=0}$, which corresponds to the so-called quasi-stationary or inductionless approximation commonly used in the liquid-metal MHD (Roberts Reference Roberts1967). Then we have 
$\boldsymbol {E}=-\boldsymbol {\nabla }\varphi$, where 
$\varphi$ is the electric potential.
In addition, the perturbation of the electric field induced by the flow across the magnetic field, 
$\boldsymbol {u}\times \boldsymbol {B}$, is assumed to negligible relative to the horizontal electric field perturbation caused by the interface deformation. Using the charge conservation 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {j}=0$, the latter can be estimated as 
$\sim j_{0}L/(\sigma H),$ whereas for the former we have 
$\sim B_{0}c,$ which thus leads to 
$HB_{0}c\sigma /(Lj_{0})\ll 1.$ For 
$B_{0}\sim \mu _{0}j_{0}L,$ this reduces to 
$Rm\ll 1,$ which, as noted above, is typically the case, and confirms the applicability of the aforementioned inductionless approximation. Note that in this approximation, also the magnetic field induced by the fluid flow is negligible relative to the external magnetic field. On the other hand, the perturbation of the vertical magnetic field component produced by the interface deformation with amplitude 
$\sim H$ can be evaluated as 
$\sim B_{0}H/L.$ Thus, we have 
$\boldsymbol {B}=B_{0}\boldsymbol {e}_{z}+O(\varepsilon ,Rm).$
2.2. Reduction of the electric potential equation
Taking into account the above estimates, Ohm's law (2.3) reduces to 
$\boldsymbol {j}=-\sigma \boldsymbol {\nabla }\varphi$ and, thus, the charge conservation 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {j}=0$ results in
\begin{equation} \nabla^{2}\varphi=0, \end{equation}
which governs the electric potential in both layers. At the interface 
$z=\zeta (\boldsymbol {x},t)$, we have the continuity of both the potential and the normal component of electric current:
\begin{equation} [\varphi]=[j_{n}]=0, \end{equation}
where 
$[f]$ stands for a jump of 
$f$ across the interface. The top electrode (anode) is assumed to be perfectly conducting relative to the electrolyte, which means that the interface between both is effectively equipotential. Since the electric potential is defined up an additive constant, we can set
\begin{equation} \left.\varphi_{-}\right|_{z=H}=0. \end{equation}As the electrical conductivity of the bottom electrode (cathode) is much lower than that of aluminium, the former may be regarded as effectively insulating with respect to internal current perturbations in the two-layer system. It means that the current density at the cathode is fixed
\begin{equation} \left.\sigma_{+}\partial_{z}\varphi_{+}\right|_{z=0}=j_{0}. \end{equation}
Since 
$H/L=\varepsilon \ll 1$, the potential distribution in the aluminium layer can be approximated by the power series expansion in 
$z$ (Bojarevics & Romerio Reference Bojarevics and Romerio1994):
\begin{equation} \varphi_{+}(\boldsymbol{r})=\phi_{+}^{(0)}(\boldsymbol{x})+z\phi_{+}^{(1)}(\boldsymbol{x})+\frac{z^{2}}{2}\phi_{+}^{(2)}(\boldsymbol{x})+O(\varepsilon^{3}), \end{equation}
where 
$\boldsymbol {x}=\boldsymbol {r}-z\boldsymbol {e}_{z}=x\boldsymbol {e}_{x}+y\boldsymbol {e}_{y}$ and 
$\phi _{+}^{(k)}(\boldsymbol {x})\equiv \partial _{z}^{k}\varphi _{+}|_{z=0}$. Firstly, requiring (2.8) to satisfy (2.7) and (2.4), we find
\begin{gather} \phi_{+}^{(1)} = j_{0}/\sigma_{+}, \end{gather}
\begin{gather}\phi_{+}^{(2)} ={-}\nabla^{2}\phi_{+}^{(0)}+O(\varepsilon). \end{gather}Secondly, using the interface normal
\begin{equation} \boldsymbol{n}=\frac{\boldsymbol{e}_{z}-\boldsymbol{\nabla}h_{+}}{|\boldsymbol{e}_{z}-\boldsymbol{\nabla}h_{+}|}=\boldsymbol{e}_{z}-\boldsymbol{\nabla}h_{+}+\boldsymbol{e}_{z}O(\varepsilon^{2}), \end{equation}the normal current at the interface can be expressed as
\begin{equation} \left.j_{n}\right|_{z=h_{+}}={-}\sigma_{+}\left.\partial_{n}\varphi_{+}\right|_{z=h_{+}} \approx{-}\sigma_{+}(\phi_{+}^{(1)}+h_{+}\phi_{+}^{(2)}-\boldsymbol{\nabla}h_{+}\boldsymbol{\cdot} \boldsymbol{\nabla}\phi_{+}^{(0)})+O(\varepsilon^{2}).\end{equation}Finally, substituting (2.9) and (2.10) into (2.12), after a few rearrangements, we obtain
\begin{equation} \sigma_{+}\boldsymbol{\nabla}\boldsymbol{\cdot}(h_{+}\boldsymbol{\nabla}\phi_{+}^{(0)})=j_{0}+\left.j_{n}\right|_{z=h_{+}}+O(\varepsilon^{2}), \end{equation}which is a reduced two-dimensional potential equation for the aluminium layer.
In order to determine the normal current 
$j_{n}|_{z=h_{+}}$ on the right-hand side of (2.13), we need to consider the potential distribution in the electrolyte layer, which can be sought as in the aluminium layer:
\begin{equation} \varphi_{-}(\boldsymbol{r})=\phi_{-}^{(0)}(\boldsymbol{x})+(z-H)\phi_{-}^{(1)}(\boldsymbol{x})+\tfrac{1}{2}(z-H)^{2}\phi_{-}^{(2)}(\boldsymbol{x})+O(\varepsilon^{3}), \end{equation}
where 
$\phi _{-}^{(k)}(\boldsymbol {x})\equiv \partial _{z}^{k}\varphi _{-}|_{z=H}$. Then applying (2.6) and (2.4), we find
\begin{equation} \phi_{-}^{(0)}=\phi_{-}^{(2)}=0. \end{equation}Similar to the top electrode, the aluminium layer is a much better conductor than the electrolyte layer. It means that the aluminium is effectively equipotential with respect to the electrolyte:
\begin{equation} \left.\varphi_{-}\right|_{z=h_{+}}=(h_{+}-H)\phi_{-}^{(1)}(\boldsymbol{x})+O(\varepsilon^{3})=\varphi_{0}, \end{equation}
where 
$\varphi _{0}$ is an unknown constant potential of the aluminium layer. As a result, we have
\begin{equation} j_{z}={-}\sigma_{-}\left.\partial_{z}\varphi_{-}\right|_{z=H}=\sigma_{-}\varphi_{0}/h_{-}+O(\varepsilon^{2}), \end{equation}
where 
$h_{-}=H-h_{+}$ is the depth of the electrolyte layer, and thus
\begin{equation} \left.j_{n}\right|_{z=h_{+}}=j_{z}/n_{z}=j_{z}+O(\varepsilon^{2}) =\sigma_{-}\varphi_{0}/h_{-}+O(\varepsilon^{2}).\end{equation}Finally, substituting (2.18) into (2.13), we obtain
\begin{equation} \sigma_{+}\boldsymbol{\nabla}\boldsymbol{\cdot}(h_{+}\boldsymbol{\nabla}\phi_{+}^{(0)}) =j_{0}+\sigma_{-}\varphi_{0}h_{-}^{{-}1},\end{equation}
which is accurate up to 
$O(\varepsilon ^{2})$.
The unknown potential of the aluminium layer 
$\varphi _{0}$ is determined by the solvability condition of (2.19) which is required by the Neumann boundary condition at the insulating sidewalls where we have
\begin{equation} \left.\partial_{n}\varphi_{+}\right|_{\varGamma} =\left.\partial_{n}\phi_{+}^{(0)}\right|_{\varGamma}+O(\varepsilon^{2})=0.\end{equation}
Integrating (2.19) over the horizontal cross-section area 
$S$ and using this boundary condition, we obtain
\begin{equation} \int_{S}(\,j_{0}+\sigma_{-}\varphi_{0}h_{-}^{{-}1})\,\mathrm{d}^{2}\boldsymbol{x}=0, \end{equation}
which is the solvability condition of (2.19). This condition, which requires the constancy of the total current 
$j_{0}S=I_{0}$, defines the potential of the aluminium layer:
\begin{equation} \varphi_{0}={-}\sigma_{-}^{{-}1}I_{0}\left/\int_{S}h_{-}^{{-}1}\,\mathrm{d}^{2}\boldsymbol{x}\right. . \end{equation}Substituting this into (2.19), we obtain
\begin{equation} \sigma_{+}\boldsymbol{\nabla}\boldsymbol{\cdot}(h_{+}\boldsymbol{\nabla}\phi_{+}^{(0)}) =j_{0}\left(1-Sh_{-}^{{-}1}\left/\int_{S}h_{-}^{{-}1}\,\mathrm{d}^{2}\boldsymbol{x}\right.\right), \end{equation}which is the final form of the reduced electric potential equation for the aluminium layer with the boundary condition (2.20).
Note that (2.23) is fully nonlinear and, thus, applicable not only to small but also to large amplitude waves. Since the present linear stability analysis is concerned with the former, only the linearised version of (2.23) will be required in the following. Although the linear form of the potential equation has been derived and used in many previous studies, the full nonlinear equation (2.23) is new to the best of our knowledge. For example, Bojarevics (Reference Bojarevics1998) considers nonlinear waves in aluminium electrolysis cells but uses linear potential equation (17) derived earlier by Bojarevics & Romerio (Reference Bojarevics and Romerio1994). Equation (2.23) differs also from analogous equation (11) of Zikanov et al. (Reference Zikanov, Thess, Davidson and Ziegler2000). Instead of expanding the potential in powers of 
$z$, which is essential in this case (Bojarevics & Romerio Reference Bojarevics and Romerio1994), they integrate a three-dimensional potential equation over the depth of the aluminium layer. As this operation cannot consistently be carried through, Zikanov et al. (Reference Zikanov, Thess, Davidson and Ziegler2000) effectively postulate that the integral of the horizontal current over the depth of the aluminium layer can be written as 
$\boldsymbol {J}=-\sigma _{+}\boldsymbol {\nabla }\varPsi ,$ whereas we have 
$\boldsymbol {J}=-\sigma _{+}h_{+}\boldsymbol {\nabla }\phi _{+}^{(0)},$ which follows from Ohm's law and the 
$z$-invariance of 
$\phi _{+}^{(0)}.$ Both expressions are equivalent only for small-amplitude waves when 
$h_{+}\approx \bar {h}_{+}=\text {const.}$ and, thus, 
$\varPsi =\bar {h}_{+}\phi _{+}^{(0)}.$ Consequently, equation (11) of Zikanov et al. (Reference Zikanov, Thess, Davidson and Ziegler2000) is limited to small-amplitude waves.
2.3. Linearised shallow-water model
The distribution of the electric potential in the aluminium layer (2.8) indicates that the respective electromagnetic force
\begin{equation} \boldsymbol{j_{+}}\times\boldsymbol{B_{0}}=\sigma_{+}\boldsymbol{B_{0}}\times\boldsymbol{\nabla}\phi_{+}^{(0)}(\boldsymbol{x})+O(\varepsilon^{2}) \end{equation}
is depth invariant up to 
$O(\varepsilon ^{2})$. Since the curl of this force has zero horizontal components, it preserves zero horizontal vorticity of the flow. In the hydrostatic shallow-water approximation, when the vertical velocity component is negligible, this is equivalent to the depth invariance of the horizontal velocity: 
$\partial _{z}\boldsymbol {u}\equiv 0$. Consequently, this electromagnetic force is compatible with the conservation of the depth invariance which is a generic property of the two-layer shallow-water system (Priede Reference Priede2020). In the electrolyte layer, according to (2.14) and (2.15), the horizontal component of the current perturbation and, thus, also the associated electromagnetic force is 
$O(\varepsilon )$, i.e. negligible in the hydrostatic shallow-water approximation. Therefore, the conservation of the depth invariance holds also in the electrolyte layer.
The system of shallow-water equations is closed by adding the conservation of mass in each layer:
\begin{equation} \partial_{t}h+\boldsymbol{\nabla}\boldsymbol{\cdot}(h\boldsymbol{u})=0, \end{equation}
which is obtained by integrating the incompressibility constraint 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}=0$ over the depth of the respective layer and using the kinematic constraint at the interface 
$z=\zeta (\boldsymbol {x},t)$:
\begin{equation} \frac{\textrm{d}\zeta}{\textrm{d}t}=\partial_{t}\zeta+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\zeta=w. \end{equation}In the following, we focus on the interfacial waves of small amplitude:
\begin{equation} \eta(\boldsymbol{x},t)=\zeta(\boldsymbol{x},t)-\bar{h}_{+}\ll H. \end{equation}It means that the nonlinear terms in the governing equations are higher-order small and, thus, negligible relative to the linear terms. After linearisation, (2.2) and (2.25) take the form
\begin{gather} \rho_{{\pm}}\partial_{t}\boldsymbol{u}_{{\pm}} ={-}\boldsymbol{\nabla}({\varPi}+\rho_{{\pm}}g\eta)+\sigma_{{\pm}}\boldsymbol{B}_{0}\times\boldsymbol{\nabla}\phi_{{\pm}}^{(0)}, \end{gather}
\begin{gather}\partial_{t}\eta ={\mp}\bar{h}_{{\pm}}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_{{\pm}}. \end{gather}Taking the divergence of the difference of (2.28) for the top and bottom layers, which correspond to the plus and minus indices, and using (2.29), we obtain the classical wave equation
\begin{equation} \partial_{t}^{2}\eta=c^{2}\nabla^{2}\eta, \end{equation}for interfacial gravity waves propagating with the speed
\begin{equation} c=\left(\frac{g(\rho_{+}-\rho_{-})}{\rho_{+}/\bar{h}_{+} +\rho_{-}/\bar{h}_{-}}\right)^{1/2}.\end{equation} Note that since the electromagnetic force is entirely solenoidal, i.e. 
$\boldsymbol {B}_{0}\times \boldsymbol {\nabla }\phi _{\pm }^{(0)}=-\boldsymbol {\nabla }\times \phi _{\pm }^{(0)}\boldsymbol {B}_{0}$, it can drive a flow without disturbing pressure as long as the flow is not obstructed by the sidewalls. It implies that the electromagnetic force can affect interfacial waves only at the sidewalls. The boundary condition for (2.30) follows from 
$u_{n}|_{{\varGamma }}=0$ which holds at the impermeable sidewall 
${\varGamma }$. Taking again the difference of (2.28) for the top and bottom layers to eliminate the pressure and then applying the impermeability condition, we obtain
\begin{equation} \left.g[\rho]\partial_{n}\eta\right|_{{\varGamma}} =\left.B_{0}\partial_{\tau}[\sigma\phi^{(0)}]\right|_{{\varGamma}}, \end{equation}
where 
$\partial _{n}\equiv \boldsymbol {n}\boldsymbol {\cdot }\boldsymbol {\nabla }$ and 
$\partial _{\tau }\equiv \boldsymbol {\tau }\boldsymbol {\cdot }\boldsymbol {\nabla }$ with 
$\boldsymbol {n}$ and 
$\boldsymbol {\tau }=\boldsymbol {n}\times \boldsymbol {e}_{z}$ standing for the outward normal and tangent vectors to the lateral cell boundary 
${\varGamma };$ the square brackets denote the difference of the enclosed quantity between the top and bottom layers.
After linearisation, that is, applying
\begin{gather} h_{+} = \bar{h}_{+}+\eta\approx\bar{h}_{+}, \end{gather}
\begin{gather}h_{-}^{{-}1} = (\bar{h}-\eta)^{{-}1}\approx\bar{h}_{-}^{{-}1}+\eta\bar{h}_{-}^{{-}2}, \end{gather}
and taking into account that 
$\int _{s}\eta \,\mathrm {d}^{2}\boldsymbol {x}=0$ due to the mass conservation, the potential equation (2.23) reduces to
\begin{equation} \nabla^{2}\phi_{+}^{(0)}={-}\frac{j_{0}\eta}{\sigma_{+}\bar{h}_{+}\bar{h}_{-}}, \end{equation}
which satisfies the solvability condition automatically owing to the mass conservation. This equation originally formulated by Urata et al. (Reference Urata, Mori and Ikeuchi1976) differs from the analogous equation derived later by Bojarevics & Romerio (Reference Bojarevics and Romerio1994) by the absence of terms containing the ratio of electric conductivities of electrolyte and aluminium, 
$\sim \sigma _{-}/\sigma _{+}\sim 10^{-4},$ which is 
$\sim (H/L)^{2},$ i.e. comparable to the higher-order-small terms usually neglected in the hydrostatic shallow-water approximation (Lukyanov et al. Reference Lukyanov, El and Molokov2001).
Henceforth, we change to dimensionless variables by using 
$L=\sqrt {L_{x}L_{y}}=S^{1/2}$, 
$\tau _{0}=L/c$ and 
$\phi _{0}=I_{0}/(\sigma _{+}\bar {h}_{+})$ as the length, time and electric potential scales, respectively. The dimensionless governing equations (2.30) and (2.35) read as
\begin{gather} \partial_{t}^{2}\eta+\gamma\partial_{t}\eta = \nabla^{2}\eta, \end{gather}
\begin{gather}-\eta = \nabla^{2}\phi, \end{gather}
where 
$\phi$ is the dimensionless counterpart of 
$\phi _{+}^{(0)}$. Note that following Lukyanov et al. (Reference Lukyanov, El and Molokov2001), we have added a linear damping (viscous friction) term to (2.36) defined by a phenomenological friction parameter 
$\gamma$. The latter is subsequently assumed to be small and, thus, ignored, unless stated otherwise. The dimensionless boundary conditions corresponding to (2.20) and (2.32) are
\begin{gather} \left.\partial_{n}\eta\right|_{{\varGamma}} = \beta\left.\partial_{\tau}\phi\right|_{{\varGamma}}, \end{gather}
\begin{gather}\left.\partial_{n}\phi\right|_{{\varGamma}} = 0, \end{gather}where
\begin{equation} \beta=\frac{I_{0}B_{0}}{g(\rho_{+}-\rho_{-})\bar{h}_{+}\bar{h}_{-}} \end{equation}is the key dimensionless parameter for this problem, referred to by Gerbeau et al. (Reference Gerbeau, Le Bris and Lelièvre2006) as the Sele number (Sele Reference Sele1977), which defines the strength of electromagnetic force relative to that of gravity.
3. Linear stability analysis
3.1. Eigenvalue perturbation solution for 
$\beta \ll 1$
Let us now turn to the linear stability problem posed by (2.36) and (2.37) for a rectangular cell of aspect ratio 
$\alpha =L_{x}/L_{y}$ which is bounded laterally by four sidewalls. Although the problem does not admit an exact analytical solution as for the semi-infinite and channel geometries (see Appendix A), it can still be solved approximately for sufficiently small 
$\beta$ using the classical eigenvalue perturbation method (Hinch Reference Hinch1991, § 1.6). To this end, the eigenmode
\begin{equation} \{\eta,\phi\}(\boldsymbol{x},t)=\{\hat{\eta},\hat{\phi}\}(\boldsymbol{x})\textrm{e}^{-\textrm{i}\omega t}+\text{c.c.}, \end{equation}
where c.c. stands for the complex conjugate, is sought by expanding the eigenvalue 
$\lambda =\omega ^{2}$ and the respective amplitude distribution in the power series of 
$\beta :$
\begin{gather} \omega^{2} = \lambda^{(0)}+\beta\lambda^{(1)}+\ldots, \end{gather}
\begin{gather}\{\hat{\eta},\hat{\phi}\}(\boldsymbol{x}) = \{\hat{\eta},\hat{\phi}\}^{(0)}(\boldsymbol{x})+\beta\{\hat{\eta},\hat{\phi}\}^{(1)}(\boldsymbol{x})+\ldots. \end{gather}
At the leading order, which corresponds to 
$\beta =0$, (2.36)–(2.39) reduce to
\begin{gather} \lambda^{(0)}\hat{\eta}^{(0)}+\nabla^{2}\hat{\eta}^{(0)}=0,\quad \left.\partial_{n}\hat{\eta}^{(0)}\right|_{{\varGamma}}=0, \end{gather}
\begin{gather}\hat{\eta}^{(0)}+\nabla^{2}\hat{\phi}^{(0)}=0, \quad \left.\partial_{n}\hat{\phi}^{(0)}\right|_{{\varGamma}}=0. \end{gather}The respective leading-order solution is
\begin{gather} \lambda_{\boldsymbol{k}}^{(0)} = \boldsymbol{k}^{2}, \end{gather}
\begin{gather}\{\hat{\eta},\hat{\phi}\}^{(0)}(\boldsymbol{x}) = \{\boldsymbol{k}^{2},1\}\hat{\phi}{}_{\boldsymbol{k}}^{(0)}\varPsi_{\boldsymbol{k}}(\boldsymbol{x}), \end{gather}where
\begin{equation} \varPsi_{\boldsymbol{k}}(\boldsymbol{x})=\cos(xk_{x})\cos(yk_{y}) \end{equation}is the gravity wave mode with the wave vector
\begin{equation} \boldsymbol{k}=(k_{x},k_{y})={\rm \pi}(m/\sqrt{\alpha},n\sqrt{\alpha}),\quad m,n=0,1,2,\ldots. \end{equation} It is important to note that some leading-order eigenmodes (3.7) may be degenerate, namely, to have the same eigenfrequency (3.6) for different wave vectors. In this case, the leading-order solution can be a superposition of the eigenmodes (3.7) with the same magnitude of the wave vector 
$\boldsymbol {k}'^{2}=\boldsymbol {k}^{2}:$
\begin{equation} \{\hat{\eta},\hat{\phi}\}^{(0)}(\boldsymbol{x}) =\sum_{\boldsymbol{\boldsymbol{k}'^{2}=\boldsymbol{k}^{2}}}\{\hat{\eta},\hat{\phi}\}_{\boldsymbol{k'}}^{(0)}\varPsi_{\boldsymbol{k'}}(\boldsymbol{x}). \end{equation} Following the standard approach, the first-order correction 
$\{\hat {\eta },\hat {\phi }\}^{(1)}$ to (3.10) is sought as an expansion in the leading-order eigenmodes
\begin{equation} \{\hat{\eta},\hat{\phi}\}^{(1)}(\boldsymbol{x})= \sum_{\boldsymbol{k}}\{\hat{\eta},\hat{\phi}\}_{\boldsymbol{k}}^{(1)}\varPsi_{\boldsymbol{k}}(\boldsymbol{x}), \end{equation}where the summation is over all wave vectors. Substituting this expansion into the first-order problem
\begin{gather} \lambda^{(0)}\hat{\eta}^{(1)}+\nabla^{2}\hat{\eta}^{(1)}={-}\lambda^{(1)}\hat{\eta}^{(0)}, \quad \left.\partial_{n}\hat{\eta}^{(1)}\right|_{{\varGamma}}=\left. \partial_{\tau}\hat{\phi}^{(0)}\right|_{\varGamma}, \end{gather}
\begin{gather}\hat{\eta}^{(1)}+\nabla^{2}\hat{\phi}^{(1)}=0, \quad \left.\partial_{n}\hat{\phi}^{(1)}\right|_{{\varGamma}}=0, \end{gather}
and using the orthogonality of leading-order eigenmodes by projecting the result onto 
$\varPsi _{\boldsymbol {k}}$, after a few rearrangements, we obtain the solvability condition of the first-order problem:
\begin{equation}\lambda_{\boldsymbol{k}}^{(1)}\boldsymbol{k}^{2}\langle \varPsi_{\boldsymbol{k}}^{2}\rangle \hat{\phi}_{\boldsymbol{k}}^{(0)}-\sum_{\boldsymbol{k}'^{2}=\boldsymbol{k}^{2}} F_{\boldsymbol{k},\boldsymbol{k}'}\hat{\phi}_{\boldsymbol{k}'}^{(0)} =(\boldsymbol{k}^{2}-\lambda_{\boldsymbol{k}}^{(0)})\boldsymbol{k}^{2}\langle \varPsi_{\boldsymbol{k}}^{2}\rangle \hat{\phi}_{\boldsymbol{k}}^{(1)}=0, \end{equation}
where the right-hand side is zero owing to the leading-order solution (3.6) and, thus, the left-hand side defines the eigenvalue perturbation 
$\lambda _{\boldsymbol {k}}^{(1)}.$ The summation on the left-hand side involving the electromagnetic interaction matrix 
$F_{\boldsymbol{k},\boldsymbol{k}'}$ is over the degenerate modes as in (3.10); for non-degenerate modes, it reduces to a single term with 
$\boldsymbol {k}'=\boldsymbol {k}$. The angle brackets denote integral over 
$S=L_{x}\times L_{y}:$
\begin{equation} \langle \varPsi_{\boldsymbol{k}}^{2}\rangle = \int_{S}\varPsi_{\boldsymbol{k}}^{2}\,\mathrm{d}^{2}\boldsymbol{x}=c_{k_{x}}^{{-}1}c_{k_{y}}^{{-}1}, \end{equation}
where 
$c_{0}=1$ and 
$c_{k}=2$ for 
$k\not =0$. The second left-hand side term of (3.14) is due to Green's first identity
\begin{equation} \langle \varPsi_{\boldsymbol{k}}\nabla^{2}\hat{\eta}^{(1)}+\boldsymbol{\nabla}\varPsi_{\boldsymbol{k}}\boldsymbol{\cdot}\boldsymbol{\nabla}\hat{\eta}^{(1)}\rangle ={\oint\limits_\varGamma}\varPsi_{\boldsymbol{k}}\partial_{n}\hat{\eta}^{(1)}\,\mathrm{d}{\varGamma}, \end{equation}where the boundary integral can be transformed using (3.12b) and Green's theorem as follows:
\begin{equation}{\oint\limits_\varGamma}\varPsi_{\boldsymbol{k}}\partial_{\tau}\hat{\phi}^{(0)}\,\mathrm{d}{\varGamma}={-}\langle \boldsymbol{e}_{z}\boldsymbol{\cdot}\boldsymbol{\nabla}\varPsi_{\boldsymbol{k}}\times\boldsymbol{\nabla}\hat{\phi}^{(0)}\rangle . \end{equation}This results in the electromagnetic interaction matrix
\begin{equation} F_{\boldsymbol{k},\boldsymbol{k}'}=\langle \boldsymbol{e}_{z}\boldsymbol{\cdot}\boldsymbol{\nabla}\varPsi_{\boldsymbol{k}}\times\boldsymbol{\nabla}\varPsi_{\boldsymbol{k}'}\rangle =G_{\boldsymbol{k},\boldsymbol{k}'}-G_{\boldsymbol{k}',\boldsymbol{k}}, \end{equation}
where 
$G_{\boldsymbol {k},\boldsymbol {k}'}=H_{k_{x},k_{x}'}H_{k_{y}',k_{y}}$ and
\begin{equation} H_{k_{m},k_{n}}=k_{m}\left\langle \sin(xk_{m})\cos(xk_{n})\right\rangle =\frac{2m^{2}}{m^{2}-n^{2}}\text{mod}(m+n,2). \end{equation}Expressions (3.18) and (3.19) differ from those used by Sneyd & Wang (Reference Sneyd and Wang1994) only by a more compact form.
Since the electromagnetic interaction matrix (3.18) is anti-symmetric, as it is well known, there is no electromagnetic back-reaction on separate gravity wave modes, i.e. 
$F_{\boldsymbol {k},\boldsymbol {k}}=0$. Thus, for a non-degenerate (single) mode, we have
\begin{equation} \lambda_{\boldsymbol{k}}^{(1)}\langle \varPsi_{\boldsymbol{k}}^{2}\rangle k^{2}\hat{\phi}_{\boldsymbol{k}}=F_{\boldsymbol{k},\boldsymbol{k}}\hat{\phi}_{\boldsymbol{k}}=0, \end{equation}
which means no electromagnetic effect at the first order in 
$\beta$. For a degenerate mode consisting of a superposition of two eigenmodes with the same frequency,
\begin{equation} \lambda_{\boldsymbol{k}}^{(0)}=\boldsymbol{k}_{1}^{2}=\boldsymbol{k}_{2}^{2}, \end{equation}(3.14) can be written in the matrix form as
\begin{equation} \lambda_{\boldsymbol{k}}^{(1)}k^{2} \begin{pmatrix} \langle \varPsi_{\boldsymbol{k}_{1}}^{2}\rangle \hat{\phi}_{\boldsymbol{k}_{1}}\\ \langle \varPsi_{\boldsymbol{k}_{2}}^{2}\rangle \hat{\phi}_{\boldsymbol{k}_{2}} \end{pmatrix}=\begin{pmatrix} 0 & F_{\boldsymbol{k}_{1},\boldsymbol{k}_{2}}\\ -F_{\boldsymbol{k}_{1},\boldsymbol{k}_{2}} & 0 \end{pmatrix}\begin{pmatrix} \hat{\phi}_{\boldsymbol{k}_{1}}\\ \hat{\phi}_{\boldsymbol{k}_{2}} \end{pmatrix}.\end{equation}The solution of this second-order matrix eigenvalue problem yields
\begin{equation} \lambda_{\boldsymbol{k}}^{(1)}k^{2}\sqrt{\langle \varPsi_{\boldsymbol{k}_{1}}^{2}\rangle \langle \varPsi_{\boldsymbol{k}_{2}}^{2}\rangle}={\pm} \textrm{i}F_{\boldsymbol{k}_{1},\boldsymbol{k}_{2}}, \end{equation}
which indicates that 
$\lambda _{\boldsymbol {k}}^{(1)}$ is imaginary and, thus, the system is unstable provided that 
$F_{\boldsymbol {k}_{1},\boldsymbol {k}_{2}}\not =0$. According to (3.19), this is the case only if both components of the two wave vectors 
$\boldsymbol {k}_{1}={\rm \pi} (m_{1}/\sqrt {\alpha },n_{1}\sqrt {\alpha })$ and 
$\boldsymbol {k}_{2}={\rm \pi} (m_{2}/\sqrt {\alpha },n_{2}\sqrt {\alpha })$ have opposite parities; namely, 
$m_{1}\pm m_{2}$ and 
$n_{1}\pm n_{2}$ are odd numbers. Then the degeneracy condition (3.21) yields
\begin{equation} \alpha_{c}^{2}={-}\frac{(m_{1}+m_{2})(m_{1}-m_{2})}{(n_{1}+n_{2})(n_{1}-n_{2})}=\frac{m}{n}, \end{equation}
where 
$m$ and 
$n$ can be any odd numbers. Consequently, all cells with aspect ratios squared equal to the ratio of two odd numbers are inherently unstable, i.e. they become unstable at infinitesimal 
$\beta >\beta _{c}=0$.
For 
$\alpha ^{2}=1$, which corresponds to a square cell, the unstable wavenumbers are 
$m_{1}=l, n_{1}=0$ and 
$m_{2}=0, n_{2}=l$, where 
$l=1,3,5,\ldots$. In this case, (3.23) yields 
$\lambda _{l}^{(1)}=\pm \textrm {i}8/(l{\rm \pi} )^{2}$ and hence the instability growth rate resulting from (3.2) is
\begin{equation} \mathrm{Im}[\omega_{l}^{(1)}]={\pm}\frac{4}{(l{\rm \pi})^{3}},\end{equation}
which means that the most unstable is the mode with 
$l=1$. For a general 
$\alpha ^{2}=m/n$ defined by arbitrary odd numbers 
$m$ and 
$n$, the lowest unstable wavenumbers can be written as
\begin{gather} m_{(3\pm1)/2} = (m\pm1)/2, \end{gather}
\begin{gather}n_{(3\pm1)/2} = (n\mp1)/2. \end{gather}
In this case, (3.23) yields 
$\lambda _{m,n}^{(1)}=\pm \textrm {i}8c_{m-1}^{1/2}c_{n-1}^{1/2}/({\rm \pi} ^{2}\sqrt {mn})$ and, respectively, (3.2) results in
\begin{equation} \mathrm{Im}[\omega_{m,n}^{(1)}]={\pm}\frac{8}{{\rm \pi}^{3}}\left(\frac{c_{m-1}c_{n-1}}{(m+n)(1+mn)\sqrt{mn}}\right)^{1/2}, \end{equation}
which reduces to (3.25) with 
$l=1$ when 
$m=n=1$. It is important to note that the growth rate (3.28) drops off with the increase of 
$m$ and 
$n$ as 
$\sim (m+n)^{-1/2}(mn)^{-3/4}.$
For an aspect ratio 
$\alpha$ sufficiently close to the critical value 
$\alpha _{c}=\sqrt {m/n},$ the stability of the system is expected to be determined by the interaction of two critical modes with the wavenumbers ((3.26) and (3.27)) which correspond to the wavenumbers
\begin{equation} \boldsymbol{k}_{(3\pm1)/2}=\frac{\rm \pi}{2}((m\pm1)/\sqrt{\alpha},(n\mp1)\sqrt{\alpha}). \end{equation}Searching the amplitude distribution in (3.1) as
\begin{equation} \{\hat{\eta},\hat{\phi}\}(\boldsymbol{x})=\{\hat{\eta},\hat{\phi}\}_{\boldsymbol{k}_{1}}\varPsi_{\boldsymbol{k}_{1}}+\{\hat{\eta},\hat{\phi}\}_{\boldsymbol{k}_{2}}\varPsi_{\boldsymbol{k}_{2}}, \end{equation}
without assuming 
$\beta$ to be small, we obtain the following second-order matrix eigenvalue problem for 
$\omega ^{2}$:
\begin{equation} \begin{pmatrix} (\omega^{2}-\boldsymbol{k}_{1}^{2})\boldsymbol{k}_{1}^{2} & \beta F_{\boldsymbol{k}_{1},\boldsymbol{k}_{2}}c_{n-1}\\ -\beta F_{\boldsymbol{k}_{1},\boldsymbol{k}_{2}}c_{m-1} & (\omega^{2}-\boldsymbol{k}_{2}^{2})\boldsymbol{k}_{2}^{2}\end{pmatrix}\begin{pmatrix} \hat{\phi}_{\boldsymbol{k}_{1}}\\ \hat{\phi}_{\boldsymbol{k}_{2}} \end{pmatrix}=\boldsymbol{\mathsf{0}}, \end{equation}
where 
$F_{\boldsymbol {k}_{1},\boldsymbol {k}_{2}}=2(m+n)(mn+1)/(mn)$. As before, for the system to be stable, the eigenvalue 
$\omega ^{2}$ has to be real. This is the case if 
$\beta \le \beta _{m,n}(\alpha ^{2})$, where
\begin{align} \beta_{m,n}(\alpha^{2}) &= \frac{{\rm \pi}^{4}}{16}\frac{\left|\alpha_{c}^{2}/\alpha{}^{2}-1\right|mn^{2}}{(m+n)(mn+1)}\nonumber\\ &\quad \times \left(\frac{(m^{2}-1)^{2}+\alpha^{4}(n^{2}-1)^{2}+2\alpha^{2}((m+n)^{2}+(mn+1)^{2})}{c_{m-1}c_{n-1}}\right)^{1/2}. \end{align}
For aspect ratios with 
$\alpha ^{2}=m$, where 
$m$ is an even number, using the two-mode approximation (3.32) for the preceding major critical point, i.e. 
$\alpha _{c}^{2}=m-1$, we have
\begin{equation} \beta_{m-1,1}(m)=\frac{{\rm \pi}^{4}}{16}\frac{(m-1)}{m^{2}}\sqrt{\frac{m^{2}+4}{c_{m-2}}}. \end{equation}
The lowest 
$\beta _{m-1,1}=({3{\rm \pi} ^{4}})/{128}\sqrt {{5}/{2}}\approx 3.61$ occurs at 
$m=4$, whereas the highest 
$\beta _{m-1,1}={{\rm \pi} ^{4}}/{16\sqrt {2}}\approx 4.30$ at 
$m=2$, which is then approached asymptotically at large even 
$m$.
The intersection of instability thresholds associated with the first two major critical points, which is defined by 
$\beta _{1,1}(\alpha ^{2})=\beta _{3,1}(\alpha ^{2})$, yields 
$\alpha _{c}^{2}\approx 2.125$ and, thus, in the two-mode approximation, we have 
$\beta _{c}\approx 4.699$. As shown in the next section, this is the optimal aspect ratio and the corresponding highest stability threshold which can be achieved with a small viscous friction.
3.2. Numerical solution of the matrix eigenvalue problem
For arbitrary 
$\beta$, (3.1) leads to
\begin{gather} \lambda\hat{\eta}+\nabla^{2}\hat{\eta}=0, \quad \left.\partial_{n}\hat{\eta}\right|_{{\varGamma}}=\beta\left.\partial_{\tau}\hat{\phi}\right|_{{\varGamma}}, \end{gather}
\begin{gather}\hat{\eta}+\nabla^{2}\hat{\phi}=0, \quad \left.\partial_{n}\hat{\phi}\right|_{{\varGamma}}=0, \end{gather}
which is an eigenvalue problem for 
$\lambda =\omega ^{2}+\mathrm {i}\gamma$, where 
$\gamma$ is a viscous friction coefficient. The problem can be discretised using the Galerkin method with the gravity wave modes (3.8) as basis functions. This leads to the generalisation of (3.14):
\begin{equation} (\lambda-\boldsymbol{k}^{2})\boldsymbol{k}^{2}\langle \varPsi_{\boldsymbol{k}}^{2}\rangle \hat{\phi}_{\boldsymbol{k}}=\beta\sum_{\boldsymbol{k}'}F_{\boldsymbol{k},\boldsymbol{k}'}\hat{\phi}_{\boldsymbol{k}'}, \end{equation}
with the electromagnetic interaction matrix on the right-hand side defined by (3.18). This is a matrix eigenvalue problem of size 
$(M+1)^{2}\times (N+1)^{2}$, where 
$M$ and 
$N$ are the cutoff limits of the 
$x$- and 
$y$-components of the wave vectors (3.9). This method with 
$M\times N=12$ modes was used by Sneyd & Wang (Reference Sneyd and Wang1994) to compute the stability threshold for a cell of length 
$[7.7]{m}$ depending on the width which varied from 
$[\approx 0.5]{m}$ to 
$[8]{m}.$ Alternatively, (3.34a,b) and (3.35a,b) can be discretised using the Chebyshev collocation method (Boyd Reference Boyd2013) (see Appendix B). As seen in figure 2, the latter has a significantly faster convergence rate than 
${\sim }N^{-4}$ achieved by the Galerkin approximation with 
$N$ modes in each direction. As the relative accuracy of the Chebyshev collocation approximation saturates at 
$\approx 10^{-11}\text {--}10^{-12}$ when 
$N {\mathop > \limits_\sim} 24$, in the following, we use 
$16\text {--}24$ collocation points in each direction.
The relative variation of the complex frequency with the largest imaginary part (
$\omega =3.14455+\mathrm {i}0.12889)$ versus the number of modes and nodes used in the Galerkin and Chebyshev approximations for 
$\alpha =\beta =1$.
The highest growth rate 
$\omega _{i}=\mathrm {Im}[\omega ]$, which is computed using Chebyshev collocation method with 
$N_{x}=N_{y}=16$ points and plotted in figure 3 against the interaction parameter 
$\beta$ for various aspect ratios 
$\alpha$, confirms the eigenvalue perturbation solution obtained in the previous section. Namely, for 
$\alpha ^{2}$ equal to the ratio of two odd numbers, the growth rate becomes positive at 
$\beta >\beta _{c}=0$ whereas for other aspect ratios this happens only at finite 
$\beta _{c}$. The dependence of the instability threshold 
$\beta _{c}$ on 
$\alpha ^{2}$ is shown in figure 4 for various viscous friction coefficients 
$\gamma$. Note that the instability threshold cannot be computed for 
$\gamma =0$ because it becomes absolutely fractal in this limit. Due to the finite accuracy of numerical solution, the lowest friction coefficient for which the stability threshold can reliably be computed is around 
$\gamma =10^{-5}.$ For this value, the stability diagram is very rugged and exhibits both small- and large-scale patterns. The key feature are the dips in 
$\beta _{c}$ which can be seen to occur at 
$\alpha _{c}^{2}$ equal to the ratio of two odd numbers as predicted by the eigenvalue perturbation analysis in the previous section. The increase of the friction coefficient 
$\gamma$ gradually smooths out the dependence of 
$\beta _{c}$ on 
$\alpha$, especially at small scales and larger aspect ratios. However, the main feature of the stability diagram, which is the location of minima and maxima of 
$\beta _{c}$ in the vicinity of odd and even values of 
$\alpha ^{2}$, respectively, persists up to relatively large friction coefficients 
$\gamma \sim 0.1$. It is remarkable that the approximate solution (3.32) for the major critical points closely reproduces the numerical results obtained with 
$N_{x}=N_{y}=24$ collocation points. According to the approximate solution obtained in the previous section, the first four maxima of 
$\beta _{c}\approx 4.70,4.11,4.19,4.23$ are located at 
$\alpha ^{2}\approx 2.13,4.16,6.12,8.10$.
The highest growth rate 
$\mathrm {Im}[\omega ]$ depending on the interaction parameter 
$\beta$ for various aspect ratios 
$\alpha$ computed by the Chebyshev collocation method with 
$N_{x}=N_{y}=16$. For 
$\alpha ^{2}=m/n$, where 
$m$ and 
$n$ are odd numbers, numerical results are compared with the approximate analytical solution (3.28).
Instability threshold 
$\beta _{c}$ depending on the aspect ratio squared 
$(\alpha ^{2})$ computed for various viscous friction coefficients 
$\gamma$ using Chebyshev collocation method with 
$N_{x}=N_{y}=16\text {--}24$ points. Analytical solution (3.32) is plotted for the dominant critical points 
$\alpha _{c}^{2}$ equal to odd numbers as well as their thirds and fifths.
4. Summary and discussion
We have carried out a comprehensive linear stability analysis of interfacial waves in the idealised model of an aluminium reduction cell. The model consists of two stably stratified liquid layers carrying a vertical electric current in the presence of a uniform collinear external magnetic field. Using the long-wave approximation, we derived a fully nonlinear shallow-water model comprising a novel integro-differential equation for the electric potential. Our shallow-water model differs from that of Bojarevics (Reference Bojarevics1998) by being fully nonlinear and, thus, applicable not only to small but also to large amplitude waves. It should also be noted that the electric potential equation we derived differs from that postulated by Zikanov et al. (Reference Zikanov, Thess, Davidson and Ziegler2000). The nonlinear potential equation is just a side result which may be useful for numerical modelling of large amplitude shallow-water waves, but it is not essential for the linear stability analysis of infinitesimal amplitude waves considered in the present study.
In this study, we showed that rectangular cells are inherently unstable if their aspect ratio squared equals the ratio of two odd numbers: 
$\alpha ^{2}=m/n.$ In the inviscid limit, such cells can be destabilised by an infinitesimally weak electromagnetic interaction while cells with other aspect ratios have finite instability thresholds. The unstable aspect ratios form an absolutely discontinuous dense set of points which intersperse aspect ratios with finite stability thresholds. This inherent fractality of the instability threshold was revealed analytically using an eigenvalue perturbation method and confirmed by accurate numerical solution of the underlying linear stability problem.
Although the instability threshold is gradually smoothed out by viscous friction, its principal structure, which is dominated by major critical aspect ratios corresponding to moderate values of 
$m$ and 
$n$, is well preserved up to relatively large dimensionless linear damping coefficients 
$\gamma \sim 0.1$. For 
$\gamma \lesssim 1,$ cells with such aspect ratios have the lowest stability threshold with the critical electromagnetic parameter 
$\beta _{c}\sim \gamma ,$ while the most stable are cells with 
$\alpha ^{2}\approx 2.13$ which have 
$\beta _{c}\approx 4.7$.
For a reduction cell operated at a total electric current of 
$I_{0}\approx [300]{kA}$ with the electrolyte and aluminium layers of depth 
$\bar {h}_{-}\approx [5]{cm}$ and 
$\bar {h}_{+}\approx [25]{cm}$, respectively, and the density difference 
$\Delta \rho =\rho _{+}-\rho _{-}\approx 200\ \textrm {kg}\,\textrm {m}^{-3}$, this 
$\beta _{c}$ corresponds to the critical magnetic field
\begin{equation} B_{0}=\beta_{c}g\Delta\rho\bar{h}_{+}\bar{h}_{-}/I_{0}\approx[0.27]{mT}. \end{equation}This field is an order of magnitude stronger than Earth's magnetic field but still an order of magnitude weaker than the magnetic fields which can be generated by the electric current flowing in the reduction cell itself and in the adjacent cells. It should be noted that it is not the whole magnetic field but only its vertical component which is important for the metal pad instability.
It is also important to note that the above result is based on the assumption of a relatively small linear damping in the two-layer system, i.e. 
$\gamma \sim \tau _{0}/\tau _{\nu }\lesssim 0.1$, where 
$\tau _{0}$ and 
$\tau _{\nu }$ are the characteristic wave propagation and viscous damping times with the latter being the reciprocal of the linear friction coefficient. For a cell with 
$L_{x}L_{y}\approx [40]{m^{2}}$ and 
$c\approx 0.2\ \textrm {m}\,\textrm {s}^{-1}$, we have 
$\tau _{0}=L/c\sim [10]{s.}$ Viscous damping is expected to be dominated by the electrolyte layer, which is thinner and also more viscous than the aluminium layer. Using 
$\nu _{-}\sim 10^{-6}\ \textrm {m}^{2}\,\textrm {s}^{-1}$ as a characteristic kinematic viscosity of cryolite (Hertzberg, Tørklep & Øye Reference Hertzberg, Tørklep and Øye2013), we obtain 
$\tau _{\nu }\sim \bar {h}_{-}^{2}/\nu _{-}\sim [10^{3}]{s}$. On the other hand, assuming 
$B_{0}\sim [1]{mT,}$we obtain the same order-of-magnitude magnetic damping time in the aluminium layer: 
$\tau _{m}\sim \rho _{+}/\sigma _{+}B_{0}^{2}\sim [10^{3}]{s}$. Consequently, we have 
$\gamma \sim 10^{-2}$ for both the viscous and the magnetic damping. This, however, is at odds with the customary assumption that the viscous damping time in real reduction cells is 
$\tau _{\nu }\sim [10]{s}$ and, thus, 
$\gamma \sim 1$ (Moreau & Ziegler Reference Moreau and Ziegler1988; Zikanov et al. Reference Zikanov, Thess, Davidson and Ziegler2000; Molokov et al. Reference Molokov, El and Lukyanov2011).
Such a high damping may be due to turbulence which is likely to develop when the waves become sufficiently large. Damping can also be affected by the large-scale circulation which occurs in real cells due to the non-uniformity of electric current and magnetic field. The same non-uniformity can also affect the instability threshold directly by modifying the electromagnetic coupling strength of different gravity wave modes. These, however, are rather complex effects which are outside the scope of the basic Sele instability model. Nevertheless, this model is still relevant to the stability of real aluminium reduction cells as confirmed by complete 3-D MHD simulations (Gerbeau et al. Reference Gerbeau, Le Bris and Lelièvre2006, § 6.1.4.3 and references therein).
When the background circulation is negligible, which is usually assumed to be the case in this model, the waves can have a sufficiently small initial amplitude for the linear viscous friction to be relevant. In the idealised model with uniform vertical current and a strictly collinear background magnetic adopted in this study, turbulent damping is a nonlinear effect which can be modelled with a velocity-dependent turbulent friction coefficient (Bojarevics & Pericleous Reference Bojarevics and Pericleous2006). Although this effect can potentially halt the growth of instability once the waves reach a certain amplitude, it is not expected to affect the onset of instability. Moreover, if the saturation amplitude is relatively small, the waves may remain practically unnoticeable well above the linear stability threshold. Nonlinear evolution and, especially, the saturation of metal pad instability are open questions which can be addressed using the fully nonlinear shallow-water model proposed in this paper. This, however, is outside the scope of the present paper which is focused on the linear stability of a two-layer system. It is also relatively straightforward to extend this study to three-layer systems such as liquid metal batteries, which will be considered elsewhere.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Linear stability of semi-infinite and channel geometries
Since the electromagnetic interaction parameter appears only in the boundary condition (2.38), the destabilisation of small-amplitude interfacial waves in the presence of a uniform vertical magnetic field is essentially a boundary effect (Lukyanov et al. Reference Lukyanov, El and Molokov2001). It means that lateral boundaries are critical for the metal pad instability. As pointed out by Lukyanov et al. (Reference Lukyanov, El and Molokov2001), this fact is obscured by the weak mathematical formulation of the problem where boundary conditions are represented by integral volumetric terms (Bojarevics & Romerio Reference Bojarevics and Romerio1994). It is also missed by the electromechanical pendulum models where the liquid metal layer is substituted by a solid slab (Davidson & Lindsay Reference Davidson and Lindsay1998; Zikanov Reference Zikanov2015). The simplest hydrodynamic model of the metal pad instability is provided by a semi-infinite system 
$(0 < y < \infty )\times (-\infty < x < \infty )$ bounded by a single lateral wall at 
$y=0$.
Let us briefly review this elementary model which was first considered by Lukyanov et al. (Reference Lukyanov, El and Molokov2001) and later revisited by Molokov et al. (Reference Molokov, El and Lukyanov2011). Owing to the 
$x$-invariance and stationarity of the base state, a small-amplitude disturbance of the interface height 
$\eta$ and the associated electric potential 
$\phi$ can be sought as the normal mode:
\begin{equation} \{\eta,\phi\}(\boldsymbol{x},t)=\{\hat{\eta},\hat{\phi}\}(y)\textrm{e}^{\textrm{i}(kx-\omega t)}, \end{equation}
where 
$k$ is a real wavenumber describing harmonic variation of coupled interface and potential perturbation in the longitudinal 
$(x)$ direction, 
$\omega$ is a generally complex frequency, and 
$\hat {\eta }(y)$ and 
$\hat {\varphi }(y)$ are amplitude distributions in the transverse 
$(y)$ direction. Then (2.36) and (2.37) for the perturbation amplitudes take the form
\begin{gather} {-}\omega^{2}\hat{\eta} = \hat{\eta}''-k^{2}\hat{\eta}, \end{gather}
\begin{gather} {-}\hat{\eta} = \hat{\phi}''-k^{2}\hat{\phi}, \end{gather}
while the boundary conditions at 
$y=0$ read as
\begin{equation} \hat{\eta}'(y)-\textrm{i}k\beta\hat{\phi}(y)=\hat{\phi}'(y)=0. \end{equation}
The solution of (A2) and (A3) bounded at 
$y\rightarrow \infty$ can be written in the general form as
\begin{gather} \hat{\eta}(y) = \hat{\eta}_{-}\textrm{e}^{-\textrm{i}\kappa y}+\hat{\eta}_{+}\textrm{e}^{\textrm{i}\kappa y}, \end{gather}
\begin{gather}\hat{\phi}(y) = \hat{\phi}_{0}\textrm{e}^{{-}ky}+\omega^{{-}2}\hat{\eta}(y), \end{gather}
where 
$\hat {\eta }_{\pm }$ and 
$\hat {\phi }_{0}$ are unknown constants and 
$\kappa =\sqrt {\omega ^{2}-k^{2}}$ is the wavenumber of perturbation in the transverse 
$(y)$ direction. There are two types of solutions possible. The first corresponds to a real 
$\kappa$ and describes pure gravity waves with real frequency 
$\omega =\pm \sqrt {k^{2}+\kappa ^{2}}$. For this solution, boundary conditions (A4) yield
\begin{equation} \frac{\hat{\eta}_{+}-\hat{\eta}_{-}}{\hat{\eta}_{+}+\hat{\eta}_{-}} =\frac{k\beta}{\kappa(\omega^{2}-\textrm{i}\beta)},\end{equation}
which defines the ratio of amplitudes forming the eigenmode (A1). As seen, the electromagnetic effect breaks the reflection symmetry, 
$\hat {\eta }_{+}=\hat {\eta }_{-}$, which holds in the non-magnetic case 
$(\beta =0)$. With 
$\beta \not =0$, the amplitude and phase of the reflected wave are no longer the same as those of the incident wave. However, the amplification of the reflected wave, which is possible in this case, does not mean instability as suggested by Lukyanov et al. (Reference Lukyanov, El and Molokov2001). The stability of the system is defined by the frequency of the eigenmode rather than by its amplitude distribution. Namely, the eigenmode is neutrally stable as long as the frequency is real. Also, note that for real frequencies, there is no physical distinction between the incident and reflected waves which can be swapped one with the other because of the time inversion symmetry. Therefore the amplification of the reflected wave is not the mechanism behind the metal pad instability.
There is another solution branch which explicitly describes interfacial wave instability. This genuinely unstable mode, which was missed by Lukyanov et al. (Reference Lukyanov, El and Molokov2001) and considered first by Morris & Davidson (Reference Morris and Davidson2003), is defined by the complex frequencies 
$\omega$. In this case, the general solution of (A2) decaying at 
$y\rightarrow \infty$ can be written as
\begin{equation} \hat{\eta}(y)=\hat{\eta}_{0}\textrm{e}^{{-}y\kappa}, \end{equation}
where 
$\kappa =\sqrt {k^{2}-\omega ^{2}}$ represents the branch with positive real part. Substituting the solution above along with (A6) into the boundary conditions (A4), after a few rearrangements, we obtain the dispersion relation 
$\kappa (\kappa +k)+i\beta =0$, which yields
\begin{equation} \kappa=\sqrt{\frac{k^{2}}{4}-\textrm{i}\beta}-\frac{k}{2}. \end{equation}The respective frequency 
$\omega =\pm \sqrt {k^{2}-\kappa ^{2}}$ then follows from the definition of 
$\kappa$ given below (A8). The complexity of 
$\kappa$ for 
$\beta \not =0$ implies that one of the two frequencies has a positive imaginary part which, in turn, means that the respective mode is unstable at however small 
$\beta$. For 
$\beta /k^{2}=\varepsilon \ll 1$, which is close to the instability threshold and opposite to the limit of 
$\beta \gg 1$ producing the so-called edge waves (Morris & Davidson Reference Morris and Davidson2003), we have 
$\kappa \approx -\textrm {i}k(\varepsilon +\textrm {i}\varepsilon ^{2})$ and 
$\omega \approx \pm k(1+\textrm {i}\varepsilon ^{3})$. These relations describe a slightly oblique wave with the transverse wavenumber 
$\beta /k\ll k$ which decays over the distance 
$\sim k^{3}/\beta ^{2}$ from the wall and grows exponentially in time at the rate 
$\sim \beta ^{3}/k^{5}$. Namely, in this model, the interface can be destabilised by an arbitrary weak electromagnetic effect which gives rise to a nearly transverse wave travelling along the wall.
It is instructive to contrast the instability in the semi-infinite model considered above with that in the slightly more realistic model of a finite-width channel. This model was first considered by Davidson & Lindsay (Reference Davidson and Lindsay1998) and further studied in the context of hydromagnetic edge waves by Morris & Davidson (Reference Morris and Davidson2003) and Molokov et al. (Reference Molokov, El and Lukyanov2011).
For a channel bounded by two lateral walls at 
$y=\pm 1$, the general solution of (A2) and (A3) can be written as
\begin{gather} \hat{\eta}(y) = \hat{\eta}_{-}\sin(\kappa y)+\hat{\eta}_{+}\cos(\kappa y), \end{gather}
\begin{gather}\hat{\phi}(y) = \hat{\phi}_{-}\sinh(ky)+\hat{\phi}_{+}\cosh(ky)+\omega^{{-}2}\hat{\eta}(y), \end{gather}
where 
$\kappa$ is the transverse wavenumber satisfying the dispersion relation of pure interfacial gravity waves
\begin{equation} \omega^{2}=k^{2}+\kappa^{2}. \end{equation} Applying boundary conditions (A4) at 
$y=\pm 1$, after a few rearrangements, we obtain
\begin{equation} \beta(k,\kappa)={\pm}\left(k^{2}+\kappa^{2}\right)\left(\frac{k}{\kappa}\left(\frac{k}{\kappa}+\frac{\tan\kappa}{\tanh k}-\frac{\tanh k}{\tan\kappa}\right)-1\right)^{{-}1/2}, \end{equation}
which explicitly defines 
$\beta$ depending on the longitudinal and transverse wavenumbers 
$k$ and 
$\kappa .$ It is important to note that the respective eigenmode is a pure gravity wave as long as 
$\kappa$ is real. On the other hand, (A13) implicitly defines the spectrum of admitted 
$\kappa$ values for given 
$k$ and 
$\beta$. For 
$\beta =0$, this equation reduces to 
$\sin 2\kappa =0$, which describes the standard wavenumber spectrum for a channel of dimensionless width equal to 2: 
$\kappa _{n}=n{\rm \pi} /2, n=0,1,2,\ldots$. As seen in figure 5, where 
$\beta (k,\kappa )$ is plotted against 
$\kappa$ for various 
$k$, the increase of 
$\beta$ just modifies the spectrum of the wave modes admitted by the electromagnetic reflection condition (2.38). It is important to note that the eigenmodes remain pure gravity waves satisfying the dispersion relation (A12) with real wavenumbers up to the point where two branches of 
$\kappa$ merge together. The absence of coupling between the gravity wave modes below the critical point is due to the effective reduction of the electromagnetic force to the boundary condition (2.38) which affects only the reflection but not the propagation of gravity waves. At the critical point, a complex conjugate pair of wavenumbers emerges, which means instability as the frequency becomes complex (see figure 6). As found by Morris & Davidson (Reference Morris and Davidson2003), who use the reciprocal of 
$\beta$, this happens first at 
$\beta _{c}\approx 1.365$ when two purely longitudinal 
$(k=0)$ gravity wave modes with 
$\kappa _{c}=1.113$ merge.
The electromagnetic interaction parameter 
$\beta$ versus the transverse wavenumber 
$\kappa$ for various longitudinal wavenumbers 
$k$.
(a) Real 
$(\mathrm {Re})$ and imaginary 
$(\mathrm {Im})$ parts of the first three transverse wavenumbers 
$\kappa$, (b) the respective frequencies 
$\omega$ versus the electromagnetic interaction parameter 
$\beta$ for various longitudinal wavenumbers 
$k$ and (c) the critical 
$\beta$ and transverse wavenumbers 
$\kappa$ as functions of the longitudinal wavenumbers 
$k.$
Note that this critical mode is longitudinal and emerges at finite 
$\beta _{c}$, whereas the one predicted by the semi-infinite model is transverse and emerges at 
$\beta _{c}=0$. In the channel, the latter mode is recovered in the short-wave limit 
$k\gg 1$. For such waves, the marginal interaction parameter and the respective transverse wavenumber can be seen in figure 6(c) to scale as 
$\beta \sim k$ and 
$\kappa \sim 1$. The interaction parameter based on the wavelength, which is the relevant horizontal length scale in this limit, then scales as 
$\tilde{{\beta }}=\beta /k^{2}\sim k^{-1}\rightarrow 0$ so recovering the result of the semi-infinite model.
Appendix B. Chebyshev collocation method
To solve the eigenvalue problem posed by (3.34a,b) and (3.35a,b) numerically, we use a collocation method with the two-dimensional Chebyshev–Lobatto grid:
\begin{align} (x_{m},y_{n})=\frac{1}{2}\left(\sqrt{\alpha}\cos\frac{{\rm \pi} m}{M+1},\frac{1}{\sqrt{\alpha}}\cos\frac{{\rm \pi} n}{N+1}\right),(0,0)\le(m,n)\le(M+1,N+1), \end{align}
on which the discretised solution 
$\{\hat {\eta },\hat {\phi }\}(x_{m},y_{n})=\{\boldsymbol {\eta },\boldsymbol {\phi }\}_{m,n}$ and its derivatives are sought. The latter are expressed in terms of the former by substituting the first and second derivatives with the differentiation matrices (Peyret Reference Peyret2002, 393–394). Requiring the equations to be satisfied at the inner collocation points 
$(1,1)\le (m,n)\le (M,N)$, we obtain a system of discrete equations which can be written as
\begin{gather} {-}\lambda\boldsymbol{\eta}_{0} = \boldsymbol{\mathsf{A}}_{0}\boldsymbol{\eta}_{0}+\boldsymbol{\mathsf{A}}_{1}\boldsymbol{\eta}_{1}, \end{gather}
\begin{gather} {-}\boldsymbol{\eta}_{0} = \boldsymbol{\mathsf{A}}_{0}\boldsymbol{\phi}_{0}+\boldsymbol{\mathsf{A}}_{1}\boldsymbol{\phi}_{1}, \end{gather}
where the subscripts 
$0$ and 
$1$ denote the parts of the solution at the inner and boundary collocation points, respectively. The matrices 
$\boldsymbol{\mathsf{A}}_{0}$ and 
$\boldsymbol{\mathsf{A}}_{1}$ represent parts of the discretised Laplacian using the respective sets of points. The boundary conditions are applied at the boundary points 
$m=0,M+1$ and 
$n=0,N+1$, which results in
\begin{gather} \boldsymbol{\mathsf{B}}_{0}\boldsymbol{\eta}_{0}+\boldsymbol{\mathsf{B}}_{1}\boldsymbol{\eta}_{1} = \beta\boldsymbol{\mathsf{C}}\boldsymbol{\phi}_{1}, \end{gather}
\begin{gather}\boldsymbol{\mathsf{B}}_{0}\boldsymbol{\phi}_{0}+\boldsymbol{\mathsf{B}}_{1}\boldsymbol{\phi}_{1} = \boldsymbol{0}, \end{gather}
where 
$\boldsymbol{\mathsf{B}}$ and 
$\boldsymbol{\mathsf{C}}$ are the discretised normal and tangential derivative operators at the boundary. The above system of discrete equations is reduced to a matrix eigenvalue problem as follows. First, (B4) and (B5) can solved for 
$\boldsymbol {\phi }_{1}$ and 
$\boldsymbol {\eta }_{1}$ in terms of 
$\boldsymbol {\phi }_{0}$ and 
$\boldsymbol {\eta }_{0}$ as
\begin{gather} \boldsymbol{\phi}_{1} ={-}\boldsymbol{\mathsf{B}}_{1}^{{-}1}\boldsymbol{\mathsf{B}}_{0}\boldsymbol{\phi}_{0}, \end{gather}
\begin{gather}\boldsymbol{\eta}_{1} ={-}\boldsymbol{\mathsf{B}}_{1}^{{-}1}\boldsymbol{\mathsf{B}}_{0}\boldsymbol{\eta}_{0}+\beta\boldsymbol{\mathsf{B}}_{1}^{{-}1}\boldsymbol{\mathsf{C}}\boldsymbol{\phi}_{1}. \end{gather}
Substituting (B6) into (B3), we have 
$\boldsymbol {\eta }_{0}=-\widetilde {\boldsymbol{\mathsf{A}}}_{0}\boldsymbol {\phi }_{0}$, where 
$\widetilde {\boldsymbol{\mathsf{A}}}_{0}=\boldsymbol{\mathsf{A}}_{0}-\boldsymbol{\mathsf{A}}_{1}\boldsymbol{\mathsf{B}}_{1}^{-1}\boldsymbol{\mathsf{B}}_{0}$. Combining this and (B6) with (B7) and (B2), after a few rearrangements, we obtain the sought matrix eigenvalue problem
\begin{equation} \lambda\boldsymbol{\phi}_{0}=[-\widetilde{\boldsymbol{\mathsf{A}}}_{0} +\beta\widetilde{\boldsymbol{\mathsf{A}}}_{0}^{{-}1}\boldsymbol{\mathsf{A}}_{1}\boldsymbol{\mathsf{B}}_{1}^{{-}1} \boldsymbol{\mathsf{C}}\boldsymbol{\mathsf{B}}_{1}^{{-}1}\boldsymbol{\mathsf{B}}_{0}]\boldsymbol{\phi}_{0}, \end{equation}which is solved using the standard DGEEV routine of LAPACK.
