# Electromagnetically driven flow in unsupported electrolyte layers:
lubrication theory and linear stability of annular flow

Andrey Pototsky<EMAIL_ADDRESS>Swinburne University of Technology,
Department of Mathematics, Hawthorn, 3122 Victoria, Australia  Sergey A.
Suslov  Swinburne University of Technology, Department of Mathematics,
Hawthorn, 3122 Victoria, Australia

###### Abstract

We consider a thin horizontal layer of a non-magnetic electrolyte containing a
bulk solution of salt and carrying an electric current. The layer is bounded
by two deformable free surfaces loaded with an insoluble surfactant and is
placed in a vertical magnetic field. The arising Lorentz force drives the
electrolyte in the plane of the layer. We employ the long-wave approximation
to derive general two-dimensional hydrodynamic equations describing symmetric
pinching-type deformations of the free surfaces. These equations are used to
study the azimuthal flow in an annular film spanning the gap between two
coaxial cylindrical electrodes. In weakly deformed films, the base azimuthal
flow and its linear stability with respect to azimuthally invariant
perturbations are studied analytically. For relatively thick layers and weak
magnetic fields, the leading mode with the smallest decay rate is found to
correspond to a monotonic azimuthal velocity perturbation. The Marangoni
effect leads to further stabilisation of the flow while perturbations of the
solute concentration in the bulk of the fluid have no influence on the flow
stability. In strongly deformed films in the diffusion-dominated regime, the
azimuthal flow becomes linearly unstable with respect to an oscillatory mixed
mode characterised by the combination of radial and azimuthal velocity
perturbations when the voltage applied between electrodes exceeds the critical
value.

interfacial flows (free surface)

###### pacs:

##  I Introduction

Electromagnetically driven flows in shallow layers and channels of
electrically conducting fluids in the presence of deformable interfaces
attracted much attention due to their importance in plasma physics Fiflis et
al. (2016); Lunz & Howell (2019) and various microfluidic applications
including contactless manipulation of flow in magnetohydrodynamic networks
(Bau et al., 2003), liquid channels embedded into carrier fluids (Dunne et
al., 2020), droplet microfluidics (Shang et al., 2017) and electromagnetic
stirring (Bau et al., 2001; Qian & Bau, 2005).

The comprehensive theoretical description of magnetohydrodynamic (MHD) flows
of magnetic fluid films in arbitrary strong magnetic fields is technically
challenging as the hydrodynamic equations must be coupled with Maxwell’s
equations in the presence of deformable moving boundaries. However, in the
case of non-magnetic electrically conducting fluids, the description can be
greatly simplified. For this class of fluids, which comprises electrolyte
solutions and liquid metals with weak magnetic properties, the additional
stresses that typically appear at the interfaces due to externally applied
magnetic fields can be neglected. For relatively weak magnetic fields of the
order of $10^{-2}-10^{0}$ T, which can be created using conventional permanent
magnets, the magnetic Reynolds number $\,\textit{Re}_{m}=UL/\eta_{m}$
associated with the flow of fluid with the magnetic diffusivity $\eta_{m}$ in
the domain of a characteristic size $L$ with velocity $U$, is typically
small,$\,\textit{Re}_{m}\ll 1$. In this regime, the magnetic field induced by
the electric current flowing through the fluid can be neglected compared to
the external field (Müller & Bühler, 2013). With this simplification, the MHD
equations have been successfully applied to study flows of electrolitic
solutions and non-magnetic liquid metals such as mercury in various
geometries. These include liquid metal layers confined between two parallel
insulating walls (Sommeria & Moreau, 1982) and in thin horizontal films
(Sommeria, 1986) and electrolyte solutions in annular channels (Messadek &
Moreau, 2002; Figueroa et al., 2009; Pérez-Barrera et al., 2016; Suslov et
al., 2017; McCloughan & Suslov, 2020).

Geometric parameters of the system such as the depth and the aspect ration of
the layer have been shown to have the major influence on the flow
characteristics. In shallow horizontal layers of electrolyte solutions with a
depth of several milimetres and a small depth-to-width aspect ratio placed
between two coaxial vertical electrodes the flow was found to be essentially
three-dimensional even for relatively weak currents (Figueroa et al., 2009;
Pérez-Barrera et al., 2016; Suslov et al., 2017; McCloughan & Suslov, 2020).
The quasi-two-dimensional approximation, developed in Figueroa et al. (2009);
Pérez-Barrera et al. (2016) by using the depth-averaging method, could capture
some of the main features of the base azimuthal flow but was shown to be
inadequate when describing toroidal flows that lead to the formation of the
experimentally observable free-surface vortices Suslov et al. (2017).

As the depth of a horizontal layer and the aspect ratio of the system are
further decreased, the vertical component of the flow velocity is impeded by
the boundaries and the horizontal component of the flow becomes dominant. The
two-dimensional nature of the flow in very thin liquid layers was used to
study two-dimensional turbulence as pioneered around four decades ago by
Couder (1981, 1984) in experiments with soap films that were mechanically
stirred by the array of rods to produce turbulent flow. At around the same
time Sommeria (1986) studied effectively two-dimensional flows in thin mercury
films with a free upper surface and supported from below by an array of
conducting electrodes. Instead of a mechanical stirring, a contactless
electromagnetic Lorentz forcing was used to drive the flow in the presence of
highly non-uniform magnetic fields. Later, similar contactless electromagnetic
forcing was used to gain a deeper understanding of the scaling properties of
the velocity correlation function in turbulent regimes (Cardoso et al., 1994;
Marteau et al., 1995; Williams et al., 1997). A comprehensive review of two-
dimensional turbulence can be found in (Kellay & Goldburg, 2002, e.g.).

Historically, electromagnetic driving was extensively used in supported liquid
layers, but not in unsupported systems such as soap films. This, perhaps, was
due to the intrinsic instability of soap-type films and the difficulty of
controlling their curvature. In fact, to the best of our knowledge, the first
attempt to use Lorentz force in unsupported free films was made almost 20
years after Couder (1981, 1984) pioneered the film turbulence studies. It used
a soap film containing chloride salt spanning a region between two parallel
conducting electrodes placed above an array of permanent magnets (Rivera & Wu,
2000). The main advantage of using an unsupported film when studying two-
dimensional turbulence is the elimination of energy leakage at the no-slip
bottom of the container. In a recent experimental study Cruz Gómez (2016)
investigated the flow dynamics in an electromagnetically forced film in the
case of a localised source of electric current. A series of experiments with
soap films spanning the gap between two coaxial electrodes placed in an
external magnetic field is currently underway Cuevas & Figueroa (2023).

The theoretical description and modelling of electromagnetically driven flows
in supported films with a free interface are now well developed (Morley &
Abdou, 1995; Morley & Roberts, 1996; Morley & Abdou, 1997; Morley et al.,
2002; Gao & Morley, 2002; Gao et al., 2002; Miloshevsky & Hassanein, 2010;
Giannakis et al., 2009a, b; Lunz & Howell, 2019) and continue to attract
attention mainly due to applications in plasma flows and tokamaks. In the
absence of the Lorentz force, the pure hydrodynamic description of the flow in
unsupported liquid films was initiated in Prévost & Gallez (1986); Sharma &
Ruckenstein (1988) and later received a huge boost because of its relevance to
nonlinear film rupture and two-dimensional turbulence problems (Couder et al.,
1989; Chomaz & Cathalau, 1990; Erneux & Davis, 1993; Gharib & Derango, 1989;
Sharma et al., 1995; Van De Fliert et al., 1995; Wu et al., 1995) as reviewed
in (Kellay & Goldburg, 2002; Oron et al., 1997). In unsupported thin viscous
films, the hydrodynamic equations are simplified using two main assumptions.
Firstly, the long-wave approximation is applied by taking into account large-
scale flow patterns and film deformations the wavelength of which is much
larger than the average film thickness. Secondly, deformations of the free
interfaces are assumed to be mirror-symmetric with respect to the centre plane
of the layer, which corresponds to a varicose-type pinching deformation mode.
Under these assumptions, the effective two-dimensional dynamic equations were
derived at the leading order of the lubrication approximation for curved soap
films in the presence of a surfactant (Ida & Miksis, 1998; Miksis & Ida, 1998)
and for horizontally stretched free films in the presence of solute and
surfactant (Chomaz, 2001). So far, the application of the lubrication
approximation to describe MHD flows in unsupported free films of electrolyte
solutions in external magnetic fields has not been reported.

Here we build upon earlier theoretical studies (Ida & Miksis, 1998; Miksis &
Ida, 1998; Chomaz, 2001) to derive the leading order dynamic equations for an
electromagnetically driven two-dimensional flow in a thin free horizontal
layer of electrolyte solution the free surfaces of which are loaded with an
insoluble surfactant. The flow is driven by the Lorentz force generated by the
electric current flowing through the electrolyte in the presence of a
homogeneous external magnetic field normal to the layer. We show that at the
leading order of the lubrication approximation the product of the current
density and the local thickness of the layer is divergence-free reflecting the
condition of no accumulation of electric charge in the bulk. The complete set
of the derived dynamic equations is written in an invariant vector form
suitable for applications in arbitrary geometries. It consists of the dynamic
equation for the two-dimensional flow field depending on the solute
concentration in the bulk, surfactant concentration, symmetric film
deformations and the electric potential. As an example, we apply the derived
equations to study the azimuthal flow and its linear stability in annular free
films spanning the gap between two coaxial conducting electrodes.

The paper is organised as follows. In section II we present the derivation of
the leading order equations in the lubrication approximation using the
systematic expansion technique suggested earlier in Erneux & Davis (1993);
Chomaz (2001). The derived equations correctly reflect the conservation of the
total mass of the surfactant and solute as well as the continuity of electric
current under the condition of no accumulation of the electric charge in the
fluid. In section III we rewrite the derived equations in polar coordinates to
study the flow in annular free films. In section IV we study the linear
stability of the annular azimuthally invariant steady-state with respect to
perturbations that depend only on the radial coordinate. We present analytical
results for the linear stability of a flat film and vanishingly small flow
velocities first. Subsequently, we use the numerical continuation method
Krauskopf et al. (2014); Doedel et al. (1994) to study the stability of a
strongly deformed layer. The obtained theoretical and computational results
are summarised in section V.

##  II Lubrication theory of electrically conducting free films in an external
magnetic field

Consider a horizontal free film of an electrolyte solution, which can be
created by supporting the weight of the film by the pressure difference
between the regions below and above the film. Following Erneux & Davis (1993);
Chomaz (2001) we exclude film bending and only consider symmetric pinching-
type surface deformation modes $z=\pm h(x,y,t)$ with each surface being a
mirror image of the other at all times as schematically shown in figure 1(a).
Here, $x$ and $y$ are the coordinates in the horizontal plane, $z$ is a
vertical coordinate and $t$ is time. The local thickness of the film is $2h$
and the average film thickness is

$2\langle h\rangle=2S^{-1}\int_{S}h(x,y,t)\,\textrm{d}x\,\textrm{d}y\,,$ (1)

where $S$ denotes the area of the centre plane $z=0$. Each surface of the film
is loaded with an insoluble surfactant with local concentration $c_{s}$. The
addition of surfactants is particularly important in soap films that contain
fatty acid carboxylates, which are typically found at the surface. The
electrolyte solution is composed of a solvent fluid (typically pure water) and
dissociated salt molecules with bulk concentration $c_{b}$. In what follows we
assume that salt is completely soluble and does not form a molecular surface
layer. Two electrodes are immersed in the fluid so that the electric current
can flow between them through the film when external voltage is applied. One
may consider at least four possible topological configurations of a free film
spanning space between two electrodes as shown in figure 1($b$–$e$). The
surface of the electrodes $\partial\Sigma$ is assumed to be chemically inert
and impenetrable to the surfactant and the solute in the film. In addition,
the flow field ${\bm{u}}$ vanishes at $\partial\Sigma$.

Figure 1: ($a$) Symmetric deformation mode in a horizontal free liquid film
with two deformable surfaces located at $z=\pm h(x,y,t)$. The flow field
${\bm{u}}=(u,v,w)$ is mirror symmetric with respect to the centre-plane $z=0$.
($b$–$e$) The top view of the system: four possible topological configurations
of a free film spanning space between two electrodes $(1,2)$. The surface of
each electrode $\partial\Sigma$ represents a no-slip equipotential boundary
impenetrable for surfactant and electrolyte solution.

The electrical conductivity $\sigma$ of the electrolyte solution generally
depends on $c_{b}$. The solution density $\rho$ and dynamic viscosity $\mu$
(kinematic viscosity $\nu=\mu/\rho$) are assumed to be constant and
independent of $c_{b}$. The externally applied magnetic field
${\bm{B}}=(0,0,B(x,y))$ is assumed to be significantly stronger than that
induced by the motion of the fluid. In what follows we neglect gravity effects
anticipating that hydrostatic pressure in a sub-micrometre thin free film is
negligible as compared to the Laplace pressure. The motion of the
incompressible fluid with three-dimensional velocity ${\bm{u}}=(u,v,w)$ is
described by the continuity and Navier-Stokes equations with the added Lorentz
force term (Müller & Bühler, 2013)

$\displaystyle{\bm{\nabla}}\cdot{\bm{u}}$ $\displaystyle=$ $\displaystyle
0\,,$ (2)
$\displaystyle\partial_{t}{\bm{u}}+({\bm{u}}\cdot{\bm{\nabla}}){\bm{u}}$
$\displaystyle=$
$\displaystyle-\frac{1}{\rho}{\bm{\nabla}}(p-\Pi)+\frac{\mu}{\rho}\nabla^{2}{\bm{u}}+\frac{1}{\rho}{\bm{j}}\times{\bm{B}}\,,$
(3)

where $p$ is the pressure in the fluid and $\Pi=\Pi(h(x,y,t))$ represents the
disjoining pressure due to intermolecular forces that become important when
the film thickness is approximately $100$ nm or less Israelachvili (2011);
Overbeek (1960). Because of the symmetry of the pinching mode, the horizontal
flow velocity components $(u,v)$ and the vertical velocity component $w$
should be even and odd function of $z$, respectively.

The current density ${\bm{j}}=(j_{x},j_{y},j_{z})$ is related to the electric
potential $\phi$ via Ohm’s law

${\bm{j}}=\sigma(c_{b})(-{\bm{\nabla}}\phi+{\bm{u}}\times{\bm{B}})\,.$ (4)

For electrolyte solutions, we assume linear dependence between conductivity
$\sigma$ and $c_{b}$

$\sigma(c_{b})=Kc_{b}\,,$ (5)

where $K$ is an empirical constant specific to a particular salt and solvent.
Under the condition of no accumulation of electric charge in the bulk, the
electric potential $\phi$ is found from

${\bm{\nabla}}\cdot{\bm{j}}=0\,,$ (6)

which must be solved instantaneously for any given velocity field ${\bm{u}}$.
Additionally, the current continuity condition 6 must be supplemented with the
boundary conditions for $\phi$ that correspond to the equipotential surfaces
$\partial\Sigma$ of the conducting electrodes.

Note that for a magnetic field ${\bm{B}}=(0,0,B)$ orthogonal to the layer
${\bm{j}}\times{\bm{B}}=B(j_{y},-j_{x},0)$,
${\bm{j}}=\sigma(c_{b})(-\partial_{x}\phi+Bv,-\partial_{y}\phi-
Bu,-\partial_{z}\phi)$ and ${\bm{u}}\times{\bm{B}}=B(v,-u,0)$. The normal
component of the electric current vanishes at the film surfaces. Thus, at t
$z=h$ we require

$j_{n}={\bm{j}}\cdot{\bm{n}}=0\,,$ (7)

where
${\bm{n}}=(-\partial_{x}h,-\partial_{y}h,1)/\sqrt{1+(\partial_{x}h)^{2}+(\partial_{y}h)^{2}}$
is the unit normal vector to the upper film surface directed away from the
fluid.

At $z=h$, the kinematic boundary condition applies

$\partial_{t}h+({\bm{u}}\cdot{\bm{\nabla}}_{\parallel})h=w\,,$ (8)

where ${\bm{\nabla}}_{\parallel}=(\partial_{x},\partial_{y})$ is the
horizontal gradient. Note that 8 can also be written in the equivalent form as
$\partial_{t}h=\sqrt{1+({\bm{\nabla}}_{\parallel}h)^{2}}({\bm{n}}\cdot{\bm{u}})$.

To describe the dynamics of surfactant and the concentration of salt in the
bulk of the fluid, we follow Jensen & Grotberg (1993). The bulk concentration
$c_{b}$ is described by the advection-diffusion equation

$\partial_{t}c_{b}+{\bm{\nabla}}\cdot({\bm{u}}c_{b}-d_{b}{\bm{\nabla}}c_{b})=0\,,$
(9)

where $d_{b}$ is the bulk diffusion coefficient, ${\bm{u}}c_{b}$ is the
advective flux and $-d_{b}{\bm{\nabla}}c_{b}$ is the diffusive flux.

The advection-diffusion equation for the surfactant concentration at the upper
film surface $z=h(x,y,t)$ is given by

$\partial_{t}c_{s}+{\bm{\nabla}}_{s}\cdot(c_{s}{\bm{u}})=d_{s}{\bm{\nabla}}_{s}^{2}c_{s}\,,$
(10)

where $d_{s}$ is the surface diffusivity of the surfactant and
${\bm{\nabla}}_{s}={\bm{\nabla}}-{\bm{n}}({\bm{n}}\cdot{\bm{\nabla}})$ is the
surface gradient.

At $z=h$ the diffusive flux normal to the film surface must vanish

${\bm{n}}\cdot{\bm{\nabla}}c_{b}=0\,.$ (11)

It can be shown that 11 supplemented with the condition that the flow velocity
${\bm{u}}$ and the normal diffusive fluxes of surfactant and solute vanish at
the surface of the electrodes leads to the conservation of the total mass of
the solute and the surfactant.

Next, we consider the balance of the normal and tangential forces at the upper
surface $z=h(x,y,t)$

$(p+\Gamma\kappa){\bm{n}}={\bm{\nabla}}_{s}\Gamma+{\bm{T}}\cdot{\bm{n}}\,,$
(12)

where
${\bm{T}}=\mu[{\bm{\nabla}}\otimes{\bm{u}}+({\bm{\nabla}}\otimes{\bm{u}})^{T}]$
is the viscous stress tensor, $\kappa(x,y,t)$ is the local mean curvature of
the surface defined as $\kappa=-{\bm{\nabla}}\cdot{\bm{n}}$, $\Gamma(x,y,t)$
is the local surface tension, $\otimes$ is the tensor product and the
superscript $T$ denotes transposed quantities. In what follows we assume that
the liquid is non-magnetic and the applied magnetic field is relatively weak
so that the Maxwell component in the stress tensor can be completely
neglected.

The gradient of the surface tension along the interface $\partial_{s}\Gamma$
is induced by the distribution of the surfactant according to the soluto-
Marangoni effect

$\Gamma(x,y,t)=\gamma-\Gamma_{M}c_{s}(x,y,t)\,,$ (13)

where $\gamma$ is the reference surface tension in the absence of a surfactant
and $\Gamma_{M}=-d\Gamma/dc_{s}>0$ is assumed to be constant. The balance of
forces at the lower surface $z=-h(x,y,t)$ is automatically achieved for
symmetric deformation modes.

The horizontal length scale $L$ of the flow in submicrometre thin films is
several orders of magnitude larger than the average film thickness $2\langle
h\rangle$. The long-wave approximation theory of one-dimensional free films in
the absence of a surfactant and an electric current has been developed some
thirty years ago (Erneux & Davis, 1993) using systematic expansion of the
Navier-Stokes equations for small lubrication parameter $\epsilon=\langle
h\rangle/L\ll 1$. Subsequently, the theory was generalised to describe two-
dimensional flat and curved free films loaded with soluble and insoluble
surfactant agents (Ida & Miksis, 1998; Miksis & Ida, 1998; Chomaz, 2001). Here
we extend earlier results to derive the leading-order equations for long-wave
symmetric deformations of a free electrically conducting film placed in an
external magnetic field.

We scale horizontal coordinates $(x,y)$ with $L$ and the vertical coordinate
$z$ and the local interface deflection $h(x,y,t)$ with $\langle
h\rangle=\epsilon L$. The horizontal fluid velocity $(u,v)$ is scaled with
some reference velocity $U=O(1)$, the vertical velocity $w$ with $\epsilon
U=O(\epsilon)$, time with $L/U=O(1)$ and the pressure and the disjoining
pressure with $\rho U^{2}$. The magnetic field $B(x,y)$ is non-dimensionalised
using some reference value ${\tilde{B}}$ while the scaling for the electric
potential is $U\tilde{B}L$. The bulk salt and the surface surfactant
concentrations are scaled using arbitrary reference concentrations
$c_{b}^{(0)}$ and $c_{s}^{(0)}$, respectively, so that the conductivity
$\sigma(c_{b})$ is scaled with $Kc_{b}^{(0)}$. The dimensionless bulk
equations 2, 3, 6 and 9 then become

$\displaystyle\partial_{x}u+\partial_{y}v+\partial_{z}w$ $\displaystyle=$
$\displaystyle 0\,,$ (14)
$\displaystyle\partial_{t}u+(u\partial_{x}+v\partial_{y}+w\partial_{z})u$
$\displaystyle=$
$\displaystyle-\partial_{x}[p-\Pi]+\,\textit{Re}^{-1}\left(\partial_{x}^{2}+\partial_{y}^{2}+\epsilon^{-2}\partial_{z}^{2}\right)u$
(15)
$\displaystyle-\,\textit{Ha}^{2}\,\textit{Re}^{-1}c_{b}B(Bu+\partial_{y}\phi)\,,$
$\displaystyle\partial_{t}v+(u\partial_{x}+v\partial_{y}+w\partial_{z})v$
$\displaystyle=$
$\displaystyle-\partial_{y}[p-\Pi]+\,\textit{Re}^{-1}\left(\partial_{x}^{2}+\partial_{y}^{2}+\epsilon^{-2}\partial_{z}^{2}\right)v$
(16)
$\displaystyle-\,\textit{Ha}^{2}\,\textit{Re}^{-1}c_{b}B(Bv-\partial_{x}\phi)\,,$
$\displaystyle\partial_{t}w+(u\partial_{x}+v\partial_{y}+w\partial_{z})w$
$\displaystyle=$
$\displaystyle-\epsilon^{-2}\partial_{z}p+\,\textit{Re}^{-1}\left(\partial_{x}^{2}+\partial_{y}^{2}+\epsilon^{-2}\partial_{z}^{2}\right)w\,,$
(17) $\displaystyle\partial_{x}[c_{b}(Bv-\partial_{x}\phi)]$ $\displaystyle-$
$\displaystyle\partial_{y}[c_{b}(Bu+\partial_{y}\phi)]-\epsilon^{-2}\partial_{z}[c_{b}\partial_{z}\phi]=0\,,$
(18)
$\displaystyle\partial_{t}c_{b}+u\partial_{x}c_{b}+v\partial_{y}c_{b}+w\partial_{z}c_{b}$
$\displaystyle=$
$\displaystyle\,\textit{Sc}^{-1}\,\textit{Re}^{-1}\left(\partial_{x}^{2}+\partial_{y}^{2}+\epsilon^{-2}\partial_{z}^{2}\right)c_{b}\,,$
(19)

where we used the same notations for the dimensionless quantities and
introduced the Reynolds ( Re), Hartmann ( Ha), Peclet ( Pe), and Schmidt ( Sc)
numbers

$\,\textit{Re}=\frac{UL\rho}{\mu}\,,\quad\,\textit{Ha}^{2}=\frac{{\tilde{B}}^{2}L^{2}Kc_{b}^{(0)}}{\mu}\,,\quad\,\textit{Pe}=\frac{LU}{d_{s}}\,,\quad\,\textit{Sc}=\frac{\mu}{\rho
d_{b}}\,.$ (20)

From equation 12, the dimensionless normal and the tangential balances of
stresses at $z=h$ are given by

$\displaystyle
p+(\,\textit{Ca}\,\textit{Re})^{-1}\left(\partial_{x}^{2}h+\partial_{y}^{2}h\right)$
$\displaystyle=$ $\displaystyle
2\,\textit{Re}^{-1}\left(\partial_{z}w-\partial_{x}h\,\partial_{z}u-\partial_{y}h\,\partial_{z}v\right)+O(\epsilon^{2})\,,$
(21) $\displaystyle-\,\textit{Ma}\,\textit{Pe}^{-1}\partial_{x}c_{s}$
$\displaystyle=$
$\displaystyle\epsilon^{-2}\partial_{z}u+\left[-2\partial_{x}h\,\partial_{x}u-\partial_{y}h(\partial_{y}u+\partial_{x}v)+\partial_{x}w\right.$
(22)
$\displaystyle-\left.(\partial_{x}h)^{2}\partial_{z}u-\partial_{x}h\,\partial_{y}h\,\partial_{z}v+2\partial_{x}h\,\partial_{z}w\right]+O(\epsilon^{2})\,,$
$\displaystyle-\,\textit{Ma}\,\textit{Pe}^{-1}\partial_{y}c_{s}$
$\displaystyle=$
$\displaystyle\epsilon^{-2}\partial_{z}v+\left[-2\partial_{y}h\,\partial_{y}u-\partial_{x}h(\partial_{x}v+\partial_{y}u)+\partial_{y}\,w\right.$
(23)
$\displaystyle-\left.(\partial_{y}h)^{2}\partial_{z}v-\partial_{x}h\,\partial_{y}h\,\partial_{z}v+2\partial_{y}h\,\partial_{z}w\right]+O(\epsilon^{2})\,,$

with the capillary ( Ca) and Marangoni ( Ma) numbers defined as

$\,\textit{Ca}=\frac{\mu UL}{\langle
h\rangle\gamma}\,,\quad\,\textit{Ma}=\frac{c_{s}^{(0)}\Gamma_{M}L^{2}}{\langle
h\rangle\mu d_{s}}\,.$ (24)

The scaled boundary condition for the electric current 7, the kinematic
condition 8, the surfactant equation 10 and 11 multiply by
$\sqrt{1+(\partial_{x}h)^{2}+(\partial_{y}h)^{2}}$ are given by

$\displaystyle 0$ $\displaystyle=$ $\displaystyle(\partial_{x}\phi-
Bv)\partial_{x}h+(\partial_{y}\phi-
Bu)\partial_{y}h+\epsilon^{-2}\partial_{z}\phi\,,$ (25)
$\displaystyle\partial_{t}h$ $\displaystyle=$
$\displaystyle-u\partial_{x}h-v\partial_{y}h+w\,,$ (26)
$\displaystyle\partial_{t}c_{s}$ $\displaystyle=$
$\displaystyle-\partial_{x}[uc_{s}]-\partial_{y}[vc_{s}]+\,\textit{Pe}^{-1}\left(\partial_{x}^{2}+\partial_{y}^{2}\right)c_{s}+O\left(\epsilon^{2}\right)\,,$
(27) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\partial_{x}h\,\partial_{x}c_{b}+\partial_{y}h\,\partial_{y}c_{b}-\epsilon^{-2}\partial_{z}c_{b}\,.$
(28)

Crucial for further analysis is to determine the order of magnitude of all
dimensionless parameters appropriate for the physical regime of interest.
Following Chomaz (2001) we assume that for free liquid films the inertial
effects play an essential role implying that $\,\textit{Re}=O(1)$. The leading
contribution to pressure in the fluid is anticipated to come from the Laplace
pressure, which implies that $\,\textit{Ca}=O(1)$. Note that, for example, for
slipper bearing flows and liquid films on a solid substrate, the capillary
number scales as $\,\textit{Ca}=(U\mu/\gamma)\epsilon^{-3}=O(1)$ (Oron et al.,
1997). We assume that the Marangoni effect is weak so that at the leading
order the film surfaces can be considered stress-free (Chomaz, 2001). This can
be achieved by setting $\,\textit{Ma}=O(1)$. An additional assumption must be
made regarding the strength of the magnetic field and the induced electric
current. Here we consider weakly conducting electrolytes in weak magnetic
fields and assume that the Lorentz force is of the same order of magnitude as
the viscous force in the absence of vertical shear, i.e.
${\bm{\nabla}}_{\parallel}^{2}u\sim\,\textit{Ha}^{2}{\bm{\nabla}}_{\parallel}\phi$.
This implies that $\,\textit{Ha}=O(1)$.

All fields in equations (14)–(28) are then expanded into a series in powers of
$\epsilon^{2}$, e.g. $u=u_{0}+\epsilon^{2}u_{1}+\ldots$ with $u_{i}=O(1)$, and
the leading zero-order equations are derived by following the procedure
outlined in Erneux & Davis (1993) and Chomaz (2001). At the zeroth order, the
equations for $h_{0}$, $u_{0}$, $v_{0}$, $w_{0}$, $p_{0}$, and $(c_{s})_{0}$
are identical to those derived in Chomaz (2001) as the Lorentz force enters
the equations only at the next order. Therefore, $u_{0}$, $v_{0}$ and $p_{0}$
are independent of $z$ and $w_{0}=-(\partial_{x}u_{0}+\partial_{y}v_{0})z$.

At the leading order, from equations 18 and 19 we obtain for the electric
potential $\phi_{0}$ and the bulk concentration $(c_{b})_{0}$

$\partial_{z}[(c_{b})_{0}\,\partial_{z}\phi_{0}]=0\,,\quad\partial_{z}^{2}(c_{b})_{0}=0\,.$
(29)

Since $\partial_{z}\phi_{0}=\partial_{z}(c_{b})_{0}=0$ at $z=h_{0}$, from (25)
and (28) we conclude that both $(c_{b})_{0}$ and $\phi_{0}$ are independent of
$z$.

The surfactant concentration $(c_{s})_{0}$ satisfies the two-dimensional
advection-diffusion equation

$\partial_{t}(c_{s})_{0}{\bm{\nabla}}_{\parallel}\cdot[(c_{s})_{0}{\bm{u}}_{0}]=\,\textit{Pe}^{-1}{\bm{\nabla}}_{\parallel}^{2}(c_{s})_{0}\,.$
(30)

where ${\bm{u}}_{0}=(u_{0},v_{0})$ is the leading order horizontal velocity.
The kinematic equation 26 together with
$w_{0}(z=h_{0})=-({\bm{\nabla}}_{\parallel}\cdot{\bm{u}}_{0})h_{0}$ yield the
evolution equation for the local film deformation

$\partial_{t}h_{0}+{\bm{\nabla}}_{\parallel}\cdot[h_{0}u_{0}]=0\,.$ (31)

At the next order, the Navier-Stokes equations for the horizontal flow contain
the Lorentz force terms

$\displaystyle\partial_{t}u_{0}+(u_{0}\partial_{x}+v_{0}\partial_{y})u_{0}$
$\displaystyle=$
$\displaystyle-\partial_{x}[p_{0}-\Pi_{0}]+\,\textit{Re}^{-1}\left(\partial_{x}^{2}+\partial_{y}^{2}\right)u_{0}-\,\textit{Ha}^{2}\,\textit{Re}^{-1}(c_{b})_{0}B_{0}(u_{0}+\partial_{y}\phi_{0})$
(32) $\displaystyle+\,\textit{Re}^{-1}\partial_{z}^{2}u_{2}\,,$
$\displaystyle\partial_{t}v_{0}+(u_{0}\partial_{x}+v_{0}\partial_{y})v_{0}$
$\displaystyle=$
$\displaystyle-\partial_{y}[p_{0}-\Pi_{0}]+\,\textit{Re}^{-1}\left(\partial_{x}^{2}+\partial_{y}^{2}\right)v_{0}-\,\textit{Ha}^{2}\,\textit{Re}^{-1}(c_{b})_{0}B_{0}(v_{0}-\partial_{x}\phi_{0})$
(33) $\displaystyle+\,\textit{Re}^{-1}\partial_{z}^{2}v_{2}\,,$

where the electric potential $\phi_{2}$ satisfies the equation

$\partial_{x}[(c_{b})_{0}(B_{0}v_{0}-\partial_{x}\phi_{0})]-\partial_{y}[(c_{b})_{0}(B_{0}u_{0}+\partial_{y}\phi_{0})]+(c_{b})_{0}\partial_{z}^{2}\phi_{2}=0\,.$
(34)

The boundary condition for $\phi_{2}$ at the film surface $z=h_{0}$ obtained
from equation 25 is

$\partial_{z}\phi_{2}=(B_{0}u_{0}+\partial_{y}\phi_{0})\partial_{y}h_{0}-(B_{0}v_{0}-\partial_{x}\phi_{0})\partial_{x}h_{0}\,.$
(35)

Equation 34 shows that $\phi_{2}$ is a quadratic function of $z$

$\phi_{2}=A(x,y)z^{2}+B(x,y)z+C(x,y)\,,$ (36)

Since for symmetric deformations the potential $\phi$ must be an even function
of $z$, we set $B(x,y)=0$ and determine $A(x,y)$ from the boundary condition
35 to obtain

$\phi_{2}=-\frac{[(B_{0}v_{0}-\partial_{x}\phi_{0})\partial_{x}h_{0}-(B_{0}u_{0}+\partial_{y}\phi_{0})]z^{2}\partial_{y}h_{0}}{2h_{0}}+C(x,y)\,.$
(37)

Finally, substituting 37 into 34 we obtain

$\displaystyle\partial_{x}[(c_{b})_{0}(B_{0}v_{0}-\partial_{x}\phi_{0})]$
$\displaystyle-$
$\displaystyle\partial_{y}[(c_{b})_{0}(B_{0}u_{0}+\partial_{y}\phi_{0})]$ (38)
$\displaystyle+$
$\displaystyle\frac{(c_{b})_{0}[(B_{0}v_{0}-\partial_{x}\phi_{0})\partial_{x}h_{0}-(B_{0}u_{0}+\partial_{y}\phi_{0})\partial_{y}h_{0}]}{h_{0}}=0\,.$

Multiplying equation 38 by $h_{0}$ and introducing the current
${\bm{j}}_{0}=(c_{b})_{0}({\bm{u}}_{0}\times{\bm{B}}_{0}-{\bm{\nabla}}_{\parallel}\phi_{0})$
we rewrite equation 38 in an invariant vector form

${\bm{\nabla}}_{\parallel}\cdot(h_{0}{\bm{j}}_{0})=0\,.$ (39)

Equation 39 represents the continuity equation for the electric current per
unit length of the cross-section of the film.

Next, we eliminate $u_{2}$ and $v_{2}$ from equations 33 by taking into
account the boundary conditions for the tangential and normal components of
the stress tensor at $z=h_{0}$. Because the electromagnetic component of the
viscous stress tensor is neglected here, the result of the elimination
procedure is identical to that of Chomaz (2001). Consequently, we arrive at
the leading-order dynamic equation for ${\bm{u}}_{0}$ including the Lorentz
force

$\displaystyle\partial_{t}{\bm{u}}_{0}$ $\displaystyle+$
$\displaystyle({\bm{u}}_{0}\cdot{\bm{\nabla}}_{\parallel}){\bm{u}}_{0}=\frac{\textrm{d}\Pi_{0}}{\textrm{d}h_{0}}{\bm{\nabla}}_{\parallel}h_{0}+(\,\textit{Ca}\,\textit{Re})^{-1}{\bm{\nabla}}_{\parallel}\nabla_{\parallel}^{2}h_{0}+3\,\textit{Re}^{-1}{\bm{\nabla}}_{\parallel}({\bm{\nabla}}_{\parallel}\cdot{\bm{u}}_{0})$
(40) $\displaystyle+$
$\displaystyle\,\textit{Re}^{-1}\nabla_{\parallel}^{2}{\bm{u}}-\frac{(\,\textit{Re}\,\textit{Pe})^{-1}\,\textit{Ma}}{h_{0}}{\bm{\nabla}}_{\parallel}(c_{s})_{0}+\frac{\,\textit{Re}^{-1}}{h_{0}}{\bm{V}}+\,\textit{Ha}^{2}\,\textit{Re}^{-1}{\bm{j}}_{0}\times{\bm{B}}_{0}\,,$

where we replaced the vector
$(c_{b})_{0}B_{0}[(-\partial_{y}\phi_{0},\partial_{x}\phi_{0})-B_{0}{\bm{u}}_{0}]$
with ${\bm{j}}_{0}\times{\bm{B}}_{0}$ and introduced an additional flow field

${\bm{V}}=2({\bm{\nabla}}_{\parallel}h_{0}\cdot{\bm{\nabla}}_{\parallel}){\bm{u}}_{0}+{\bm{\nabla}}_{\parallel}h_{0}\times({\bm{\nabla}}_{\parallel}\times{\bm{u}}_{0})+2{\bm{\nabla}}_{\parallel}h_{0}({\bm{\nabla}}_{\parallel}\cdot{\bm{u}}_{0})$
(41)

that can be associated with the so-called extensional Trouton viscosity.

To close the system of leading order dynamic equations we derive the equation
for the bulk concentration $(c_{b})_{0}$. At the leading order, from 29 we see
that $(c_{b})_{0}$ is independent of $z$. At the next order, from 19 we obtain

$\partial_{t}(c_{b})_{0}+u_{0}\partial_{x}(c_{b})_{0}+v_{0}\partial_{y}(c_{b})_{0}=\,\textit{Sc}^{-1}\,\textit{Re}^{-1}(\partial_{x}^{2}+\partial_{y}^{2})(c_{b})_{0}+\,\textit{Sc}^{-1}\,\textit{Re}^{-1}\partial_{z}^{2}(c_{b})_{2}\,.$
(42)

The boundary condition at $z=h_{0}$ following from 28 reads

$\partial_{z}(c_{b})_{2}=\partial_{x}h_{0}\partial_{x}(c_{b})_{0}+\partial_{y}h_{0}\partial_{y}(c_{b})_{0}\,.$
(43)

For a symmetric mode and according to 42 the field $(c_{b})_{2}$ is a
quadratic function of $z$: $(c_{b})_{2}=a(x,y,t)z^{2}+c(x,y,t)$. Applying 43
we write

$(c_{b})_{2}=\frac{\partial_{x}h_{0}\,\partial_{x}(c_{b})_{0}+\partial_{y}h_{0}\,\partial_{y}(c_{b})_{0}}{2h_{0}}z^{2}+c(x,y,t)\,.$
(44)

Substituting 44 into 42 we obtain

$\displaystyle\partial_{t}(c_{b})_{0}$ $\displaystyle+$ $\displaystyle
u_{0}\partial_{x}(c_{b})_{0}+v_{0}\partial_{y}(c_{b})_{0}=\,\textit{Sc}^{-1}\,\textit{Re}^{-1}(\partial_{x}^{2}+\partial_{y}^{2})(c_{b})_{0}$
(45) $\displaystyle+$
$\displaystyle\,\textit{Sc}^{-1}\,\textit{Re}^{-1}\frac{\partial_{x}h_{0}\,\partial_{x}(c_{b})_{0}+\partial_{y}h_{0}\,\partial_{y}(c_{b})_{0}}{h_{0}}\,.$

Multiplying 45 by $h_{0}$ and using the kinematic condition 31 we arrive at

$\partial_{t}[h_{0}(c_{b})_{0}]+{\bm{\nabla}}_{\parallel}\cdot[h_{0}{\bm{u}}_{0}(c_{b})_{0}]=\,\textit{Sc}^{-1}\,\textit{Re}^{-1}{\bm{\nabla}}_{\parallel}\cdot[h_{0}{\bm{\nabla}}_{\parallel}(c_{b})_{0}]\,.$
(46)

Equation 46 is identical to the transport equation for the solute
concentration obtained at the leading order of the lubrication approximation
after averaging over the film cross-section that was derived in Jensen &
Grotberg (1993). It generalises the leading order equation for the bulk
concentration derived in Chomaz (2001) to the case of large film deformations.

We summarise our results by writing the complete set of dimensional governing
equations using physical variables and parameters

$\displaystyle\rho(\partial_{t}{\bm{u}}+({\bm{u}}\cdot{\bm{\nabla}}_{\parallel}){\bm{u}})$
$\displaystyle=$ $\displaystyle
g(h){\bm{\nabla}}_{\parallel}h+\gamma{\bm{\nabla}}_{\parallel}{\bm{\nabla}}_{\parallel}^{2}h+\mu{\bm{\nabla}}_{\parallel}^{2}{\bm{u}}+3\mu{\bm{\nabla}}_{\parallel}({\bm{\nabla}}_{\parallel}\cdot{\bm{u}})$
$\displaystyle-\frac{\Gamma_{M}}{h}{\bm{\nabla}}_{\parallel}c_{s}+\frac{\mu{\bm{V}}}{h}+{\bm{j}}\times{\bm{B}}\,,$
$\displaystyle{\bm{V}}$ $\displaystyle=$ $\displaystyle
2({\bm{\nabla}}_{\parallel}h\cdot{\bm{\nabla}}_{\parallel}){\bm{u}}+{\bm{\nabla}}_{\parallel}h\times({\bm{\nabla}}_{\parallel}\times{\bm{u}})+2{\bm{\nabla}}_{\parallel}h({\bm{\nabla}}_{\parallel}\cdot{\bm{u}})\,,$
$\displaystyle{\bm{\nabla}}_{\parallel}\cdot(h{\bm{j}})$ $\displaystyle=$
$\displaystyle 0\,,$ $\displaystyle{\bm{j}}$ $\displaystyle=$ $\displaystyle
Kc_{b}[{\bm{u}}\times{\bm{B}}-{\bm{\nabla}}_{\parallel}\phi]\,,$
$\displaystyle\partial_{t}h$ $\displaystyle=$
$\displaystyle-{\bm{\nabla}}_{\parallel}\cdot(h{\bm{u}})\,,$
$\displaystyle\partial_{t}(hc_{b})$ $\displaystyle=$
$\displaystyle-{\bm{\nabla}}_{\parallel}\cdot(hc_{b}{\bm{u}})+d_{b}{\bm{\nabla}}_{\parallel}\cdot(h{\bm{\nabla}}_{\parallel}c_{b})\,,$
$\displaystyle\partial_{t}c_{s}$ $\displaystyle=$
$\displaystyle-{\bm{\nabla}}_{\parallel}\cdot(c_{s}{\bm{u}})+d_{s}{\bm{\nabla}}_{\parallel}^{2}c_{s}\,,$
(47)

where we introduced function $g(h)=\textrm{d}\Pi(h)/\textrm{d}h$ and dropped
subscript $0$. Note that equations 47 are written in a compact vector form,
which is invariant with respect to the choice of a coordinate system. This is
especially important in applications with non-rectangular geometry as
exemplified in the next section.

##  III Azimuthal flow in an annular free film between two coaxial cylinders

Electromagnetically driven flows of electrolytes in annular bounded by
cylindrical vertical electrodes and a solid bottom have been extensively
studied experimentally and theoretically (Pérez-Barrera et al., 2016; Suslov
et al., 2017; McCloughan & Suslov, 2020, .e.g.). It was found that the steady
azimuthal flow may become unstable giving rise to free-surface vortices
developing close to the outer cylindrical wall. The steady flow field has a
three-dimensional toroidal structure while the deformation of the upper free
surface is negligible.

However, in thin liquid layers the film deformation can no longer be neglected
as has demonstrated in (Wu et al., 1995) using a mechanically driven flow in
an unsupported soap film spanning the gap between two thin coaxial discs. If
the outer disc is fixed and the inner one is rotated, the fluid is set in
motion in an azimuthal direction similar to a Couette cell flow. Centrifugal
forces push the liquid towards the outer disc making the film thinner near the
inner disc. Rather unexpectedly, the flow was found to be laminar and the
onset of turbulence was not observed even at the linear rotation speed of up
to 3 ms-1.

Inspired by these experiments, we consider an electromagnetically driven flow
in an unsupported free film between two conducting coaxial cylindrical
electrodes with radii $R_{1}$ and $R_{2}>R_{1}$ and placed in a vertical
uniform magnetic field ${\bm{B}}=(0,0,B)$ as schematically shown in figure 2.
The potential difference between the inner and outer electrodes is $V$.

Figure 2:  The top (left) and side (right) views of the radial cross-section
of an annular free film of electrically conducting fluid spanning the gap
between two coaxial cylindrical electrodes with radii $R_{1}$ and
$R_{2}>R_{1}$. The film is placed in a vertical uniform magnetic field
${\bm{B}}=(0,0,B)$. Electric current flowing through the film between the two
electrodes generates Lorentz force that drives the flow azimuthally.

In what follows we scale the radial coordinate $r$ with $R_{2}-R_{1}$, the
film thickness $h$ with its average value $\langle h\rangle$ and choose the
flow velocity scale in such a way that $\,\textit{Ca}\,\textit{Re}=1$ in 20,
that is $U=\sqrt{\gamma\langle h\rangle/\left(\rho(R_{2}-R_{1})^{2}\right)}$.
With $B$ used as the magnetic field scaling, the Hartmann number becomes
$\,\textit{Ha}^{2}=B^{2}(R_{2}-R_{1})^{2}Kc_{b}^{(0)}/\mu$, where
$c_{b}^{(0)}$ is the average solute concentration in the bulk. The electric
potential is scaled with $UB(R_{2}-R_{1})$, the pressure with $\rho U^{2}$ and
all other dimensionless parameters Ma, Ha, Ca, Pe and Sc are obtained from 20
by setting $L=R_{2}-R_{1}$. The dimensionless inner and outer radii are given
by $\alpha$ and $1+\alpha$, respectively, where $\alpha=R_{1}/(R_{2}-R_{1})$.
The concentration of a surfactant is scaled using the average value
$c_{s}^{(0)}$. Using the same symbols for non-dimensionless fields, we convert
the invariant form of equations 47 to polar coordinates $(r,\theta)$ and
obtain

$\displaystyle\partial_{t}u_{r}$ $\displaystyle+$ $\displaystyle
u_{r}\partial_{r}u_{r}+\frac{u_{\theta}}{r}\partial_{\theta}u_{r}-\frac{u_{\theta}^{2}}{r}=g(h)\partial_{r}h+\partial_{r}\left(\frac{\partial_{r}h}{r}+\partial_{r}^{2}h+\frac{\partial_{\theta}^{2}h}{r^{2}}\right)$
(48) $\displaystyle+$
$\displaystyle\,\textit{Ca}\left(\frac{\partial_{r}u_{r}}{r}+\partial_{r}^{2}u_{r}+\frac{\partial_{\theta}^{2}u_{r}}{r^{2}}-\frac{u_{r}}{r^{2}}-\frac{2\partial_{\theta}u_{\theta}}{r^{2}}\right)-\,\textit{Ca}\,\textit{Ha}^{2}u_{r}-\,\textit{Ca}\,\textit{Ha}^{2}\frac{\partial_{\theta}\phi}{r}+\frac{\,\textit{Ca}}{h}V_{r}$
$\displaystyle+$ $\displaystyle
3\,\textit{Ca}\partial_{r}\left(\frac{u_{r}}{r}+\partial_{r}u_{r}+\frac{\partial_{\theta}u_{\theta}}{r}\right)-\frac{\,\textit{Ca}\,\textit{Ma}}{\,\textit{Pe}}\frac{\partial_{r}c_{s}}{h}\,,$
$\displaystyle\partial_{t}u_{\theta}$ $\displaystyle+$ $\displaystyle
u_{r}\partial_{r}u_{\theta}+\frac{u_{\theta}}{r}\partial_{\theta}u_{\theta}+\frac{u_{\theta}u_{r}}{r}=\frac{1}{r}g(h)\partial_{\theta}h+\frac{1}{r}\partial_{\theta}\left(\frac{\partial_{r}h}{r}+\partial_{r}^{2}h+\frac{\partial_{\theta}^{2}h}{r^{2}}\right)$
(49) $\displaystyle+$
$\displaystyle\,\textit{Ca}\left(\frac{\partial_{r}u_{\theta}}{r}+\partial_{r}^{2}u_{\theta}+\frac{\partial_{\theta}^{2}u_{\theta}}{r^{2}}-\frac{u_{\theta}}{r^{2}}+\frac{2\partial_{\theta}u_{r}}{r^{2}}\right)+\,\textit{Ca}\,\textit{Ha}^{2}\partial_{r}\phi+\frac{\,\textit{Ca}}{h}V_{\theta}-\,\textit{Ca}\,\textit{Ha}^{2}u_{\theta}$
$\displaystyle+$
$\displaystyle\frac{3\,\textit{Ca}}{r}\partial_{\theta}\left(\frac{u_{r}}{r}+\partial_{r}u_{r}+\frac{\partial_{\theta}u_{\theta}}{r}\right)-\frac{\,\textit{Ca}\,\textit{Ma}}{\,\textit{Pe}}\frac{\partial_{\theta}c_{s}}{hr}\,,$

with

$\displaystyle V_{r}$ $\displaystyle=$ $\displaystyle
2\left(\partial_{r}h\,\partial_{r}u_{r}+\frac{1}{r^{2}}\partial_{\theta}h\,\partial_{\theta}u_{r}-\frac{1}{r^{2}}u_{\theta}\,\partial_{\theta}h\right)+\frac{1}{r^{2}}(\partial_{r}[ru_{\theta}]-\partial_{\theta}u_{r})\,\partial_{\theta}h$
(50)
$\displaystyle+\frac{2}{r}(\partial_{r}[ru_{r}]+\partial_{\theta}u_{\theta})\partial_{r}h\,,$
$\displaystyle V_{\theta}$ $\displaystyle=$ $\displaystyle
2\left(\partial_{r}h\,\partial_{r}u_{\theta}+\frac{1}{r^{2}}\partial_{\theta}h\,\partial_{\theta}u_{\theta}+\frac{1}{r^{2}}u_{r}\,\partial_{\theta}h\right)-\frac{1}{r}(\partial_{r}[ru_{\theta}]-\partial_{\theta}u_{r})\,\partial_{r}h$
(51)
$\displaystyle+\frac{2}{r^{2}}(\partial_{r}[ru_{r}]+\partial_{\theta}u_{\theta})\,\partial_{\theta}h\,.$

The dynamic equations for the salt and surfactant concentrations and the
kinematic condition are given by

$\displaystyle\partial_{t}[hc_{b}]$ $\displaystyle+$
$\displaystyle\frac{1}{r}(\partial_{r}[rhc_{b}u_{r}]+\partial_{\theta}[hc_{b}u_{\theta}])=\frac{\,\textit{Sc}^{-1}\,\textit{Re}^{-1}}{r^{2}}[r\partial_{r}[rh\partial_{r}c_{b}]+\partial_{\theta}[h\partial_{\theta}c_{b}]]\,,$
(52) $\displaystyle\partial_{t}c_{s}$ $\displaystyle+$
$\displaystyle\frac{1}{r}(\partial_{r}[rc_{s}u_{r}]+\partial_{\theta}[c_{s}u_{\theta}])=\frac{\,\textit{Pe}^{-1}}{r^{2}}[r\partial_{r}[r\partial_{r}c_{b}]+\partial_{\theta}^{2}c_{s}]\,,$
(53) $\displaystyle\partial_{t}h$ $\displaystyle+$
$\displaystyle\frac{1}{r}(\partial_{r}[rhu_{r}]+\partial_{\theta}[hu_{\theta}])=0\,.$
(54)

The system is completed by the continuity equation for the electric current

$\partial_{r}[c_{b}rh(u_{\theta}-\partial_{r}\phi)]-\partial_{\theta}\left[\frac{c_{b}h}{r}(ru_{r}+\partial_{\theta}\phi)\right]=0\,.$
(55)

System of equations 49–55 admits a steady solution that corresponds to the
azimuthal flow field $u_{r}=0$, $u_{\theta}=f(r)$ induced by axisymmetric
electric potential $\phi(r)$ in a film with axisymmetric profile $h(r)$
containing uniformly dissolved salt with constant bulk concentration $c_{b}=1$
and covered by a uniformly distributed surfactant with surface concentration
$c_{s}=1$. The functions $f(r)$, $h(r)$ and $\phi(r)$ are found from

$\displaystyle\left(\frac{(rh^{\prime})^{\prime}}{r}\right)^{\prime}+g(h)h^{\prime}+\frac{f^{2}}{r}=0\,,$
(56) $\displaystyle
f^{\prime\prime}+\frac{f^{\prime}}{r}-\frac{f}{r^{2}}-\,\textit{Ha}^{2}f+\,\textit{Ha}^{2}\phi^{\prime}+\frac{h^{\prime}}{h}\left(f^{\prime}-\frac{f}{r}\right)=0\,,$
(57) $\displaystyle(rh(\phi^{\prime}-f))^{\prime}=0\,,$ (58)

where primes denote the radial derivative $\frac{\textrm{d}}{\,\textrm{d}r}$.

Equation 58 is integrated once to yield

$\phi^{\prime}=f-\frac{e}{rh}\,,$ (59)

where constant $e$ is linked to the potential difference
$\phi(1+\alpha)-\phi(\alpha)=\Delta\phi$ between the electrodes via

$\Delta\phi=\int_{\alpha}^{1+\alpha}f(r)\,\textrm{d}r-e\int_{\alpha}^{1+\alpha}\frac{\,\textrm{d}r}{rh}\,.$
(60)

The term $e/(rh)$ in equation 59 represents the radial current density

$j_{r}=f-\phi^{\prime}=\frac{e}{rh}\,.$ (61)

The total steady current through the vertical cylindrical section of the film
at any radial location $\alpha\leq r\leq 1+\alpha$ is independent of the
radius of a cross-section and is given by $4\pi rhj_{r}=4\pi e$.

Eliminating $\phi$ from equation 57 we obtain two coupled equations for $f$
and $h$

$\displaystyle\left(\frac{(rh^{\prime})^{\prime}}{r}\right)^{\prime}+g(h)h^{\prime}+\frac{f^{2}}{r}=0\,,$
(62) $\displaystyle
f^{\prime\prime}+\frac{f^{\prime}}{r}-\frac{f}{r^{2}}-\frac{\,\textit{Ha}^{2}e}{rh}+\frac{h^{\prime}}{h}\left(f^{\prime}-\frac{f}{r}\right)=0$
$\displaystyle\,.$ (63)

Fluid velocity vanishes at the surface of the electrodes leading to
$f(\alpha)=f(1+\alpha)=0$. The boundary condition for the film deformation
must be compatible with the long-wave approximation used here, which only
takes into account relatively small film slopes $h^{\prime}\ll 1$. In what
follows we assume $h^{\prime}(\alpha)=h^{\prime}(1+\alpha)=0$.

We are looking for a solution of the boundary value problems 62 and 63 that
corresponds to a given average film half-thickness

$\frac{2}{1+2\alpha}\int_{\alpha}^{1+\alpha}h(r)r\,\textrm{d}r=1$ (64)

and satisfies additional integral condition 60 for any given value of the
applied voltage $\Delta\phi$.

For reference, we derive the approximate analytic solution that corresponds to
the flow in a flat undeformed film. By setting $h=1$ and neglecting the
centrifugal acceleration $f^{2}/r$ we find from equations 63 and 60

$f=\alpha\,\textit{Ha}^{2}e\left[C\left(\frac{\alpha}{r}-\frac{r}{\alpha}\right)+\frac{r}{2\alpha}\ln\frac{r}{\alpha}\right]\,,$
(65)

where $C=\frac{(1+\alpha)^{2}\ln{(1+\alpha^{-1})}}{2(1+2\alpha)}$,
$e=\frac{\Delta\phi}{\alpha^{2}\,\textit{Ha}^{2}D-\ln(1+\alpha^{-1})}$ and
$D=C\ln\left(1+\alpha^{-1}\right)-\frac{1+2\alpha}{8\alpha^{2}}$.

To find non-trivial solutions of the boundary value problem 62, 63 with the
integral condition 64 we use a numerical continuation package AUTO (Krauskopf
et al., 2014; Doedel et al., 1994). The trivial solution $f=0$ and $h=1$ that
exists for $\Delta\phi=0$ is used as a starting point for numerical
continuation with $\Delta\phi$ being gradually increased.

In the absence of the disjoining pressure, that is for $g(h)=0$, the boundary
value problem 62, 63 has an important scaling property. Namely, it contains
three dimensionless parameters, Ha, $\alpha$ and $\Delta\phi$, but only one of
them, $\Delta\phi$, depends explicitly on the dimensional layer thickness
$2\langle h\rangle$ via $V=\Delta\phi B\sqrt{\gamma\langle h\rangle/\rho}$.
This implies that for any fixed Ha and $\alpha$ there exists a universal
branch of solutions parameterised by $\Delta\phi$. A solution that corresponds
to an arbitrary value of the average film half-thickness $\langle h\rangle$
and an arbitrary applied voltage $V$ is found on the universal branch for
$\Delta\phi=(V/B)\sqrt{\rho/(\gamma\langle h\rangle)}$.

Taking into account the above scaling property we consider a free film with an
arbitrary average half-thickness $\langle h\rangle$ spanning the gap between
cylindrical electrodes with radii $R_{1}=1$ cm and $R_{2}=2$ cm ($\alpha=1$).
As an example we chose fluid properties and magnetic field strength similar to
those used in Pérez-Barrera et al. (2016); Suslov et al. (2017):
$Kc_{b}^{(0)}=\sigma=5$ Ohm m-1, $\mu=0.001$ kg m s-1, $\rho=1000$ kg m-3,
$B=0.05$ T. Additionally, we use $\gamma=0.03$ N m-1 for the surface tension
coefficient. These correspond to $\,\textit{Ha}^{2}\approx 1.3\times 10^{-3}$.
As the solution measure we take the maximum flow velocity $f_{\max}$ and the
minimum film thickness $h_{\min}$ attained at the inner cylinder ($r=\alpha$).
The numerically obtained maximum velocity $f_{\max}$ is shown in figure 3($a$)
by the solid line while the value found from the flat-film approximation 65 is
depicted by the dashed line. The minimum film thickness is shown in the log-
log scale in figure 3($b$) indicating that $h_{\min}$ asymptotically
approaches zero as a power law function $\sim(\Delta\phi)^{-2.4}$. This
implies that the solution exists for any fixed value of $\Delta\phi$ no matter
how large it is. However, the film thickness becomes vanishingly small at the
inner electrode indicating that the film is likely to rapture there. The film
profile and the flow velocity at points $1$ and $2$ are shown in figures
3($c$) and 3($d$), respectively. The dashed line in figure 3($d$) corresponds
to the approximate solution 65.

Figure 3: ($a$) Maximum fluid velocity in a free annular film as a function of
the applied voltage $\Delta\phi$ for $\,\textit{Ha}^{2}=1.3\times 10^{-3}$ and
$\alpha=1$. The dashed and solid lines correspond to the flat film
approximation 65 and the numerical solution, respectively. ($b$) Minimum film
thickness $h_{\min}$ at the inner cylinder as a function of the applied
voltage. The dashed line depicts the power law function
$\sim(\Delta\phi)^{-2.4}$. ($c$) Film thickness $h(r)$ for the values of
$\Delta\phi$ at points $1$ and $2$ labelled in panel ($a$). ($d$) Velocity
$f(r)$ for solutions at points $1$ and $2$ in panel ($a$) (solid lines) and
corresponding to the flow in the flat-film approximation 65 (dashed lines).

##  IV Linear stability of the azimuthal flow

In this section, we study the linear stability of the base azimuthal flow
$(u_{r},u_{\theta})=(0,f(r))$ in a free film with the half-thickness
$h_{0}(r)$ in the absence of the disjoining pressure ($\Pi(h)=0$). It is
anticipated that the least stable perturbations are azimuthally invariant due
to the stabilising effect of the surface tension. Indeed, any perturbation
that varies azimuthally must be periodic in $\theta$ and, consequently, be
proportional to $e^{\textrm{i}n\theta}$, $n=0,1,2\ldots$. Therefore, the
magnitude of the stabilising surface tension terms in equation 49 increases as
$n^{3}$ and is the smallest for $n=0$.

Equations 49–55 are linearised about the base flow by writing

$\displaystyle u_{r}=e^{\lambda t}{\tilde{u}}_{r}(r)\,,\quad
u_{\theta}=f(r)+e^{\lambda t}{\tilde{u}}_{\theta}(r)\,,\quad
h=h_{0}+e^{\lambda t}{\tilde{h}}\,,$ $\displaystyle c_{b}=1+e^{\lambda
t}{\tilde{c}}_{b}\,,\quad c_{s}=1+e^{\lambda
t}{\tilde{c}}_{s}\,,\quad\phi^{\prime}=f-\frac{e}{rh_{0}}+e^{\lambda
t}{\tilde{\phi}}^{\prime}(r)\,,$ (66)

where $\lambda$ is the perturbation growth rate, $h_{0}=h_{0}(r)$ is the
steady film profile and the tilded variables represent small-amplitude
perturbations, substituting these in the equations and neglecting the products
of perturbations. After dropping the tildas we obtain

$\displaystyle\lambda u_{r}$ $\displaystyle=$
$\displaystyle-\frac{2fu_{\theta}}{r}+\left(\frac{(rh^{\prime})^{\prime}}{r}\right)^{\prime}+4\,\textit{Ca}\left(\frac{(ru_{r})^{\prime}}{r}\right)^{\prime}-\,\textit{Ca}\,\textit{Ha}^{2}u_{r}$
(67)
$\displaystyle+\frac{2\,\textit{Ca}h_{0}^{\prime}}{h_{0}}\left(\frac{u_{r}}{r}+2u_{r}^{\prime}\right)-\frac{\,\textit{Ca}\,\textit{Ma}}{\,\textit{Pe}h_{0}}c^{\prime}_{s}\,,$
$\displaystyle\lambda u_{\theta}$ $\displaystyle=$ $\displaystyle-
u_{r}f^{\prime}-\frac{fu_{r}}{r}+\,\textit{Ca}\left(\frac{(ru_{\theta})^{\prime}}{r}\right)^{\prime}-\,\textit{Ca}\,\textit{Ha}^{2}\left(u_{\theta}-\phi^{\prime}\right)$
(68)
$\displaystyle+\frac{\,\textit{Ca}h_{0}^{\prime}}{h_{0}}\left(u_{\theta}^{\prime}-\frac{u_{\theta}}{r}\right)+\,\textit{Ca}\left(\frac{h}{h_{0}}\right)^{\prime}\left(f^{\prime}-\frac{f}{r}\right)\,,$
$\displaystyle\lambda h$ $\displaystyle=$
$\displaystyle-\frac{1}{r}\left(rh_{0}u_{r}\right)^{\prime}\,,$ (69)
$\displaystyle\lambda c_{b}h_{0}+\lambda h$ $\displaystyle=$
$\displaystyle-\frac{1}{r}(rh_{0}u_{r})^{\prime}+\,\textit{Sc}^{-1}\,\textit{Re}^{-1}\frac{(rh_{0}(c_{b})^{\prime})^{\prime}}{r}\,,$
(70) $\displaystyle\lambda c_{s}$ $\displaystyle=$
$\displaystyle-\frac{1}{r}\left(ru_{r}\right)^{\prime}+\,\textit{Pe}^{-1}\frac{(r(c_{s})^{\prime})^{\prime}}{r}\,,$
(71) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\left(rh_{0}(u_{\theta}-\phi^{\prime})\right)^{\prime}+\left(\frac{e\left(h+c_{b}h_{0}\right)}{h_{0}}\right)^{\prime}\,.$
(72)

The flow velocity vanishes at $r=\alpha$ and $r=1+\alpha$ so that
$u_{r}=u_{\theta}=0$ there. This leads to the automatic conservation of the
total volume of the fluid

$\int_{\alpha}^{1+\alpha}2\pi
hr\,\textrm{d}r=-\frac{1}{\lambda}\int_{\alpha}^{1+\alpha}\frac{2\pi}{r}\left(rh_{0}u_{r}\right)^{\prime}r\,\textrm{d}r=0\,.$
(73)

Additionally, we require $h^{\prime}(\alpha)=h^{\prime}(1+\alpha)=0$, which
with the help of 69 and conditions
$h_{0}^{\prime}(\alpha)=h_{0}^{\prime}(1+\alpha)=0$ and
$u_{r}(\alpha)=u_{\theta}(\alpha)=u_{r}(1+\alpha)=u_{\theta}(1+\alpha)=0$
translates into two boundary conditions for the radial velocity $u_{r}$

$u_{r}^{\prime}(\alpha)+u_{r}^{\prime\prime}(\alpha)=\frac{u_{r}^{\prime}(1+\alpha)}{1+\alpha}+u_{r}^{\prime\prime}(1+\alpha)=0\,.$
(74)

Integrating the equation for the perturbation of the electric potential 72
once we obtain

$u_{\theta}-\phi^{\prime}=\frac{c}{rh_{0}}-\frac{eh}{rh_{0}^{2}}-\frac{ec_{b}}{rh_{0}}\,,$
(75)

where $c$ is some constant. Integrating 75 and imposing the condition
$\phi(\alpha)=\phi(1+\alpha)=0$ we find

$\int_{\alpha}^{1+\alpha}u_{\theta}\,\textrm{d}r=c\int_{\alpha}^{1+\alpha}\frac{\,\textrm{d}r}{rh_{0}}-e\int_{\alpha}^{1+\alpha}\frac{h}{rh_{0}^{2}}\,\textrm{d}r-e\int_{\alpha}^{1+\alpha}\frac{c_{b}}{rh_{0}}\,\textrm{d}r\,.$
(76)

Multiplying equations 67 and 68 by $\lambda$, using equations 69 and 75 to
eliminate $h$ and $\phi^{\prime}$ and, subsequently, differentiating 70 and 71
with respect to the radius $r$ we arrive at a nonlinear eigenvalue problem

$\displaystyle\lambda^{2}u_{r}-\frac{2\lambda fu_{\theta}}{r}$
$\displaystyle=$
$\displaystyle-\left(\frac{(r(r^{-1}\left(rh_{0}u_{r}\right)^{\prime})^{\prime})^{\prime}}{r}\right)^{\prime}+4\lambda\,\textit{Ca}\left(\frac{(ru_{r})^{\prime}}{r}\right)^{\prime}-\lambda\,\textit{Ca}\,\textit{Ha}^{2}u_{r}$
(77) $\displaystyle+\frac{2\,\textit{Ca}\lambda
h_{0}^{\prime}}{h_{0}}\left(\frac{u_{r}}{r}+2u_{r}^{\prime}\right)-\frac{\,\textit{Ca}\,\textit{Ma}}{\,\textit{Pe}}\frac{\lambda
c_{s}^{\prime}}{h_{0}}\,,$ $\displaystyle\lambda^{2}u_{\theta}+\lambda
u_{r}f^{\prime}$ $\displaystyle=$ $\displaystyle-\frac{\lambda
fu_{r}}{r}+\lambda\,\textit{Ca}\left(\frac{(ru_{\theta})^{\prime}}{r}\right)^{\prime}-\,\textit{Ca}\,\textit{Ha}^{2}\left(\frac{e(rh_{0}u_{r})^{\prime}}{r^{2}h_{0}^{2}}-\frac{\lambda
ec_{b}}{rh_{0}}+\frac{\lambda c}{rh_{0}}\right)$
$\displaystyle+\frac{\lambda\,\textit{Ca}h_{0}^{\prime}}{h_{0}}\left(u_{\theta}^{\prime}-\frac{u_{\theta}}{r}\right)-\,\textit{Ca}\left(\frac{(rh_{0}u_{r})^{\prime}}{rh_{0}}\right)^{\prime}\left(f^{\prime}-\frac{f}{r}\right)\,,$
$\displaystyle\lambda c_{b}^{\prime}$ $\displaystyle=$
$\displaystyle\,\textit{Sc}^{-1}\,\textit{Re}^{-1}\left(\frac{(rh_{0}c_{b}^{\prime})^{\prime}}{rh_{0}}\right)^{\prime}\,,$
(78) $\displaystyle\lambda c_{s}^{\prime}$ $\displaystyle=$
$\displaystyle\,\textit{Pe}^{-1}\left(\frac{(rc_{s}^{\prime})^{\prime}}{r}\right)^{\prime}-\left(\frac{\left(ru_{r}\right)^{\prime}}{r}\right)^{\prime}$
(79)

that must be solved in conjunction with the integral condition 76. The set of
boundary conditions 74 and
$u_{r}(\alpha)=u_{\theta}(\alpha)=u_{r}(1+\alpha)=u_{\theta}(1+\alpha)=0$ must
be extended with

$\displaystyle
c_{s}^{\prime}(\alpha)=c_{b}^{\prime}(\alpha)=c_{s}^{\prime}(1+\alpha)=c_{b}^{\prime}(1+\alpha)=0\,,$
(80)

which accounts for the chemically passive impenetrable boundaries of the two
electrodes. Note that unlike 62 and 63 the eigenvalue problem 79 contains the
average film thickness $\langle h\rangle$ as a part of the capillary number
Ca.

###  IV.1 Linear stability of a flat film with no applied voltage

In this section, we consider the stability of a flat film in the absence of an
electric current. With $h_{0}=1$, $f=0$ and $e=0$ the eigenvalue problem 77–79
and the integral condition 76 reduce to three decoupled eigenvalue problems:
one for the azimuthal velocity perturbations $u_{\theta}$, one for the radial
velocity $u_{r}$ and surfactant $c_{s}$ perturbations and one for the
perturbation of the solute concentration $c_{b}$:

$\displaystyle\lambda u_{\theta}$ $\displaystyle=$
$\displaystyle\,\textit{Ca}\mathcal{L}\\{{[}\\}u_{\theta}]-\frac{\,\textit{Ca}\,\textit{Ha}^{2}}{r\ln\left(1+\alpha^{-1}\right)}\int_{\alpha}^{1+\alpha}u_{\theta}\,\textrm{d}r\,,$
(81) $\displaystyle\lambda^{2}u_{r}$ $\displaystyle=$
$\displaystyle-\mathcal{L}\\{\mathcal{L}\\{{[}\\}\\}u_{r}]+4\lambda\,\textit{Ca}\mathcal{L}\\{{[}\\}u_{r}]-\lambda\,\textit{Ca}\,\textit{Ha}^{2}u_{r}-\,\textit{Ca}\,\textit{Ma}\,\textit{Pe}^{-1}\lambda
c_{s}^{\prime}\,,$ (82) $\displaystyle\lambda c_{s}^{\prime}$ $\displaystyle=$
$\displaystyle\,\textit{Pe}^{-1}\mathcal{L}\\{{[}\\}c_{s}^{\prime}]-\mathcal{L}\\{{[}\\}u_{r}]\,,$
(83) $\displaystyle\lambda c_{b}^{\prime}$ $\displaystyle=$
$\displaystyle\,\textit{Sc}^{-1}\,\textit{Re}^{-1}\mathcal{L}\\{{[}\\}c_{b}^{\prime}]\,,$
(84)

where $\mathcal{L}\\{{[}\\}u]\equiv\left(r^{-1}(ru)^{\prime}\right)^{\prime}$.

The eigenvalue problem 81–84 can be solved analytically in terms of the
auxiliary eigenvalue problem for the operator $\mathcal{L}$

$\mathcal{L}\\{{[}\\}u]=-\Lambda u\,,$ (85)

which coincides with Bessel’s differential equation. The solution of 85 that
satisfies the Dirichlet boundary conditions $u(\alpha)=u(1+\alpha)=0$ exists
only for real positive $\Lambda$ and is given by

$u(r)=C_{1}J_{1}\left(r\sqrt{\Lambda}\right)+C_{2}Y_{1}\left(r\sqrt{\Lambda}\right)\,,$
(86)

where $C_{1}$ and $C_{2}$ are arbitrary constants and $J_{1}$ and $Y_{1}$ are
the Bessel functions of order one of the first and second kind, respectively.
Applying the boundary conditions we obtain the solvability condition for
constants $C_{1}$ and $C_{2}$ that determines the entire spectrum of discrete
eigenvalues $\Lambda$

$J_{1}\left(\alpha\sqrt{\Lambda}\right)Y_{1}\left((1+\alpha)\sqrt{\Lambda}\right)-J_{1}\left((1+\alpha)\sqrt{\Lambda}\right)Y_{1}\left(\alpha\sqrt{\Lambda}\right)=0\,.$
(87)

It follows from equations 85–87 that the spectrum of eigenvalues
$\lambda_{c_{b}}$ of the bulk solute concentration perturbations 84 is real
and negative

$\lambda_{c_{b}}=-(\,\textit{Sc}\,\textit{Re})^{-1}\Lambda$ (88)

and the corresponding eigenfunction is given by

$c_{b}^{\prime}(r)=C\left[J_{1}\left(r\sqrt{\Lambda}\right)Y_{1}\left(\alpha\sqrt{\Lambda}\right)-Y_{1}\left(r\sqrt{\Lambda}\right)J_{1}\left(\alpha\sqrt{\Lambda}\right)\right]\,,$
(89)

where $C$ is an arbitrary constant.

The eigenvalue problem 82 and 83 for $u_{r}$ and $c_{s}^{\prime}$ can be
solved in a similar way using the eigenfunctions of the operator
$\mathcal{L}$. The solution that satisfies $u_{r}(\alpha)=u_{r}(1+\alpha)=0$
and $c_{s}^{\prime}(\alpha)=c_{s}^{\prime}(1+\alpha)=0$ is given by

$\displaystyle u_{r}(r)$ $\displaystyle=$ $\displaystyle
C_{1}\left[J_{1}\left(r\sqrt{\Lambda}\right)Y_{1}\left(\alpha\sqrt{\Lambda}\right)-Y_{1}\left(r\sqrt{\Lambda}\right)J_{1}\left(\alpha\sqrt{\Lambda}\right)\right]\,,$
$\displaystyle c_{s}^{\prime}(r)$ $\displaystyle=$ $\displaystyle
C_{2}u_{r}\,,$ (90)

where $C_{1}$ and $C_{2}$ are some constants.

Substituting 90 into 82 and 83 we obtain

$\displaystyle\lambda^{2}$ $\displaystyle=$
$\displaystyle-\Lambda^{2}-\,\textit{Ca}(4\Lambda+\,\textit{Ha}^{2})\lambda-
C_{2}\,\textit{Ca}\,\textit{Ma}\,\textit{Pe}^{-1}\lambda\,,$ (91)
$\displaystyle\lambda C_{2}$ $\displaystyle=$
$\displaystyle\Lambda-\,\textit{Pe}^{-1}\Lambda C_{2}$ (92)

and then eliminating $C_{2}$ from 91 we arrive at a cubic equation for
$\lambda$

$(\,\textit{Pe}\lambda+\Lambda)\left(\lambda^{2}+\,\textit{Ca}\left(4\Lambda+\,\textit{Ha}^{2}\right)\lambda+\Lambda^{2}\right)+b\Lambda\lambda\,\textit{Pe}=0\,,$
(93)

where $b=\,\textit{Ca}\,\textit{Ma}\,\textit{Pe}^{-1}>0$ characterises the
strength of the Marangoni flow.

For any admissible value of $\Lambda$ from 87 one needs to solve 93 to find
the eigenvalue $\lambda$. Then the corresponding eigenfunctions 90 only
contain one arbitrary scaling factor $C_{1}$ with
$C_{2}=\Lambda/(\lambda+\,\textit{Pe}^{-1}\Lambda)$. Here we consider two
physically distinct situations: the diffusion-dominated regime, when the
surface diffusivity is large, i.e. $d_{s}\to\infty$, and the advection-
dominated regime, when $d_{s}\to 0$. Since the Hartmann number Ha, the
capillary number Ca and parameter
$b=\,\textit{Ca}\,\textit{Ma}\,\textit{Pe}^{-1}$ in equation 93 do not depend
on $d_{s}$, the diffusion-dominated regime corresponds to $\,\textit{Pe}\to 0$
with Ha, Ca and $b$ remaining finite. In this case equation 93 has three
distinct solutions

$\displaystyle\lambda_{1,2}$ $\displaystyle=$
$\displaystyle-\frac{\,\textit{Ca}\left(4\Lambda+\,\textit{Ha}^{2}\right)+b\,\textit{Pe}}{2}\pm\frac{\Omega}{2}\left(1+b\,\textit{Pe}\frac{\,\textit{Ca}(4\Lambda+\,\textit{Ha}^{2})}{2\Omega^{2}}\right)+O(\,\textit{Pe}^{2})\,,$
(94) $\displaystyle\lambda_{3}$ $\displaystyle=$
$\displaystyle-\,\textit{Pe}^{-1}\Lambda+b\,\textit{Pe}+O(\,\textit{Pe}^{2})\,,$
(95)

where
$\Omega=\sqrt{\,\textit{Ca}^{2}(4\Lambda+\,\textit{Ha}^{2})^{2}-4\Lambda^{2}}$.

Note that $\lambda_{3}$ is real and negative and its magnitude is always much
larger than $|\lambda_{1,2}|$. Since $\Lambda$ is real and positive, we
conclude that the real parts of $\lambda_{1,2}$ are always negative. Moreover,
$\lambda_{1,2}$ become complex if
$\,\textit{Ca}\left(4\Lambda+\,\textit{Ha}^{2}\right)<2\Lambda$ implying that
the radial perturbation mode undergoes the transition from a monotonic to
oscillatory decay. In physical variables, the condition for the oscillatory
decay of the radial mode is

$\sqrt{\rho\gamma\langle
h\rangle}>2\mu+\frac{B^{2}(R_{2}-R_{1})^{2}\sigma}{2\Lambda}\,,$ (96)

where the electric conductivity is $\sigma=Kc_{b}^{(0)}$ and $\Lambda$ depends
on the ratio of the radii $R_{2}/R_{1}$ via parameter $\alpha$. It is
instructive to compare condition 96 with the linear stability of an unbounded
flat free layer in the absence of the magnetic field. Neglecting the
disjoining pressure, the growth rate $\omega(k)$ of the least stable mode with
the wave vector $k$ in an unbounded horizontal free film can be obtained from
equation (30) in Erneux & Davis (1993):

$\frac{\mu}{\rho}\omega(k)=-2k^{2}+k^{2}\sqrt{4-\frac{\rho\langle
h\rangle\gamma}{\mu^{2}}}\,.$ (97)

It follows from 97 that the critical thickness $2\langle h\rangle_{c}$ of the
layer above which the relaxation dynamics is oscillatory is given by

$2\langle h\rangle_{c}=\frac{8\mu^{2}}{\rho\gamma}\,.$ (98)

This coincides with our result 96 for $B=0$.

It is seen from equations 94 that in the absence of the Marangoni flow, that
is when $b=0$, the eigenvalue of the radial velocity perturbation with the
largest real part is given by

$\lambda_{1}=-\frac{\,\textit{Ca}\left(4\Lambda+\,\textit{Ha}^{2}\right)}{2}+\frac{\Omega}{2}\,.$
(99)

The perturbation of the surfactant is decoupled from that of the radial flow
and has a negative real eigenvalue $\lambda_{3}=-\Lambda\,\textit{Pe}^{-1}$.
In this regime, any perturbation of the surfactant distribution relaxes on the
time scale $|\lambda_{3}|^{-1}$, which is much shorter than the characteristic
decay time $-\Re({\lambda_{1,2}}^{-1})$ of fluid motion. In the presence of
the Marangoni flow when $b\neq 0$, the perturbations of the radial velocity
and surfactant concentrations are coupled and their dynamics is characterised
by the leading eigenvalues $\lambda_{1,2}$. The Marangoni flow leads to a
further stabilisation of the leading perturbation mode as follows from
equation 94. Indeed, by analysing the real parts of the leading eigenvalues we
conclude that

$\Re(\lambda_{1,2})|_{b>0}<\Re(\lambda_{1,2})|_{b=0}<0$

regardless of the sign of $\Omega^{2}$. As a consequence, the characteristic
decay time of the flow perturbation decreases in the presence of Marangoni
effect.

As follows from equation 83, in the advection-dominated regime
$\,\textit{Pe}\to\infty$ the dynamics of the surfactant perturbation is
governed by the radial flow $u_{r}$. Indeed in this limit $\lambda
c_{s}^{\prime}=-\mathcal{L}\\{{[}\\}u_{r}]$, which shows that the surfactant
plays the role of an active scalar field advected by the flow while its
gradient influences the stability of the flow. By letting
$\,\textit{Pe}\to\infty$ in 92 and then substituting $C_{2}=\Lambda/\lambda$
in equation 91 we obtain

$\displaystyle\lambda_{\pm}$ $\displaystyle\approx$
$\displaystyle\frac{1}{2}\left[-\,\textit{Ca}\left(4\Lambda+\,\textit{Ha}^{2}\right)\pm\sqrt{\,\textit{Ca}^{2}\left(4\Lambda+\,\textit{Ha}^{2}\right)^{2}-4\Lambda(\Lambda+b)}\right]\,,$
(100)

From 100 we see that the Marangoni flow has no effect on the stability of the
base flow if the expression under the radical is negative. However, if

$\,\textit{Ca}^{2}\left(4\Lambda+\,\textit{Ha}^{2}\right)^{2}-4\Lambda(\Lambda+b)>0\,,$

the leading eigenvalue $\lambda_{+}$ becomes real and negative, and the
presence of surfactant has a stabilising effect since
$|\lambda_{+}|_{b>0}>|\lambda_{+}|_{b=0}$. These analytical results extend an
earlier study on the linear stability of planar soap films (Miksis & Ida,
1998) by including the Lorentz force effects. The observation of the
stabilising role of the Marangoni flow in the limit of diffusion-dominated and
advection-dominated regimes is in agreement with Miksis & Ida (1998), where it
was also found that for long-wavelength perturbations the presence of the
Marangoni flow decreases the growth rate of the dominant perturbation mode.

To illustrate the structure of the oscillatory radial mode we consider the
diffusion-dominated regime in a film with the average thickness $2\langle
h\rangle=20\,\mu$m and take all other parameters as in figure 3. For
$\alpha=1$, the leading eigenvalue of the operator $\mathcal{L}$ is
$\Lambda\approx 10.218$, which corresponds to the complex leading eigenvalue
of the radial mode $\lambda_{r}\approx-1.179\pm 10.149\textrm{i}$. The
corresponding eigenfunction $u_{r}$ is shown in figure 4($a$).

Figure 4: ($a$) Leading eigenfunction
$u_{r}=Y_{1}\left(\alpha\sqrt{\Lambda}\right)J_{1}\left(r\sqrt{\Lambda}\right)-J_{1}\left(\alpha\sqrt{\Lambda}\right)Y_{1}\left(r\sqrt{\Lambda}\right)$
of the radial perturbation mode for the same parameters as in figure 3. ($b$)
Streamlines of the radially perturbed flow field in a vertical crosssection of
the film.

Figure 4($b$) depicts the instantaneous streamlines of the corresponding flow
field in a vertical cross-section of the film $(r,z)$. The streamlines are
obtained by recalling that the vertical flow velocity $w$ is given by
$w(r,z)=-r^{-1}(ru_{r})^{\prime}z$ so that the kinematic equation 69 can be
written in the form $\partial_{t}h=-w(r,1)=-r^{-1}(ru_{r})^{\prime}$. The film
interface oscillates about $h_{0}=1$ with a decreasing amplitude while the
fluid flows from the inner to the outer cylinder and back.

The eigenvalue problem (81) for the azimuthal velocity perturbation
$u_{\theta}$ can only be found analytically for a weak magnetic field at
$\,\textit{Ha}\to 0$. When $\,\textit{Ha}=0$, the spectrum of the azimuthal
velocity perturbations is real and is given by

$\lambda_{u_{\theta}}=-\Lambda\,\textit{Ca}\,.$ (101)

Comparing 101 with 94 we observe that in the absence of the Marangoni flow
($b=0$) and when $\,\textit{Ha}=0$ the absolute value of the eigenvalue
$|\lambda_{u_{\theta}}|$ of the monotonically stable azimuthal velocity
perturbation is exactly half of the real part of the eigenvalue
$\lambda_{u_{r}}$ of the radial velocity mode.

###  IV.2 Linear stability of a deformed film in the presence of electric
current

In this section, we study the linear stability of the azimuthal flow in a
deformed film in the diffusion-dominated regime, when perturbations of the
surfactant and solute fields relax instantaneously. The azimuthal and radial
velocity perturbation modes discussed in the previous section become coupled
in the presence of an electric current. This implies that for any
infinitesimal applied voltage $\Delta\phi$ the spectrum of the generalised
eigenvalue problem 79 is discrete and contains the eigenvalues that originate
from each of the possible radial and azimuthal velocity modes existing in the
absence of current.

To visualise how the stability of the azimuthal flow changes with the applied
voltage we choose $\langle h\rangle=10\,\mu$m and keep all other parameters as
in figure 3. The numerical continuation method is then employed to track each
mode and the corresponding eigenvalue using $\Delta\phi$ as the continuation
parameter. To this end, the generalised eigenvalue problem 79 is solved
simultaneously with equations 62 and63. To quantify the relative coupling
strength between the perturbations of the radial and the azimuthal velocity we
introduce the parameter

$\chi=\arctan{\left(\frac{\int_{\alpha}^{1+\alpha}|u_{r}|^{2}\,\textrm{d}r}{\int_{\alpha}^{1+\alpha}|u_{\theta}|^{2}\,\textrm{d}r}\right)}\,.$
(102)

Thus, $\chi=\pi/2$ corresponds to the pure radial and $\chi=0$ to pure
azimuthal modes, respectively.

Figure 5: ($a$) Real part of the first three leading eigenvalues as a function
of the applied voltage $\Delta\phi$ for the same parameters as in figure 3.
The solid (dashed) lines correspond to complex (purely real) eigenvalues,
respectively. Labels bp and ns mark the loci of the branching point of the
azimuthal velocity mode and the point of neutral stability of the radial
velocity mode, respectively; ($b$) the relative strength of the radial to
azimuthal velocity mode $\chi$; ($c$) the neutrally stable film profile at
point ns in panel ($a$); ($d$) the azimuthal velocity $f_{0}$ in a neutrally
stable film; the magnitudes of the neutrally stable radial ($e$) and azimuthal
($f$) velocity perturbations.

The real parts of the first three leading eigenvalues labelled by $r_{1}$,
$\theta_{1}$ and $\theta_{2}$ are shown in figure 5($a$). At $\Delta\phi=0$
the leading mode is the azimuthal velocity mode $\theta_{1}$ with the real
eigenvalue $\lambda_{\theta_{1}}\approx-0.589$. The second least stable mode
$r_{1}$ corresponds to the radial velocity. It has the complex eigenvalues
$\lambda_{r_{1}}=-1.179\pm 10.149\textrm{i}$. The third mode $\theta_{2}$ is
again the azimuthal velocity mode with the real eigenvalue
$\lambda_{\theta_{2}}\approx-2.298$. As $\Delta\phi$ increases, the two
leading azimuthal velocity modes remain monotonically stable until the
branching point bp is reached, where the two modes form a complex conjugate
pair with a negative real part collide. A further increase of $\Delta\phi$
leads to the destabilisation of the leading radial mode $r_{1}$, which becomes
neutrally stable at the point ns, where $\Im\\{\lambda_{r_{1}}\\}\not=0$. The
corresponding film thickness profile $h_{0}(r)$ and the base azimuthal flow
$f_{0}(r)$ are shown in figures 5($c,d$), respectively.

Oscillatory neutrally stable perturbation is a mixture of radial and azimuthal
flows as quantified by $\chi$ in figure 5($b$). Since the perturbation fields
$u_{r}$ and $u_{\theta}$ are both complex, the real physical radial and
azimuthal perturbations of the base flow are given by the real and imaginary
parts of $|u_{r}|e^{\textrm{i}(\Im\\{\lambda_{r_{1}}\\}t+\Psi_{r}(r)})$ and
$|u_{\theta}|e^{\textrm{i}(\Im\\{\lambda_{r_{1}}\\}t+\Psi_{\theta}(r))}$ with
spatially varying phase shifts
$\Psi_{r}=\arctan(\Im\\{u_{r}\\}/\Re\\{u_{r}\\})$ and
$\Psi_{\theta}=\arctan(\Im\\{u_{\theta}\\}/\Re\\{u_{\theta}\\})$. The
spatially varying amplitudes of the neutrally stable radial and azimuthal
velocity perturbations $|u_{r}|$ and $|u_{\theta}|$ are shown in figure
5($e,f$).

To gain a deeper understanding of how the stability of the base flow changes
with other parameters we plot neutral stability curves in figure 6($a$) for
three selected values of $\,\textit{Ha}^{2}$ in the plane
$(\Delta\phi_{c},\,\textit{Ca})$, where $\Delta\phi_{c}$ is the critical value
of the applied voltage above which the base flow is linearly unstable. It is
noteworthy that for each value of $\,\textit{Ha}^{2}$ there exists a finite
limiting value of $\Delta\phi_{c}$ as $\,\textit{Ca}\to 0$. Such a limiting
behaviour of the neutral stability curves confirms the dynamic nature of the
flow instability. Indeed, with the scaling used here the capillary number is
given by $\,\textit{Ca}=\mu(\gamma\langle h\rangle\rho)^{-1/2}$ and the
physical value of the applied voltage is

$V_{c}=UB(R_{2}-R_{1})\Delta\phi_{c}=B\sqrt{\frac{\gamma\langle
h\rangle}{\rho}}\Delta\phi_{c}\,.$

For the fixed $\langle h\rangle$, viscosity $\mu$ and density $\rho$ the limit
of $\,\textit{Ca}\rightarrow 0$ corresponds to a large surface tension
$\gamma\to\infty$. The existence of a finite value of $\Delta\phi_{c}$ implies
that regardless of the strength of stabilising surface tension forces the
instability can always be induced if the applied voltage exceeds $V_{c}$.

We set $\,\textit{Ca}=10^{-5}$ and trace the locus of the neutral stability in
the $(\Delta\phi_{c},\,\textit{Ha}^{2})$ plane in figure 6($b$). The dashed
line corresponds to $\,\textit{Ha}^{2}=124\Delta\phi_{c}^{-1}$ and is obtained
by fitting the lower part of the curve (small $\,\textit{Ha}^{2}$) with a
power law function. This simple power law relationship between
$\,\textit{Ha}^{2}$ and $\Delta\phi_{c}$ represents a universal stability
threshold in the limit of small capillary and Hartmann numbers. Reintroducing
dimensional variables we obtain the critical value of the applied voltage
$V_{c}$

$V_{c}=\frac{124\mu}{B\sigma(R_{2}-R_{1})^{2}}\sqrt{\frac{\gamma\langle
h\rangle}{\rho}}\,.$ (103)

The asymptotic result 103 is valid for the selected radii ratio
$R_{2}/R_{1}=2$ and $\left(\,\textit{Ca},\,\textit{Ha}^{2}\right)\to(0,0)$. At
the critical value of $V_{c}$ the film is strongly deformed and has a shape
similar to the solution (2) in figure 3($c$). The minimum value of the film
thickness is attained at the inner cylinder and is approximately 6% of
$\langle h\rangle$.

Figure 6: ($a$) Neutral stability curves in the
$(\Delta\phi_{c},\,\textit{Ca})$ plane for $\,\textit{Ha}_{1}^{2}=0.0038$,
$\,\textit{Ha}_{2}^{2}=0.0013$ and $\,\textit{Ha}_{3}^{2}=0.00065$. ($b$)
Neutral stability curve in the $\left(\Delta\phi_{c},\,\textit{Ha}^{2}\right)$
plane for $\,\textit{Ca}=10^{-5}$. The dashed line corresponds to the power
law function $\,\textit{Ha}^{2}=124\Delta\phi_{c}^{-1}$.

##  V Conclusions

We derived a set of leading order equations in the lubrication approximation
that describe the electromagnetically driven flow of electrolyte solutions
and/or non-magnetic liquid metals in a thin horizontal free layer with
deformable surfaces placed in a uniform magnetic field normal to the layer.
The equations are written in a generic coordinate-invariant form and can be
used to study electromagnetically driven flows in an arbitrary geometry
dictated by the shape and size of the deployed electrodes. The equations
account for the presence of surfactants and chemical species dissolved in the
bulk of the fluid. Inspired by recent studies of the electromagnetically
driven flows in supported shallow annular layers between two coaxial cylinders
Figueroa et al. (2009); Pérez-Barrera et al. (2016); Suslov et al. (2017);
McCloughan & Suslov (2020), we choose a similar geometry and apply the derived
model to investigate the flow and its stability in annular free film spanning
the gap between two coaxial cylindrical electrodes.

Similar to mechanically driven flows in soap films (Wu et al., 1995), a steady
azimuthal electromagnetically driven flow can only exist in a deformed layer,
where the radial component of the gradient of the Laplace pressure balances
the centrifugal force. This is in contrast to supported thicker layers, where
the free-surface deformation is typically negligible compared to the thickness
of the layer while the flow speed does not exceed several centimetres per
second. The other important feature distinguishing the flows in free layers
and films from those in supported layers is that in the former the Laplace
pressure dominates while in the latter the hydrostatic pressure gradient
defines the fluid trajectory in the plane of the layer.

If the intermolecular forces are neglected, the azimuthally invariant steady-
state flow can be found for an arbitrarily large electric current flowing
through the film. The minimal film thickness is achieved at the inner
electrode. It decreases as a power law function of the applied voltage
remaining non-zero so that the point of a true film rupture is never reached.

We determined that the azimuthal flow in approximately flat free films with a
small velocity is linearly unconditionally stable. For relatively weak
magnetic fields of the order of $10^{-2}$ T and fluid parameters corresponding
to a weak electrolyte solution used in the experiments of Figueroa et al.
(2009); Pérez-Barrera et al. (2016), we found that the azimuthal velocity
perturbations decay monotonically while the decay of radial velocity
perturbations of the base flow in layers with the average thickness in the
micrometre range is oscillatory. However, as the film thickness is decreased
below a certain critical value given by condition 96 ($\lesssim 100$ nm for a
film existing between coaxial cylinders with the inner and outer radii of 1
and 2 cm, respectively), the relaxation dynamics is found to be dominated by
viscous damping with a monotonic decay of a radial flow and the associated
surface deformation. This result contrasts the observations of instabilities
in supported thicker annular layers, where the primary flow becomes unstable
with respect to three-dimensional perturbations of the velocity without a
noticeable variation of the layer depth. The fluctuations of the solute
concentration in the bulk has no effect on the growth rate of the leading
modes. In the presence of a surfactant, the Marangoni effect leads to further
stabilisation of the base flow, which is in agreement with earlier studies of
planar soap films in the absence of magnetic fields (Miksis & Ida, 1998).

By following the branch of steady-state solutions into the regime of large
applied voltage and, consequently, large electric currents, we find that the
steady azimuthal flow eventually becomes unstable with respect to a mixture of
oscillatory azimuthal and radial velocity perturbations at a certain critical
value of the applied voltage, at which the deformation amplitude of the layer
is of the order of the average film thickness. This suggests that
electromagnetically driven flows in free films may not be ideal candidates for
studying two-dimensional turbulence since the primary instability of the base
flow sets in only when the film is already strongly deformed. We note that a
similar conclusion was made earlier in Rivera & Wu (2000), where it was
hypothesised that strong damping caused by the interactions of the flow within
a film with a surrounding gas layer that are enhanced by the surface
deformation may be the cause for the energy leakage responsible for stronger
than expected decays of the velocity correlation function.

Finally, we briefly mention the action of intermolecular forces characterised
by the disjoining pressure. They play a major role in the instability and
subsequent stabilisation of free soap films that are thinner than $\sim 100$
nm. It is well known Israelachvili (2011); Overbeek (1960) that long-range van
der Waals forces destabilise free layers sandwiched between dielectric media
with identical properties. These forces alongside the gravity-induced drainage
of the fluid constitute the primary source of film instability. As the film
thickness decreases below$10-50$ nm, an electric double layer is typically
formed consisting of the monolayers of soap ions adsorbed at each interface.
Strong electric double-layer forces lead to the repulsion between the film
surfaces and the formation of highly stable black soap films with a thickness
under $50$ nm. Therefore, it is important to study the role of the disjoining
pressure on the steady azimuthal flow and its stability in the strong-current
regimes when the deformation of the layer is significant. Intermolecular
forces, black soap films and film rupture will become the topic of our future
investigations.

## References

  * Bau et al. (2001) Bau, H. H., Zhong, J. & Yi, M. 2001 A minute magneto hydro dynamic (MHD) mixer. Sens. Actuators, B 79 (2), 207–215.
  * Bau et al. (2003) Bau, H. H., Zhu, J., Qian, S. & Xiang, Y. 2003 A magneto-hydrodynamically controlled fluidic network. Sens. Actuators, B 88 (2), 205–216.
  * Cardoso et al. (1994) Cardoso, O., Marteau, D. & Tabeling, P. 1994 Quantitative experimental study of the free decay of quasi-two-dimensional turbulence. Phys. Rev. E 49, 454–461.
  * Chomaz (2001) Chomaz, J.-M. 2001 The dynamics of a viscous soap film with soluble surfactant. J. Fluid Mech. 442, 387–409.
  * Chomaz & Cathalau (1990) Chomaz, J. M. & Cathalau, B. 1990 Soap films as two-dimensional classical fluids. Phys. Rev. A 41, 2243–2245.
  * Couder (1981) Couder, Y. 1981 The observation of a shear flow instability in a rotating system with a soap membrane. J. Phys. Let. 42 (19), 429–431.
  * Couder (1984) Couder, Y. 1984 Two-dimensional grid turbulence in a thin liquid film. J. Phys. Let. 45, 353–360.
  * Couder et al. (1989) Couder, Y., Chomaz, J.-M. & Rabaud, M. 1989 On the hydrodynamics of soap films. Physica D 37 (1), 384–405.
  * Cruz Gómez (2016) Cruz Gómez, R. 2016 Response of a soap film to a continuous electromagnetic forcing. J. Visualization 20, 87–95.
  * Cuevas & Figueroa (2023) Cuevas, S. & Figueroa, A. 2023 Personal communication.
  * Doedel et al. (1994) Doedel, E. J., Wang, X. & Fairgrieve, T. 1994 AUTO94: Software for continuation and bifurcation problems in ordinary differential equations. Applied Mathematics Report. California Institute of Technology, Pasadena, U.S.A.
  * Dunne et al. (2020) Dunne, P., Adachi, T., Dev, A. A., Sorrenti, A., Giacchetti, L., Bonnin, A., Bourdon, C., Mangin, P. H., Coey, J. M. D., Doudin, B. & Hermans, T. M. 2020 Liquid flow and control without solid walls. Nature 581, 58–62.
  * Erneux & Davis (1993) Erneux, T. & Davis, S. H. 1993 Nonlinear rupture of free films. Phys. Fluids A 5 (5), 1117–1122.
  * Fiflis et al. (2016) Fiflis, P., Christenson, M., Szott, M., Kalathiparambil, K. & Ruzic, D. N. 2016 Free surface stability of liquid metal plasma facing components. Nucl. Fusion 56 (10), 106020.
  * Figueroa et al. (2009) Figueroa, A., Demiaux, F., Cuevas, S. & Ramos, E. 2009 Electrically driven vortices in a weak dipolar magnetic field in a shallow electrolytic layer. J. Fluid Mech. 641, 245–261.
  * Gao & Morley (2002) Gao, D. & Morley, N. B. 2002 Equilibrium and initial linear stability analysis of liquid metal falling film flows in a varying spanwise magnetic field. Magnetohydrodynamics 38 (4), 359–375.
  * Gao et al. (2002) Gao, D., Morley, N. B. & Dhir, V. 2002 Numerical study of liquid metal film flows in a varying spanwise magnetic field. Fusion Eng. Des. 63, 369–374.
  * Gharib & Derango (1989) Gharib, M. & Derango, P. 1989 A liquid film (soap film) tunnel to study two-dimensional laminar and turbulent shear flows. Physica D 37 (1), 406–416.
  * Giannakis et al. (2009a) Giannakis, D., Fischer, P. F. & Rosner, R. 2009a A spectral Galerkin method for the coupled Orr–Sommerfeld and induction equations for free-surface MHD. J. Comp. Phys. 228 (4), 1188–1233.
  * Giannakis et al. (2009b) Giannakis, D., Rosner, R. & Fischer, P. F. 2009b Instabilities in free-surface Hartmann flow at low magnetic Prandtl numbers. J. Fluid Mech. 636, 217–277.
  * Ida & Miksis (1998) Ida, M. P. & Miksis, M. J. 1998 The dynamics of thin films I: General theory. SIAM J. Appl. Math. 58 (2), 456–473.
  * Israelachvili (2011) Israelachvili, J. N. 2011 Intermolecular and Surface Forces, 3rd edn. Elsevier Science.
  * Jensen & Grotberg (1993) Jensen, O. E. & Grotberg, J. B. 1993 The spreading of heat or soluble surfactant along a thin liquid film. Phys. Fluids A 5 (1), 58–68.
  * Kellay & Goldburg (2002) Kellay, H. & Goldburg, W. I. 2002 Two-dimensional turbulence: a review of some recent experiments. Rep. Prog. Phys. 65 (5), 845–894.
  * Krauskopf et al. (2014) Krauskopf, B., Osinga, H. M. & Galan-Vioque, J. 2014 Numerical Continuation Methods for Dynamical Systems: Path following and boundary value problems. Springer Netherlands.
  * Lunz & Howell (2019) Lunz, D. & Howell, P. D. 2019 Flow of a thin liquid-metal film in a toroidal magnetic field. J. Fluid Mech. 867, 835–876.
  * Marteau et al. (1995) Marteau, D., Cardoso, O. & Tabeling, P. 1995 Equilibrium states of two-dimensional turbulence: An experimental study. Phys. Rev. E 51, 5124–5127.
  * McCloughan & Suslov (2020) McCloughan, J. & Suslov, S. A. 2020 Linear stability and saddle-node bifurcation of electromagnetically driven electrolyte flow in an annular layer. J. Fluid Mech. 887, A23.
  * Messadek & Moreau (2002) Messadek, K. & Moreau, R. 2002 An experimental investigation of MHD quasi-two-dimensional turbulent shear flows. J. Fluid Mech. 456, 137–159.
  * Miksis & Ida (1998) Miksis, M. J. & Ida, M. P. 1998 The dynamics of thin films II: Applications. SIAM J. Appl. Math. 58 (2), 474–500.
  * Miloshevsky & Hassanein (2010) Miloshevsky, G. V. & Hassanein, A. 2010 Modelling of Kelvin-Helmholtz instability and splashing of melt layers from plasma-facing components in tokamaks under plasma impact. Nucl. Fusion 50 (11), 115005.
  * Morley & Abdou (1995) Morley, N. B. & Abdou, M. A. 1995 Modeling of fully-developed, liquid metal, thin film flows for fusion divertor applications. Fusion Eng. Des. 30 (4), 339–356.
  * Morley & Abdou (1997) Morley, N. B. & Abdou, M. A. 1997 Study of fully developed, liquid-metal, open-channel flow in a nearly coplanar magnetic field. Fusion Technol. 31 (2), 135–153.
  * Morley & Roberts (1996) Morley, N. B. & Roberts, P. H. 1996 Solutions of uniform, open‐channel, liquid metal flow in a strong, oblique magnetic field. Phys. Fluids 8 (4), 923–935.
  * Morley et al. (2002) Morley, N. B., Smolentsev, S. & Gao, D. 2002 Modeling infinite/axisymmetric liquid metal magnetohydrodynamic free surface flows. Fusion Eng. Des. 63–64, 343–351.
  * Müller & Bühler (2013) Müller, U. & Bühler, L. 2013 Magnetofluiddynamics in Channels and Containers. Springer Berlin, Heidelberg.
  * Oron et al. (1997) Oron, A., Davis, S. H. & Bankoff, S. G. 1997 Long-scale evolution of thin liquid films. Rev. Mod. Phys. 69, 931–980.
  * Overbeek (1960) Overbeek, J. Th. G. 1960 Black soap films. J. Phys. Chem. 64 (9), 1178–1183.
  * Pérez-Barrera et al. (2016) Pérez-Barrera, J., Ortiz, A. & Cuevas, S. 2016 Analysis of an annular MHD stirrer for microfluidic applications. In Recent advances in fluid dynamcis with environmental applications (ed. J. Klapp, L. D. G. Sigalotti, A. Medina Ovando, A. López Villa & G. Ruíz Chavarría), pp. 275–288. Springer.
  * Prévost & Gallez (1986) Prévost, M. & Gallez, D. 1986 Nonlinear rupture of thin free liquid films. J. Chem. Phys. 84 (7), 4043–4048.
  * Qian & Bau (2005) Qian, S. & Bau, H. H. 2005 Magneto-hydrodynamic stirrer for stationary and moving fluids. Sens. Actuators, B 106, 859–870.
  * Rivera & Wu (2000) Rivera, M. & Wu, X. L. 2000 External dissipation in driven two-dimensional turbulence. Phys. Rev. Let. 85, 976–979.
  * Shang et al. (2017) Shang, L., Cheng, Y. & Zhao, Y. 2017 Emerging droplet microfluidics. Chem. Rev. 117 (12), 7964–8040.
  * Sharma et al. (1995) Sharma, A., Kishore, C. S., Salaniwal, S. & Ruckenstein, E. 1995 Nonlinear stability and rupture of ultrathin free films. Phys. Fluids 7 (8), 1832–1840.
  * Sharma & Ruckenstein (1988) Sharma, A. & Ruckenstein, E. 1988 Effects of surfactants on wave-induced drainage of foam and emulsion films. Colloid Polym. Sci. 266, 60–69.
  * Sommeria (1986) Sommeria, J. 1986 Experimental study of the two-dimensional inverse energy cascade in a square box. J. Fluid Mech. 170, 139–168.
  * Sommeria & Moreau (1982) Sommeria, J. & Moreau, R. 1982 Why, how, and when, mhd turbulence becomes two-dimensional. J. Fluid Mech. 118, 507––518.
  * Suslov et al. (2017) Suslov, S. A., Pérez-Barrera, J. & Cuevas, S. 2017 Electromagnetically driven flow of electrolyte in a thin annular layer: axisymmetric solutions. J. Fluid Mech. 828, 573–600.
  * Van De Fliert et al. (1995) Van De Fliert, B. W., Howell, P. D. & Ockenden, J. R. 1995 Pressure-driven flow of a thin viscous sheet. J. Fluid Mech. 292, 359–376.
  * Williams et al. (1997) Williams, B. S., Marteau, D. & Gollub, J. P. 1997 Mixing of a passive scalar in magnetically forced two-dimensional turbulence. Phys. Fluids 9 (7), 2061–2080.
  * Wu et al. (1995) Wu, X.-L., Martin, B., Kellay, H. & Goldburg, W. I. 1995 Hydrodynamic convection in a two-dimensional Couette cell. Phys. Rev. Let. 75, 236–239.