2 Problem formulation

Figure 1.

Schematic representation of a viscous liquid film flow down an inclined plane wall. Throughout the study, the waves travel in the positive $x$ -direction.

We consider a viscous liquid film that flows down a plane wall inclined at angle $\unicode[STIX]{x1D6FD}$ to the horizontal. The film is exposed to an electric field which acts in the direction normal to the wall and which is uniform with strength $E_{0}$ at infinity, as shown in figure 1. The fluid is taken to be either a perfect conductor or a perfect dielectric with relative permittivity  $\unicode[STIX]{x1D700}_{p}$ . We use Cartesian coordinates $(x,y)$ with $x$ and $y$ measuring the distance along the wall and normal to it (pointing into the liquid), respectively. The dimensionless Navier–Stokes and continuity equations are

(2.1a,b ) $$\begin{eqnarray}\displaystyle Re(\boldsymbol{u}_{t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u})=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\boldsymbol{G},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{G}=2(\sin \unicode[STIX]{x1D6FD},-\cos \,\unicode[STIX]{x1D6FD})$ is the dimensionless gravity, $\boldsymbol{u}=(u,v)$ is the film velocity and $p$ is the film pressure; without loss of generality, we assume that the pressure in the air is zero. Distances have been made dimensionless using the undisturbed film thickness $h_{0}$ as the length scale, and $2\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C}gh_{0}$ has been used as the time scale, where $\unicode[STIX]{x1D707}$ and $\unicode[STIX]{x1D70C}$ are, respectively, the dynamic viscosity and density of the liquid, and $g$ is gravity. We have used $\unicode[STIX]{x1D70C}gh_{0}/2$ as the pressure scale and, below, we will use $E_{0}h_{0}$ as the scale for the electric potential. The Reynolds number is defined to be $Re=\unicode[STIX]{x1D70C}h_{0}U_{0}/\unicode[STIX]{x1D707}$ , where $U_{0}=\unicode[STIX]{x1D70C}gh_{0}^{2}/2\unicode[STIX]{x1D707}$ is the surface speed of a flat Nusselt film of thickness $h_{0}$ flowing down a vertical wall.

The no-slip condition at the wall is

(2.2) $$\begin{eqnarray}\displaystyle \boldsymbol{u}=\mathbf{0}\quad \text{at}~y=0, & & \displaystyle\end{eqnarray}$$

and the kinematic condition at the free surface $y=h(x,t)$ is

(2.3) $$\begin{eqnarray}\displaystyle v=h_{t}+uh_{x}. & & \displaystyle\end{eqnarray}$$

The tangential and normal dynamic stress conditions at the free surface are

(2.4a,b ) $$\begin{eqnarray}\displaystyle \boldsymbol{t}\boldsymbol{\cdot }\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}=0,\quad \boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}=-\frac{\unicode[STIX]{x1D705}}{Bo}+2We\,\boldsymbol{n}\boldsymbol{\cdot }(\unicode[STIX]{x1D648}_{2}-\unicode[STIX]{x1D648}_{1})\boldsymbol{\cdot }\boldsymbol{n}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{t}$ and $\boldsymbol{n}$ are the unit tangent and normal vectors at the free surface, respectively, the latter pointing into the liquid, $\unicode[STIX]{x1D705}$ is the free-surface curvature taken to be positive when the surface is concave downwards, and $\unicode[STIX]{x1D748}=-p\unicode[STIX]{x1D644}+(\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\text{T}})/2$ is the Newtonian stress tensor in the liquid. The Maxwell stress tensor is defined to be $\unicode[STIX]{x1D648}_{j}=\boldsymbol{E}_{j}\boldsymbol{E}_{j}-\unicode[STIX]{x1D644}|\boldsymbol{E}_{j}|^{2}/2$ , where $\boldsymbol{E}_{j}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}_{j}$ is the electric field with electric potential $\unicode[STIX]{x1D711}_{j}$ in region $j=1$ (the film) and $j=2$ (the air). The Bond number $Bo$ and electric Weber number $We$ are given by

(2.5a,b ) $$\begin{eqnarray}\displaystyle Bo=\frac{\unicode[STIX]{x1D70C}gh_{0}^{2}}{2\unicode[STIX]{x1D6FE}},\quad We=\left(1-\frac{1}{\unicode[STIX]{x1D700}_{p}}\right)^{2}\frac{\unicode[STIX]{x1D700}_{a}E_{0}^{2}}{\unicode[STIX]{x1D70C}gh_{0}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}$ is the surface tension coefficient, and $\unicode[STIX]{x1D700}_{p}=\unicode[STIX]{x1D700}_{1}/\unicode[STIX]{x1D700}_{2}$ , where $\unicode[STIX]{x1D700}_{j}$ is the electric permittivity of the film ( $j=1$ ) and the air ( $j=2$ ).

In the film and in the air above the film, the electric potentials satisfy Laplace’s equation, $\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D711}_{j}=0$ , and at the film surface they satisfy the boundary conditions,

(2.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D711}_{1}=\unicode[STIX]{x1D711}_{2},\quad \unicode[STIX]{x1D700}_{p}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}_{1}=\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}_{2}. & & \displaystyle\end{eqnarray}$$

At the wall $\unicode[STIX]{x1D711}_{1}=0$ , and far above the film the electric field is uniform so that $\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}_{2}\rightarrow -\boldsymbol{j}$ as $y\rightarrow \infty$ , where $\boldsymbol{j}$ is a unit vector pointing in the $y$ -direction.

3 Long-wave model

Assuming that the interfacial deformation wavelength $\unicode[STIX]{x1D706}$ is long compared to the undisturbed film thickness $h_{0}$ , so that the thin film parameter $\unicode[STIX]{x1D6FF}\equiv h_{0}/\unicode[STIX]{x1D706}\ll 1$ , we can derive the model equation

(3.1) $$\begin{eqnarray}\displaystyle h_{t}+\left[h^{3}\left(\frac{2\sin \unicode[STIX]{x1D6FD}}{3}-\frac{2\cos \unicode[STIX]{x1D6FD}}{3}h_{x}+\frac{1}{3Bo}h_{xxx}+\frac{2}{3}We\,{\mathcal{H}}[h_{xx}]\right)\right]_{x}=0. & & \displaystyle\end{eqnarray}$$

Here ${\mathcal{H}}$ is the Hilbert transform operator defined by ${\mathcal{H}}(f)=(1/\unicode[STIX]{x03C0})\text{p.v.}\int _{-\infty }^{\infty }f(\tilde{x})/(x-\tilde{x})$ $\text{d}\tilde{x}$ , where p.v. denotes the principal value. For any finite value of $\unicode[STIX]{x1D700}_{p}$ , the film behaves as a perfect dielectric. Mathematically, the perfect-conductor case is recovered in the limit $\unicode[STIX]{x1D700}_{p}\rightarrow \infty$ . The equation is valid provided that $Bo=O(\unicode[STIX]{x1D6FF}^{2})$ , $We=O(\unicode[STIX]{x1D6FF}^{-1})$ , and then either both $\cot \unicode[STIX]{x1D6FD}=O(1)$ and $Re=O(\unicode[STIX]{x1D6FF})$ or both $\cot \unicode[STIX]{x1D6FD}=O(\unicode[STIX]{x1D6FF}^{-1})$ and $Re=O(1)$ . It is derived by a systematic asymptotic procedure in the limit $\unicode[STIX]{x1D6FF}\rightarrow 0$ as described in Kalliadasis et al. (Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012) (see also Tseluiko & Papageorgiou Reference Tseluiko and Papageorgiou2006; Tseluiko, Blyth & Papageorgiou Reference Tseluiko, Blyth and Papageorgiou2013). Note that for $\cot \unicode[STIX]{x1D6FD}=O(\unicode[STIX]{x1D6FF}^{-1})$ the equation is a leading-order approximation, whereas for $\cot \unicode[STIX]{x1D6FD}=O(1)$ it combines both leading-order and first-order contributions.

If it is additionally assumed that $h(x,t)=1+\unicode[STIX]{x1D703}(x,t)$ , where $\unicode[STIX]{x1D703}=O(\unicode[STIX]{x1D6FF})$ , and a Galilean transformation is made to a moving frame of reference, namely $x\mapsto x+2(\sin \unicode[STIX]{x1D6FD})t$ , we obtain the non-local KS equation (see also Lin et al. Reference Lin, Pradas, Kalliadasis, Papageorgiou and Tseluiko2015),

(3.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703}_{t}+4(\sin \unicode[STIX]{x1D6FD})\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}_{x}-\frac{2}{3}(\cos \unicode[STIX]{x1D6FD})\,\unicode[STIX]{x1D703}_{xx}+\frac{2}{3}We{\mathcal{H}}[\unicode[STIX]{x1D703}_{xxx}]+\frac{1}{3Bo}\unicode[STIX]{x1D703}_{xxxx}=0. & & \displaystyle\end{eqnarray}$$

3.1 Travelling-wave solutions

To analyse travelling-wave solutions, which often characterise falling film flows, we introduce in (3.1) a moving frame coordinate via the mapping $x\mapsto x+c^{\ast }t$ , where $c^{\ast }$ is the wave speed, and look for stationary solutions $h^{\ast }$ in this frame. For solitary pulses, we will demand that the free-surface height approaches the Nusselt flat-film solution away from the centre of the pulse, i.e. $h^{\ast }\rightarrow 1$ as $x\rightarrow \pm \infty$ . Solutions can be obtained using continuation techniques. The numerical computations have spectral accuracy and were done in Matlab. We consider solutions on a periodic domain, $x\in [-L,L]$ , and can obtain solitary-pulse solutions in the limit $L\rightarrow \infty$ (when it is possible to continue solution branches to infinitely large $L$ ). We note that we need to additionally impose conditions breaking the translational and ‘volume’ symmetries. (The ‘volume’ symmetry is associated with the fact that, by slightly changing the volume $\int _{-L}^{L}h^{\ast }\,\text{d}x$ , we obtain a different solution.) To break the translational symmetry, we impose the condition $h_{x}^{\ast }=0$ at $x=0$ , and to break the ‘volume’ symmetry, we impose the condition $h^{\ast }=1$ at $x=L$ .

To start the continuation, we use a small-amplitude nearly sinusoidal solution of a nearly cutoff wavelength that can be found by a standard linear stability analysis. For acute inclination angles, gravity is stabilising, and the flat-film solution is linearly unstable only when the electric field is sufficiently strong, and in particular can be shown that this occurs when $We>We_{c}=\sqrt{2\cos \unicode[STIX]{x1D6FD}/Bo}$ . The range of unstable wavenumbers is $|k|\in (\unicode[STIX]{x1D705}_{-},\unicode[STIX]{x1D705}_{+})$ , where

(3.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D705}_{\pm }=Bo\,We\pm \sqrt{Bo^{2}We^{2}-2Bo\cos \unicode[STIX]{x1D6FD}}. & & \displaystyle\end{eqnarray}$$

For obtuse inclination angles, i.e. $\unicode[STIX]{x1D6FD}\in (\unicode[STIX]{x03C0}/2,\unicode[STIX]{x03C0})$ , $\cos \unicode[STIX]{x1D6FD}<0$ , the flat-film solution is linearly unstable due to the effect of gravity. The range of unstable wavenumbers is $|k|\in (0,\unicode[STIX]{x1D705}_{+})$ .

Figure 2.

Case of $\unicode[STIX]{x1D6FD}=0.95\unicode[STIX]{x03C0}$ , $Bo=0.005$ . Solid lines are for $We=0$ and dashed lines are for $We=3$ : (a) bifurcation diagram for the norm $\Vert h^{\ast }-1\Vert$ against domain size $L$ . The circles indicate the starting points for the continuations; (b) pulse profiles for large $L$ ; (c) shifted profiles $h^{\ast }-1$ for large $L$ shown on a logarithmic scale. Also included are dot-dashed lines showing inverse cube decay in $|x|$ for comparison purposes.

Figure 2 shows the results of numerical continuation for an obtuse inclination angle, $\unicode[STIX]{x1D6FD}=0.95\unicode[STIX]{x03C0}$ , when $Bo=0.005$ and $We=0$ and $3$ . A Stuart–Landau analysis suggests that the bifurcation at $L=L_{+}=\unicode[STIX]{x03C0}/\unicode[STIX]{x1D705}_{+}$ , both for $We=0$ and for $We=3$ , is supercritical, and the numerical results corroborate this. Both the norm and the speed converge to constants as $L$ increases. Also, for $We=3$ both the norm and the speed are greater than those for $We=0$ . Both travelling waves are single-hump solitary pulses preceded by capillary ripples of decaying amplitudes. The left tail of the pulse for $We=0$ decays monotonically to 1 as $L\rightarrow -\infty$ . However, for $We=3$ , the left tail exhibits a depression just upstream of the main hump before monotonically approaching 1. Also, the amplitude of the pulse is larger when the electric field is switched on, while the amplitude of the ripples is not so significantly affected by the electric field. The semi-log plots show that the electric field dramatically affects the decay rate of the tails, namely the decay is exponential when $We=0$ , but for $We=3$ it becomes much slower in the far field.

For acute inclination angles, the solutions of branches starting at the cutoff domain half-sizes, $L_{\pm }=\unicode[STIX]{x03C0}/\unicode[STIX]{x1D705}_{\pm }$ do not always converge to single-hump solitary pulses. In fact, we find that only the solutions of the branch starting at $L_{-}$ converge to a single-hump solitary pulse when the electric Weber number is greater than a certain value that is greater than the linear instability threshold value $We_{c}$ . We demonstrate this in figures 3 and 4 showing the results for $\unicode[STIX]{x1D6FD}=0.25\unicode[STIX]{x03C0}$ and $Bo=0.01$ and for $We=12.5$ and $13.5$ , respectively. Both of these values of the electric Weber number are above the linear instability threshold, which is $We_{c}\approx 11.89$ in this case. The bifurcations at $L=L_{+}$ and $L=L_{-}$ produce side branches which move toward larger and smaller values of $L$ , respectively, for $We=12.5$ . This is again in agreement with a Stuart–Landau analysis. For $We=12.5$ , the points corresponding to $L_{+}\approx 19.21$ and $L_{-}\approx 36.32$ are connected by a single branch of periodic travelling-wave solutions. For $We=13.5$ , there are two distinct solution branches emanating from the points corresponding to $L_{+}\approx 15.80$ and $L_{-}\approx 44.18$ , and both branches apparently extend to infinitely large values of $L$ . Both side branches move immediately toward larger values of $L$ in agreement with a Stuart–Landau analysis. For both branches, the norm $\Vert h^{\ast }-1\Vert$ and the speed $c^{\ast }$ tend to constant values as $L$ increases. The wave profile in (b) apparently converges to a single-hump solitary pulse preceded by capillary ripples of decaying amplitude, as for the case of an obtuse inclination angle. Also, like for the obtuse inclination angle, the left tail of the pulse is not monotonically decaying. There is a depression upstream of the pulse, and only then the tail tends to 1 monotonically as $x\rightarrow -\infty$ . The wave profile in (c) does not correspond to a single-hump pulse. Instead, it has the shape of a double-hollow negative solitary pulse in the terminology of Chang and co-workers (Chang Reference Chang1994; Chang & Demekhin Reference Chang and Demekhin2002) and with capillary ripples upstream and an elevation downstream, and, moreover, it travels at a speed that is smaller than that for the single-hump pulse.

Figure 3.

Case of $\unicode[STIX]{x1D6FD}=0.25\unicode[STIX]{x03C0}$ , $Bo=0.01$ , $We=12.5$ (a) bifurcation diagram for $\Vert h^{\ast }-1\Vert$ against domain size $L$ . (b,c) Wave profiles with norm $\Vert h^{\ast }-1\Vert =0.2$ , indicated by the left and right squares in (a), respectively.

Figure 4.

Case of $\unicode[STIX]{x1D6FD}=0.25\unicode[STIX]{x03C0}$ , $Bo=0.01$ , $We=13.5$ . (a) Bifurcation diagram for $\Vert h^{\ast }-1\Vert$ against domain size $L$ . The solid and dashed lines show bifurcation branches emerging from $L_{+}=15.80$ and $L_{-}=44.18$ , respectively. (b,c) Pulse profiles for large $L$ on the branches emerging from $L_{+}$ and $L_{-}$ , respectively.

Figure 5.

Case of $\unicode[STIX]{x1D6FD}=0.25\unicode[STIX]{x03C0}$ , $Bo=0.01$ . (a) Solution branch showing the norm $\Vert h^{\ast }-1\Vert$ against electric Weber number $We$ . (b,c) Pulse profiles at $We=18$ , indicated by the empty and filled circle symbols respectively in (a). The wave speeds are $c^{\ast }=2.22$ (b), 2.15 (c).

Finally, to investigate in more detail the influence of an electric field on single-hump solitary pulses, we consider such a pulse for $\unicode[STIX]{x1D6FD}=0.25\unicode[STIX]{x03C0}$ , $Bo=0.01$ on the domain $[-300,300]$ , and perform continuation in the electric Weber number, $We$ . The results are shown in figure 5. The branches have turning points at $We=We^{\dagger }\approx 12.77$ , which indicates that single-hump solitary pulses do not exist for $We<We^{\dagger }$ . The lower part of the branch in (a) of figure 5 corresponds to single-hump pulses. The upper branch in (a) corresponds to double-hump pulses. Both for single- and double-hump pulses the norm and the speed monotonically increase as $We$ increases.

A straightforward far-field analysis for the non-electrified case reveals that if $c^{\ast }>2\sin \unicode[STIX]{x1D6FD}$ we obtain monotonic decay on the upstream side of the pulse and non-monotonic, oscillatory decay on the downstream side of the tail. Referring back to the scales we have used to non-dimensionalise the problem, this inequality corresponds to the physical pulse speed being greater than the speed of small-amplitude linear long waves, which themselves propagate at twice the surface speed of a flat film on an inclined plane, namely $2U_{0}\sin \unicode[STIX]{x1D6FD}$ (e.g. Benjamin Reference Benjamin1957). If instead $c^{\ast }<2\sin \unicode[STIX]{x1D6FD}$ , we will have oscillatory decay on the upstream side of a pulse and monotonic decay on the downstream side.

With an electric field present, the decay becomes algebraic (see also Lin et al. Reference Lin, Pradas, Kalliadasis, Papageorgiou and Tseluiko2015). Assuming that the pulse has a non-zero ‘mass’, i.e. $\int _{-\infty }^{\infty }(h^{\ast }-1)\,\text{d}x\neq 0$ , it can be shown that $(h^{\ast }-1)\propto 1/x^{3}$ as $x\rightarrow \pm \infty$ .

3.2 Weak-interaction theory for solitary pulses

We generalise the earlier work of Lin et al. (Reference Lin, Pradas, Kalliadasis, Papageorgiou and Tseluiko2015) to the case of the general one-dimensional evolution equation with translational symmetry, $\unicode[STIX]{x2202}_{t}u=\unicode[STIX]{x2202}_{x}(G[u])$ , for some $G$ such that $G[h^{\ast }-\unicode[STIX]{x1D6FC}]=0$ , where $h^{\ast }\rightarrow \unicode[STIX]{x1D6FC}$ as $|x|\rightarrow \infty$ . The work of Lin et al. (Reference Lin, Pradas, Kalliadasis, Papageorgiou and Tseluiko2015) built on earlier studies by Elphick, Meron & Spiegel (Reference Elphick, Meron and Spiegel1990), Elphick et al. (Reference Elphick, Ierley, Regev and Spiegel1991), Balmforth (Reference Balmforth1995), Chang & Demekhin (Reference Chang and Demekhin2002) for the case of local equations. Certain technical details in these studies were re-examined by Pradas, Tseluiko & Kalliadasis (Reference Pradas, Tseluiko and Kalliadasis2011) and Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2014).

We consider (3.1) in a frame moving with the speed $c^{\ast }$ of a pulse solution. We represent a solution of (3.1) as a superposition of $N$ well-separated quasi-stationary pulses located at $x_{1}(t)<\cdots <x_{N}(t)$ and a small overlap function $\hat{u}$ . Suppressing the details in the interest of brevity, we obtain the dynamical system describing the locations of the pulses (see Lin et al. Reference Lin, Pradas, Kalliadasis, Papageorgiou and Tseluiko2015)

(3.4) $$\begin{eqnarray}\displaystyle {\dot{x}}_{k}=-\unicode[STIX]{x1D6F1}_{k}\left[\unicode[STIX]{x2202}_{x}\left(G\left[\mathop{\sum }_{|x_{i}-x_{k}|<b(\unicode[STIX]{x1D716})}u_{i}\right]\right)\right],\quad k=1,\ldots ,N, & & \displaystyle\end{eqnarray}$$

where $b(\unicode[STIX]{x1D716})\gg \unicode[STIX]{x1D716}^{-1/p}$ is the range over which the interactions must be taken into account. Here, $\unicode[STIX]{x1D6F1}_{k}$ is a projection operator defined by $\unicode[STIX]{x1D6F1}_{k}[f]=\int _{-\infty }^{\infty }f\unicode[STIX]{x1D713}_{k}\,\text{d}x$ , where $\unicode[STIX]{x1D713}_{k}$ is the zero adjoint eigenfunction for the operator $\unicode[STIX]{x2202}_{x}\circ \mathscr{D}_{u_{k}}G$ normalised so that $\int _{-\infty }^{\infty }u_{kx}\unicode[STIX]{x1D713}_{k}\,\text{d}x=1$ . Note that $\mathscr{D}_{u_{k}}G$ denotes the Fréchet derivative of operator $G$ at $u_{k}$ . We note that in (3.4), for the case of an electrified film, it is necessary to include long-range interactions over more than just the immediately neighbouring pulses as a direct result of the algebraic decay of the pulse tails.

Let us consider a two-pulse system, $N=2$ . From (3.4), we may derive the evolution equation for the pulse separation

(3.5) $$\begin{eqnarray}\displaystyle \dot{l}=\tilde{P}(l),\quad \tilde{P}(l)\equiv \int _{-\infty }^{\infty }\unicode[STIX]{x2202}_{x}(G[u^{\ast }(x+l)-u^{\ast }(x-l)])\unicode[STIX]{x1D713}(x)\,\text{d}x. & & \displaystyle\end{eqnarray}$$

Depending on the initial condition, the two pulses may attract ( $\tilde{P}(l)<0$ ) or repel ( $\tilde{P}(l)>0$ ) each other. The pulses can also form bound states, the separation distances for which are zeros of $\tilde{P}(l)$ . Figures 6(a) and 6(b) show $\tilde{P}(l)$ for $\unicode[STIX]{x1D6FD}=0.95\unicode[STIX]{x03C0}$ and $Bo=0.005$ , and for $We=0$ and $2.25$ , respectively. Apparently, there exist infinitely many bound states for $We=0$ , while there exist only eight possible bound states for $We=2.25$ . To understand this, we note that for $We=0$ , it can be shown that as $l\rightarrow \infty$ , $\tilde{P}$ decays either monotonically or in an oscillatory manner. If the decay is monotonic, there are zero or a finite number of positive solutions to $\tilde{P}(l)=0$ implying zero or a finite number of two-pulse bound states. If the decay is oscillatory, there is a countably infinite number of solutions to $\tilde{P}(l)=0$ , and, hence, an infinite number of two-pulse bound states.

We note that in the case of no electric field, the travelling-wave form of (3.1) is local and can be integrated once to yield a three-dimensional dynamical system. Then, a single-pulse solution corresponds to a homoclinic orbit in the three-dimensional phase space for this system. In such a case, the same conclusion on the number of bound states in existence can be reached using Shil’nikov’s theorem (see, for example, Glendinning & Sparrow Reference Glendinning and Sparrow1984). This gives a criterion for the number of subsidiary homoclinic orbits which correspond to two-pulse bound states, which exactly coincides with that found via the monotonic/oscillatory decay argument given above. However, when $We>0$ , the travelling-wave form of (3.1) is non-local and the Shil’nikov-type approach is not applicable. Instead, the weak-interaction theory developed here can be used to analyse bound states. Then, it can be shown that $\tilde{P}(l)$ decays monotonically algebraically to zero as $l\rightarrow \infty$ , and hence there exist only a finite number of two-pulse bound states.

Figure 6.

Value of $\tilde{P}(l)$ against $l$ for $\unicode[STIX]{x1D6FD}=0.95\unicode[STIX]{x03C0}$ and $Bo=0.005$ , and for (a) $We=0$ and (b) $We=2.25$ . Bound states are marked as filled circles (stable) and empty circles (unstable). The inset in (b) shows the behaviour for large $l$ , which becomes monotonic.

To validate the theoretically predicted pulse dynamics (3.5), we numerically simulated the dynamics of a two-pulse system using a superposition of two pulses as an initial condition for various parameter values and initial separation distances. We found excellent agreement between the simulations and the theory. Similarly good agreement was found for simulations with three and four pulses. For a physical explanation of pulse attraction and repulsion, see Duprat et al. (Reference Duprat, Giorgiutti-Dauphiné, Tseluiko, Saprykin and Kalliadasis2009).

4 Fully nonlinear solutions for Stokes flow

We now study fully nonlinear pulse solutions at zero Reynolds number. We discuss the case of a perfect-conductor film, that is $\unicode[STIX]{x1D700}_{p}\rightarrow \infty$ . The extension to the case of a liquid dielectric is straightforward (e.g. Tseluiko et al. Reference Tseluiko, Blyth, Papageorgiou and Vanden-Broeck2008).

We first reformulate the problem using the boundary-integral method. The flow is assumed to be periodic in $x$ with half-period $L$ . We decompose the steady flow velocity $\boldsymbol{u}$ and stress tensor $\unicode[STIX]{x1D748}$ into a basic flow part and a disturbance part, writing $\boldsymbol{u}=\boldsymbol{u}^{B}+\boldsymbol{u}^{D}$ and $\unicode[STIX]{x1D748}=\unicode[STIX]{x1D748}^{B}+\unicode[STIX]{x1D748}^{D}$ . In the travelling frame, the basic flow details are

(4.1a,b ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}^{B}=(y(2-y)\sin \unicode[STIX]{x1D6FD}-c^{\ast },0),\quad \unicode[STIX]{x1D748}^{B}=-p_{a}\unicode[STIX]{x1D644}+2(1-y)\left(\begin{array}{@{}cc@{}}-\cos \,\unicode[STIX]{x1D6FD} & \sin \unicode[STIX]{x1D6FD}\\ \sin \unicode[STIX]{x1D6FD} & -\cos \,\unicode[STIX]{x1D6FD}\end{array}\right), & & \displaystyle\end{eqnarray}$$

where $p_{a}$ is the constant dimensionless ambient pressure above the film, and the remaining symbols have been defined in § 3. Setting $c^{\ast }=0$ in (4.1), we recover the classical Nusselt solution for unidirectional flow down an inclined plane. We note that according to the decomposition described, the disturbance velocity $\boldsymbol{u}^{D}$ vanishes at the wall, $y=0$ .

Following the boundary-integral formalism, we obtain the integral equation for the disturbance velocity and disturbance traction $\boldsymbol{f}^{D}\equiv \unicode[STIX]{x1D748}^{D}\boldsymbol{\cdot }\boldsymbol{n}$ at a point $\boldsymbol{x}_{0}=(x_{0},y_{0})$ located on the free surface,

(4.2) $$\begin{eqnarray}\displaystyle 2\unicode[STIX]{x03C0}u_{j}^{D}(\boldsymbol{x}_{0})=-\int _{F}G_{ij}(\boldsymbol{x},\boldsymbol{x}_{0})f_{j}^{D}(\boldsymbol{x})\,\text{d}l(\boldsymbol{x})+\text{p.v.}\int _{F}u_{i}^{D}(\boldsymbol{x})\unicode[STIX]{x1D61B}_{ijk}(\boldsymbol{x},\boldsymbol{x}_{0})n_{k}(\boldsymbol{x})\,\text{d}l(\boldsymbol{x}), & & \displaystyle\end{eqnarray}$$

where p.v. denotes the principal value, $F$ denotes one period of the free surface, $l$ is arc length along $F$ , $\boldsymbol{n}$ is the unit normal at the free surface pointing into the fluid, $G_{ij}(\boldsymbol{x},\boldsymbol{x}_{0})$ is the periodic Green’s function for Stokes flow, which vanishes when $\boldsymbol{x}$ is located on the wall at $y=0$ , and $\unicode[STIX]{x1D61B}_{ijk}(\boldsymbol{x},\boldsymbol{x}_{0})$ is the corresponding stress tensor. Closed form expressions for $\boldsymbol{G}$ and $\boldsymbol{T}$ can be found in Pozrikidis (Reference Pozrikidis2002). The kinematic condition at the free surface in the travelling frame requires that

(4.3) $$\begin{eqnarray}\displaystyle \boldsymbol{u}^{D}\boldsymbol{\cdot }\boldsymbol{n}=-\boldsymbol{u}^{B}\boldsymbol{\cdot }\boldsymbol{n}. & & \displaystyle\end{eqnarray}$$

Following the boundary-integral method for Laplace’s equation (e.g. Pozrikidis Reference Pozrikidis2002), we obtain the integral equation for the electric potential,

(4.4) $$\begin{eqnarray}\displaystyle \frac{1}{2}\unicode[STIX]{x1D711}_{2}(\boldsymbol{x}_{0})=-\int _{F}G(\boldsymbol{x},\boldsymbol{x}_{0})\boldsymbol{n}(\boldsymbol{x})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}_{2}(\boldsymbol{x})\,\text{d}l(\boldsymbol{x})-y_{0}+B, & & \displaystyle\end{eqnarray}$$

where $B$ is an a priori unknown constant, and $G$ is the singly periodic upward-biased Green’s function with half-period $L$ , given by (see Pozrikidis Reference Pozrikidis2002, p. 261)

(4.5) $$\begin{eqnarray}\displaystyle G(\boldsymbol{x},\boldsymbol{x}_{0})=-\frac{1}{4\unicode[STIX]{x03C0}}\log [2(\cosh [k(y-y_{0})]-\cos [k(x-x_{0})])]-\frac{y-y_{0}}{4L}. & & \displaystyle\end{eqnarray}$$

By integrating Laplace’s equation for $\unicode[STIX]{x1D711}$ over one period of the semi-infinite region above $F$ , we obtain the integral condition

(4.6) $$\begin{eqnarray}\displaystyle \int _{F}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}_{2}\,\text{d}l(\boldsymbol{x})=2L, & & \displaystyle\end{eqnarray}$$

which is to be satisfied along with (4.4).

The problem for the electric field and the problem for the fluid flow are coupled together via the dynamic stress boundary condition at the free surface,

(4.7) $$\begin{eqnarray}\displaystyle \boldsymbol{f}^{D}=-\boldsymbol{f}^{B}-\frac{\unicode[STIX]{x1D705}}{Bo}\boldsymbol{n}+We|\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}|^{2}\boldsymbol{n}. & & \displaystyle\end{eqnarray}$$

To complete the formulation, we specify two more conditions akin to those imposed for the long-wave model. First, we remove the translational invariance by imposing the condition $h_{x}^{\ast }=0$ at $x=0$ , and to break the ‘volume’ symmetry, we impose the condition $h^{\ast }=1$ at $x=L$ . By fixing the location of the wave maximum and the film height in this way, we efficiently compute pulse solutions (where they exist) by increasing $L$ .

With these final two conditions, the coupled integral formulation comprising (4.2), (4.4) and (4.6), together with the kinematic condition (4.3) and the dynamic stress condition (4.7) are solved numerically using the boundary-element method by discretising one period of the a priori unknown free surface, $F$ , with a sequence of $N$ connected straight elements and treating the unknown variables as constants over the elements. In this way, we derive a set of nonlinear algebraic equations to be solved using Newton’s method for these unknown constant element values, the free-surface profile $h^{\ast }$ and the wave speed $c^{\ast }$ . Full details of the implementation of the method are provided by Tseluiko et al. (Reference Tseluiko, Blyth, Papageorgiou and Vanden-Broeck2008) and, in the interest of brevity, we do not repeat them here.

Figure 7.

Stokes flow computations with $\unicode[STIX]{x1D6FD}=0.95\unicode[STIX]{x03C0}$ , $Bo=0.005$ , $We=0$ : (a) travelling-wave branch for Stokes flow showing the wave height, $\max (h^{\ast })-\min (h^{\ast })$ , against $L$ , computed using $N=300$ boundary elements. (b) Numerical convergence for $L=150$ showing $\max (h^{\ast })$ against $N$ . (c) Pulse profile for $L=150$ and $N=300$ (solid line) with the long-wave predictions using (3.1) and (4.8) (dashed, dot-dashed lines, respectively) and the KS prediction (thin solid line). The wave speed is $c^{\ast }=0.457$ (boundary-element solution) and $c^{\ast }=0.449$ (long-wave solution to (3.1)) and $c^{\ast }=0.392$ (KS equation (3.2)).

A standard normal-mode analysis for small-amplitude periodic waves (see Blyth Reference Blyth2008) reveals that the cutoff wavenumber coincides with that found from the long-wave analysis. The properties of the neutral mode are used to construct an initial guess for the Newton iterations and to latch onto the travelling-wave solution branch. This is followed using continuation in $L$ , and in the case where $L$ may be increased indefinitely, we are able to compute solitary-pulse solutions. The solution space is determined by the three dimensionless parameters $Bo$ , $We$ and $\unicode[STIX]{x1D6FD}$ introduced in § 3.

We begin by computing fully nonlinear pulse solutions for no electric field, setting $We=0$ . In figure 7(a), we show the results of a sequence of calculations for increasing domain size $L$ . The wave height approaches a constant value, and a pulse solution is obtained, at around $L\approx 200$ . It was found that $N=300$ boundary elements were sufficient to obtain an accurate wave solution. The numerical convergence is demonstrated in figure 7(b), which shows the variation of the pulse maximum with the number of boundary elements $N$ at $L=150$ . The pulse itself is shown in figure 7(c) with a thick solid line. The long-wave prediction based on (3.1) and the weakly nonlinear prediction based on the KS equation (3.2) are also shown in this figure. The pulse speed is found to be $c^{\ast }=0.457$ for the boundary-element calculation and $c^{\ast }=0.449$ for the long-wave calculation, and so there is reasonable agreement between the two values.

Notably, the performance of the KS equation is rather poor, even at this small value of the Bond number. Nevertheless, we have confirmed that the long-wave and KS solutions converge as the Bond number is further reduced. Evidently the long-wave and boundary-element calculations agree very well with a small discrepancy around the maximum of the pulse. The discrepancy is significantly reduced by including higher-order terms in the long-wave model equation; specifically including second-order terms in the derivation we obtain the extended form of (3.1),

(4.8) $$\begin{eqnarray}\displaystyle h_{t}+q_{x}=0, & & \displaystyle\end{eqnarray}$$

where $q=q_{0}+q_{1}+q_{2}$ , with $q_{0}=O(1)$ , $q_{1}=O(\unicode[STIX]{x1D6FF})$ , $q_{2}=O(\unicode[STIX]{x1D6FF}^{2})$ . Specifically,

(4.9) $$\begin{eqnarray}\displaystyle q_{0}=\frac{2}{3}(\sin \unicode[STIX]{x1D6FD})h^{3}-\frac{2}{3}(\cos \unicode[STIX]{x1D6FD})h^{3}h_{x}+\frac{1}{3Bo}h^{3}h_{xxx}+\frac{2We}{3}\,h^{3}\,{\mathcal{H}}[h_{xx}], & & \displaystyle\end{eqnarray}$$

and $q_{1}=F_{1}+{\mathcal{E}}_{1}$ , where $F_{1}$ contains terms all of which vanish when $Re=0$ , and

(4.10) $$\begin{eqnarray}\displaystyle {\mathcal{E}}_{1}={\textstyle \frac{1}{3}}We\,h^{3}(E_{2x}+2h_{x}h_{xx}), & & \displaystyle\end{eqnarray}$$

(see Tseluiko & Papageorgiou Reference Tseluiko and Papageorgiou2010) and

(4.11) $$\begin{eqnarray}\displaystyle q_{2} & = & \displaystyle h^{3}\left(\sin \unicode[STIX]{x1D6FD}\left(2hh_{xx}+\frac{14}{3}h_{x}^{2}\right)-2\cos \unicode[STIX]{x1D6FD}\left(4hh_{x}h_{xx}+\frac{3}{5}h^{2}h_{xxx}+\frac{7}{3}h_{x}^{3}\right)\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{1}{Bo}\left(3hh_{x}h_{xxxx}-h_{x}h_{xx}^{2}+\frac{11}{6}h_{x}^{2}h_{xxx}+hh_{xx}h_{xxx}+\frac{3}{5}h^{2}h_{xxxxx}\right)\right)+F_{2}+{\mathcal{E}}_{2},\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $F_{2}$ contains terms all of which vanish when $Re=0$ . The electric-field contribution ${\mathcal{E}}_{2}$ must be found by considering the electric-field problem at the appropriate order, but is not needed for the present non-electrified case. Equation (4.8) is valid provided that $Bo=O(\unicode[STIX]{x1D6FF}^{2})$ , $We=O(\unicode[STIX]{x1D6FF}^{-1})$ , $\cot \unicode[STIX]{x1D6FD}=O(\unicode[STIX]{x1D6FF}^{-1})$ and $Re=O(1)$ . The pulse solution to this equation is shown in figure 7(c) with a dot-dashed line which almost coincides with the thick solid line representing the boundary-element solution.

Figure 8.

Stokes flow calculation for $\unicode[STIX]{x1D6FD}=0.75\unicode[STIX]{x03C0}$ with $Bo=0.3$ and $W_{e}=0$ , and with $L=16.0$ and $N=400$ . (a) Pulse profile and streamlines. (b) Logarithm of the pulse profile (solid line) with broken lines of gradient $0.576$ in $x<0$ and $-0.741$ in $x>0$ corresponding to the decay rate predictions (see appendix A). The wave speed is $c^{\ast }=3.82$ .

The pulse profile in figure 7(c) decays monotonically on the upstream side and has an oscillatory decay on the downstream side. Since the pulse speed is greater than the speed of linear long waves, namely $2\sin \unicode[STIX]{x1D6FD}=0.313$ , this is consistent with the predictions of the decay rate calculations in appendix A. These calculations predict that when the Bond number is not small the decay is monotonic both upstream and downstream of the pulse maximum. Figure 8(a) shows the pulse profile for $Bo=0.3$ and $\unicode[STIX]{x1D6FD}=0.75\unicode[STIX]{x03C0}$ . The pulse speed is $c^{\ast }=3.82$ . Evidently the decay is monotonic on both sides of the pulse. For these parameter values the calculation described in appendix A yields the upstream and downstream decay rates $0.576$ and $-0.741$ , respectively, and these show good agreement with the profile calculated using the boundary-element method, as can be seen in figure 8(b). The figure also shows streamlines inside the film in a frame of reference travelling at the speed of the pulse. These indicate the presence of a trapped eddy in the main part of the pulse. Solitary-wave eddies have recently been observed experimentally on a gravity-driven film at non-zero Reynolds number by Reck & Aksel (Reference Reck and Aksel2015).

Figure 9.

Electrified Stokes flow computation for $\unicode[STIX]{x1D6FD}=0.95\unicode[STIX]{x03C0}$ , $Bo=0.005$ , $We=6.25$ and $L=250$ , $N=500$ . (a) Pulse profile (solid line) with the long-wave prediction (dashed line) and the KS prediction (dot-dashed line). (b) A $\log$ – $\log$ plot of the pulse profile in $x>0$ and $x<0$ showing algebraic decay in the far field with the expected slope (dashed line) according to the asymptotic theory presented in appendix A. The wave speed is $c^{\ast }=1.09$ (boundary-element solution) and $c^{\ast }=0.88$ (long-wave solution) and $c^{\ast }=0.46$ (KS solution).

An electrified solitary-pulse solution is shown in figure 9(a). According to the theory of appendix A, the far-field decay of an electrified pulse is algebraic and so a wide computational domain and a large number of boundary elements are needed for an accurate computation. The prediction of the long-wave model equation (3.1) is also shown in the figure with a broken line. Once again, we see that there is good agreement over most of the pulse profile except near to the main peak. The visible difference between the Stokes calculation and the long-wave one at the pulse maximum is exacerbated in the presence of the electric field (compare the results in figure 7 c). As for the non-electrified case studied in figure 7(b,c), we would expect the agreement to improve on using the extended long-wave model equation (4.8). However, this would require computation of the corresponding electric-field contribution ${\mathcal{E}}_{3}$ . The decay rate of the pulse tails as $|x|\rightarrow \infty$ is investigated in figure 9(b). The broken lines shown in the figure indicate the algebraic decay rate expected from the asymptotic theory of appendix A. The excellent agreement provides strong evidence of algebraic decay in the far field and lends strong credence to the decay rate predicted by the asymptotic theory. The pulse speed determined from the boundary-element solution is $c^{\ast }=1.09$ ; this compares with $c^{\ast }=0.88$ obtained from the long-wave theory. We note that the pulse travels faster in the presence of the electric field. This trend is in line with the long-wave theory.

Figure 10.

Streamline patterns in a frame moving at the pulse speed for the case $\unicode[STIX]{x1D6FD}=0.95\unicode[STIX]{x03C0}$ with $Bo=0.005$ : (a,b) $We=0$ , (c,d) $We=6.25$ . We note that (b) and (d) are results for the Stokes model while (a) and (c) are for the long-wave model. The pulse profiles are shown with a thick line and the streamlines with thin lines.

Figure 10 depicts the streamline patterns for $\unicode[STIX]{x1D6FD}=0.95\unicode[STIX]{x03C0}$ and $Bo=0.005$ in a frame moving with the pulse for $We=0$ (a,b) and $We=6.25$ (c,d). Figure 10(a,c) shows results for the long-wave model and figure 10(b,d) shows the results of boundary-element calculations. The two different models give broadly similar results. Notably the electric field generates an eddy/recirculation zone inside the hump so that a quantity of fluid is transported along with the pulse.

Next, we discuss the computation of bound states for the Stokes equations. We construct an initial guess made from a superposition of a converged pulse solution and a duplicate of the same pulse separated by a nominated distance, $l$ , which is taken to coincide with that given by the long-wave theory of § 3.2. A good initial guess is required as the convergence of the Newton iterations is sensitive to the separation distance. Figure 11 shows four bound-state solutions for the case $Bo=0.005$ and $\unicode[STIX]{x1D6FD}=0.95\unicode[STIX]{x03C0}$ for $We=0$ or 2.25. For these parameter values and for the non-electrified case, the long-wave theory predicts the existence of twelve bound states with separation distances: 47.30, 69.74, 94.13, 117.62, 141.38, 165.05, 188.75, 212.42, 236.37, 260.01, 283.71, 307.41. We computed only a subset of these shown in the figure. The accuracy of each calculation was confirmed by varying the number of boundary elements and the size of the computational domain. For each case, the wave speed $c^{\ast }$ predicted by the boundary-element computation and by the long-wave theory are given in table 1. Evidently, the two are in good agreement. The wave speed is quite a lot smaller for a bound state than for a solitary pulse travelling alone. We also note that for both the boundary-element and the long-wave calculations, the separation distances, which in the case of Stokes flow are measured as the distance from maximum to maximum, decrease as the electric-field intensity is increased.

Figure 11.

Stokes flow bound-state computations for $\unicode[STIX]{x1D6FD}=0.95\unicode[STIX]{x03C0}$ , $Bo=0.005$ (with $L=175$ , $N=500$ ) for $We=0$ (thick solid lines) and $We=2.25$ (thin solid lines). The separation distances from maximum to maximum are (with electrified values in brackets): (a) 56.5 (50.5), (b) 67.0 (60.0), (c) 96.0 (84.0) and (d) 117.0 (111.0). The long-wave predictions for non-electrified (electrified) flows are shown with thick (thin) dashed lines. The long-wave bound-state separation distances are (with electrified values in brackets): (a) 47.3 (43.4), (b) 69.7 (62.5), (c) 94.1 (83.6) and (d) 117.6 (111.2). The wave speeds are given in table 1.

Table 1.

Wave speeds for the bound states in figure 11 quoted for the Stokes flow boundary-element computations ( $c_{BEM}^{\ast }$ ) and for the long-wave model, ( $c_{LW}^{\ast }$ ). Values in brackets are for the electrified cases with electric Weber numbers $We$ as given in the caption to figure 11.

In § 3.1, using the long-wave model, we found that in the case of electrified flow at an acute inclination angle, the travelling-wave branches which emerge from the two neutral stability points at $L=L_{\pm }$ either connect to form a continuous hoop from one to the other or else each branch continues independently to large $L$ ultimately producing a solitary-pulse solution. Figures 12 and 13 demonstrate that the same behaviour is observed for Stokes flow. These two figures use the same parameter values as for the long-wave calculations in figures 3 and 4. In figure 12 ( $We=12.5$ ), we see that the neutral stability points at $L_{+}\approx 19.21$ and $L_{-}\approx 36.32$ are connected by a single continuous travelling-wave branch. Sample profiles along the branch are also shown in the figure. In figure 13 ( $We=13.5$ ), two independent branches emerge from the neutral points at $L_{+}\approx 15.80$ and $L_{-}\approx 44.18$ , and continue until pulse solutions are finally attained for large $L$ . The pulse profiles are also shown in the figure.

Figure 12.

Stokes flow computations for $\unicode[STIX]{x1D6FD}=0.25\unicode[STIX]{x03C0}$ , $Bo=0.01$ , $We=12.5$ : (a) bifurcation diagram for $\Vert h^{\ast }-1\Vert$ against $L$ . (c,d) Wave profiles with $\Vert h^{\ast }-1\Vert =0.2427$ , indicated by the left and right squares in (a), respectively.

Figure 13.

Stokes flow computations for $\unicode[STIX]{x1D6FD}=0.25\unicode[STIX]{x03C0}$ , $Bo=0.01$ , $We=13.5$ : (a) bifurcation diagram for the $\Vert h^{\ast }-1\Vert$ against $L$ with branches emerging from $L_{+}\approx 15.80$ (solid line) and $L_{-}\approx 44.18$ (dashed line). (b,c) Pulse profiles for large $L$ on the branches emerging from $L_{+}$ and $L_{-}$ , respectively.

We present in figure 14 some examples of negative solitary pulses. In the non-electrified case shown in figure 14(a), the boundary-element solution is compared with the prediction for the same wave using the long-wave model. It is noticeable that, as for the elevation pulse shown in figure 7(c), the greatest discrepancy between the boundary-element and long-wave calculations is observed at the largest peak. The pulse speed is $c^{\ast }=0.259$ for the Stokes calculation and $c^{\ast }=0.257$ for the long-wave calculation. This is smaller than the speed of linear long waves $2\sin \unicode[STIX]{x1D6FD}=0.313$ . Evidently the negative pulse is also travelling more slowly than its elevation pulse counterpart for the same Bond number and inclination angle shown in figure 7(c), whose speed is $c^{\ast }=0.457$ .

Figure 14(b) shows the effect of increasing the electric-field intensity. Similar to the elevation pulses studied above, the electric field tends to deepen the depression and heighten the amplitude of the downstream oscillations. It also tends to lower the speed of the wave (values for the wave speed are quoted in the figure caption).

Figure 14.

Negative pulse profiles for Stokes flow with $Bo=0.005$ , $\unicode[STIX]{x1D6FD}=0.95\unicode[STIX]{x03C0}$ (with $L=250$ , $N=500$ ). (a) $We=0$ : boundary-element solution (solid line) and long-wave solution (dashed line). (b) Boundary-element solutions for $We=0$ , 1.0, 4.0, 9.0 and 16.0 with the depth of the main trough increasing with $We$ . The wave speeds are, respectively, $c^{\ast }=0.259$ , 0.255, 0.242, 0.217, 0.183.

5 Absolute/convective instability of pulse solutions

We now discuss the stability of the pulse solutions calculated in the previous sections. In particular, we classify pulses as being either convectively unstable or absolutely unstable. Suppose first that the flow supports pulses which are convectively unstable. Since such pulses can tolerate localised disturbances, in time-dependent solutions we would see the flow ultimately evolve into a state with an array of such pulses undergoing weak interactions. Suppose instead that the pulses supported by the flow are absolutely unstable. These can be destroyed by localised disturbances and consequently we would expect to see highly irregular dynamics in a general time-dependent simulation. For further discussion see Chang et al. (Reference Chang, Demekhin and Kopelevich1995), Chang & Demekhin (Reference Chang and Demekhin2002).

To study stability, we linearise about a pulse state by imposing a small perturbation in a reference frame which is moving with the pulse at speed $c^{\ast }$ , writing

(5.1) $$\begin{eqnarray}\displaystyle h(x,t)=h^{\ast }(x)+\unicode[STIX]{x1D702}(x,t), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D702}(x,t)$ is a small perturbation localised in space, and $h^{\ast }\rightarrow 1$ as $|x|\rightarrow \infty$ .

We consider first long-wave pulses which are solutions of the model equation (3.1). Substituting (5.1) into (3.1), written in a frame moving at the pulse speed $c^{\ast }$ , and ignoring nonlinear terms, we obtain

(5.2a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D702}={\mathcal{L}}[\unicode[STIX]{x1D702}],\quad {\mathcal{L}}\unicode[STIX]{x1D702}=(a_{0}\unicode[STIX]{x1D702}+a_{1}\unicode[STIX]{x1D702}_{x}+a_{2}{\mathcal{H}}[\unicode[STIX]{x1D702}_{xx}]+a_{3}\unicode[STIX]{x1D702}_{xxx})_{x}, & & \displaystyle\end{eqnarray}$$

with

(5.3) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle a_{0}(x)=c^{\ast }-2(\sin \unicode[STIX]{x1D6FD})h^{\ast 2}+2(\cos \unicode[STIX]{x1D6FD})h^{\ast 2}h_{x}^{\ast }-\frac{1}{Bo}h^{\ast 2}h_{xxx}^{\ast }-2We\,h^{\ast 2}{\mathcal{H}}[h_{xx}^{\ast }],\\ \displaystyle a_{1}(x)=\frac{2\cos \unicode[STIX]{x1D6FD}}{3}h^{\ast 3},\quad a_{2}(x)=-\frac{2}{3}We\,h^{\ast 3},\quad a_{3}(x)=-\frac{1}{3Bo}h^{\ast 3}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The solution of (5.2) can be written as (e.g. Chang, Demekhin & Kopelevich Reference Chang, Demekhin and Kopelevich1996; Lin et al. Reference Lin, Pradas, Kalliadasis, Papageorgiou and Tseluiko2015)

(5.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}(x,t)=\mathop{\sum }_{i}\text{e}^{\unicode[STIX]{x1D706}_{i}t}B_{i}\unicode[STIX]{x1D719}_{i}(x)+\int _{-\infty }^{\infty }\text{e}^{\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D705})t}B(\unicode[STIX]{x1D705})\unicode[STIX]{x1D719}(x,\unicode[STIX]{x1D705})\,\text{d}\unicode[STIX]{x1D705}, & & \displaystyle\end{eqnarray}$$

where the summation is over all the isolated eigenvalues $\unicode[STIX]{x1D706}_{i}$ (the discrete spectrum) with corresponding eigenfunctions $\unicode[STIX]{x1D719}_{i}(x)$ , that is the functions belonging to the null space of $\unicode[STIX]{x1D706}_{i}I-{\mathcal{L}}$ , and $B_{i}$ are constants. In fact, for the given model, we find numerically that there is only one isolated eigenvalue, which is real and negative. The same happens, for example, for the gKS equation (Chang et al. Reference Chang, Demekhin and Kopelevich1996; Tseluiko & Kalliadasis Reference Tseluiko and Kalliadasis2014). In the second integral, $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D705})$ is the essential spectrum of ${\mathcal{L}}$ which coincides precisely with the spectrum for a flat film as is the case with the established theory for ordinary differential operators (Edmunds & Evans Reference Edmunds and Evans1987; Pego & Weinstein Reference Pego and Weinstein1992). For a flat film of unit thickness this is given by

(5.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}=(c^{\ast }-2\sin \unicode[STIX]{x1D6FD})\text{i}\unicode[STIX]{x1D705}-\frac{2\cos \unicode[STIX]{x1D6FD}}{3}\,\unicode[STIX]{x1D705}^{2}+\frac{2We}{3}\,|\unicode[STIX]{x1D705}|^{3}-\frac{\unicode[STIX]{x1D705}^{4}}{3Bo}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D705}\in \mathbb{R}$ . In (5.4), $\unicode[STIX]{x1D719}(x,\unicode[STIX]{x1D705})$ are functions in to the null space of $\unicode[STIX]{x1D70E}I-{\mathcal{L}}$ , and (5.5) implies that $\unicode[STIX]{x1D70E}(-\unicode[STIX]{x1D705})=\overline{\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D705})}$ and also $\unicode[STIX]{x1D719}(x,-\unicode[STIX]{x1D705})=\overline{\unicode[STIX]{x1D719}(x,\unicode[STIX]{x1D705})}$ , where the overline denotes complex conjugation. For a real perturbation $\unicode[STIX]{x1D702}$ , the coefficients $B(\unicode[STIX]{x1D705})$ satisfy $B(-\unicode[STIX]{x1D705})=\overline{B(\unicode[STIX]{x1D705})}$ .

Since the discrete spectrum contains only $\unicode[STIX]{x1D706}_{1}$ , and $\unicode[STIX]{x1D706}_{1}$ is negative, the first term in (5.4) decays to zero as $t\rightarrow \infty$ , and the convective/absolute nature of the instability is determined by the second, integral term, which we label $I_{c}$ . From (5.5), if $We^{2}Bo\leqslant 2\cos \unicode[STIX]{x1D6FD}$ then $\text{Re}(\unicode[STIX]{x1D70E})\leqslant 0$ and the flow is spectrally stable. However, if $We^{2}Bo>2\cos \unicode[STIX]{x1D6FD}$ then $\text{Re}(\unicode[STIX]{x1D70E})>0$ if $|\unicode[STIX]{x1D705}|\in (\unicode[STIX]{x1D705}_{1},\unicode[STIX]{x1D705}_{2})$ , where

(5.6a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D705}_{1}=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D705}_{-} & \text{if}~0<\unicode[STIX]{x1D6FD}\leqslant \unicode[STIX]{x03C0}/2\\ 0 & \text{if}~\unicode[STIX]{x03C0}/2<\unicode[STIX]{x1D6FD}<\unicode[STIX]{x03C0},\end{array}\right.\quad \unicode[STIX]{x1D705}_{2}=\unicode[STIX]{x1D705}_{+}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D705}_{\pm }$ are given in (3.3), and the flow is spectrally unstable. When the flow is unstable, we can determine the nature of the instability by examining if it is possible to deform the contour of integration for $I_{c}$ so that it is entirely contained inside the region of the complex $\unicode[STIX]{x1D705}$ plane where $\text{Re}(\unicode[STIX]{x1D70E})<0$ . If this is possible, then the instability is convective and otherwise the instability is absolute. However, $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D705})$ is not analytic in the complex $\unicode[STIX]{x1D705}$ plane because of the $|\unicode[STIX]{x1D705}|^{3}$ term in (5.5) and so the classical approach of Huerre & Monkewitz (Reference Huerre and Monkewitz1990) cannot be applied directly. To proceed, we note that we can rewrite $I_{c}$ as (see Lin et al. Reference Lin, Pradas, Kalliadasis, Papageorgiou and Tseluiko2015; Vellingiri, Tseluiko & Kalliadasis Reference Vellingiri, Tseluiko and Kalliadasis2015) $I_{c}=2\,\text{Re}(K(x,t))$ , where $K(x,t)=\int _{0}^{\infty }\text{e}^{\unicode[STIX]{x1D70E}^{+}(\unicode[STIX]{x1D705})t}B(\unicode[STIX]{x1D705})\unicode[STIX]{x1D719}(x,\unicode[STIX]{x1D705})\,\text{d}\unicode[STIX]{x1D705}$ , where $\unicode[STIX]{x1D70E}^{+}(\unicode[STIX]{x1D705})$ is given by (5.5) with $|\unicode[STIX]{x1D705}|^{3}$ replaced by $\unicode[STIX]{x1D705}^{3}$ and is an analytic function in the entire complex $\unicode[STIX]{x1D705}$ plane. Furthermore, it is sufficient to consider the range of integration where $\text{Re}(\unicode[STIX]{x1D70E}^{+})>0$ and examine the integral

(5.7) $$\begin{eqnarray}\displaystyle K_{1}(x,t)=\int _{\unicode[STIX]{x1D705}_{1}}^{\unicode[STIX]{x1D705}_{2}}\text{e}^{\unicode[STIX]{x1D70E}^{+}(\unicode[STIX]{x1D705})t}B(\unicode[STIX]{x1D705})\unicode[STIX]{x1D719}(x,\unicode[STIX]{x1D705})\,\text{d}\unicode[STIX]{x1D705}. & & \displaystyle\end{eqnarray}$$

Since the integrand is analytic in the $\unicode[STIX]{x1D705}$ plane, we may deform the contour of integration into any path $\unicode[STIX]{x1D6E4}$ connecting $\unicode[STIX]{x1D705}_{1}$ and $\unicode[STIX]{x1D705}_{2}$ . The case of convective instability, in which the deformed contour $\unicode[STIX]{x1D6E4}$ lies entirely within a region with $\text{Re}\,\unicode[STIX]{x1D70E}^{+}(\unicode[STIX]{x1D705})<0$ , is shown in figure 15(a). The case of absolute instability, for which the contour $\unicode[STIX]{x1D6E4}$ must pass through a region with $\text{Re}\,\unicode[STIX]{x1D70E}^{+}(\unicode[STIX]{x1D705})>0$ , is shown in figure 15(c). The transition from convective to absolute instability under a continuous parameter change must happen by passing through the situation illustrated in figure 15(b) where the two lines on which $\text{Re}\,\unicode[STIX]{x1D70E}^{+}(\unicode[STIX]{x1D705})=0$ have pinched together at a saddle point, where $\text{d}\unicode[STIX]{x1D70E}^{+}/\text{d}\unicode[STIX]{x1D705}=0$ (note that since $\unicode[STIX]{x1D70E}^{+}$ is analytic, points where $\text{d}\unicode[STIX]{x1D70E}^{+}/\text{d}\unicode[STIX]{x1D705}$ vanishes are necessarily saddle points).

Figure 15.

Illustration of the change from (a) convective instability to (c) absolute instability. The solid lines show where $\text{Re}\,\unicode[STIX]{x1D70E}^{+}(\unicode[STIX]{x1D705})=0$ , and the broken line shows the integration contour $\unicode[STIX]{x1D6E4}$ . Regions with $\text{Re}\,\unicode[STIX]{x1D70E}^{+}(\unicode[STIX]{x1D705})>0$ and $\text{Re}\,\unicode[STIX]{x1D70E}^{+}(\unicode[STIX]{x1D705})<0$ are shown with a $+/-$ respectively. The black dot indicates a saddle point.

Following the above discussion, to identify the transition from convective to absolute instability, we seek the parameter values for which $\text{Re}\,\unicode[STIX]{x1D70E}^{+}(\unicode[STIX]{x1D705})=0$ and $\text{d}\unicode[STIX]{x1D70E}^{+}/\text{d}\unicode[STIX]{x1D705}=0$ at the relevant saddle point. To achieve this, we first take a pulse solution and solve the cubic equation $\text{d}\unicode[STIX]{x1D70E}^{+}/\text{d}\unicode[STIX]{x1D705}=0$ to find the relevant saddle point (by scrutinising the contours) in the $\unicode[STIX]{x1D705}$ plane. We then adjust $We$ until $\text{Re}\,\unicode[STIX]{x1D70E}^{+}(\unicode[STIX]{x1D705})=0$ is satisfied. We may then continue in any desired parameter to trace out the transition boundary. Figure 16 shows the boundary between absolute and convective instability in the $(Bo,We)$ plane for a range of values of the inclination angle $\unicode[STIX]{x1D6FD}$ . Evidently for fixed Bond number $Bo$ , the flow is either convectively unstable for any value of $We\geqslant 0$ , or else as $We$ increases it undergoes a transition from absolute to convective instability, and we, therefore, expect that a sufficiently strong electric field should have a regularising effect on the dynamics. From a physical point of view, this can be explained as follows. The electric field has a destabilising effect on the flow, and results in the generation of larger-amplitude pulses. The speed of the pulses is also amplified by the electric field. If the electric field becomes sufficiently strong, the pulses become sufficiently fast so that they can escape expanding wave packets generated by localised disturbances. As a result, the pulses become convectively unstable.

Our results are consistent with those found by Lin et al. (Reference Lin, Pradas, Kalliadasis, Papageorgiou and Tseluiko2015) for the non-local KS equation (3.2), which supports pulse solutions which are absolutely unstable for any $We$ (Lin et al. Reference Lin, Pradas, Kalliadasis, Papageorgiou and Tseluiko2015). This is in agreement with the findings shown in figure 16: for any Bond number there is a transition from absolute to convective instability as the Weber number is increased, and the smaller the Bond number the larger the Weber number at which the transition occurs.

Next, we turn our attention to a travelling pulse that is a solution to the full Stokes equations, as discussed in § 4. For a disturbance $\unicode[STIX]{x1D702}(x,t)$ , which is superimposed onto a pulse solution, as for the long-wave case, the solution to the linearised stability problem may be written in the form (5.4). For the present purposes, we will assume that the discrete spectrum does not contain any eigenvalues which lie in the right half-plane; consequently the first term in (5.4) decays to zero as $t\rightarrow \infty$ and, the convective/absolute nature of the instability is determined by study of the second term. As for the long-wave case the nature of the instability hinges on the ability to deform the contour of integration of the integral term $I_{c}$ , as discussed in figure 15.

All of the elevation boundary-element pulse solutions for Stokes flow computed in § 4, both electrified and non-electrified (see figures 7 c, 8 a and 9 a), are found to be convectively unstable. The negative pulse solutions presented in figure 14 are all absolutely unstable.

Figure 16.

Regions of convective ( $C$ ) and absolute ( $A$ ) instability in the $(Bo,We)$ plane for a range of inclination angles for the long-wave model (5.2): (a) angles $\unicode[STIX]{x1D6FD}=0.5\unicode[STIX]{x03C0},0.75\unicode[STIX]{x03C0},0.95\unicode[STIX]{x03C0}$ and (b) angles $\unicode[STIX]{x1D6FD}=0.05\unicode[STIX]{x03C0},0.25\unicode[STIX]{x03C0},0.5\unicode[STIX]{x03C0}$ .

Figure 17.

Time-dependent solutions of (3.1) for $\unicode[STIX]{x1D6FD}=0.75\unicode[STIX]{x03C0}$ , $Bo=0.002$ : (a) $We=0$ and (b)  $We=62$ . The top panels show the time evolution in a frame moving at speed (a) 1.4 and (b) 2.22, and the lower panels show the evolution of the norm $\Vert h-1\Vert$ .

Figure 18.

Time-dependent solutions of (3.1) for $\unicode[STIX]{x1D6FD}=0.25\unicode[STIX]{x03C0}$ , $Bo=0.01$ : (a)  $We=12.5$ , (b)  $We=14.5$ , (c)  $We=26$ . The top panels show the time evolution of the solution in a frame moving at speed (a,b) 1.4 and (c) 2.08. The lower panels show the evolution of the norm $\Vert h-1\Vert$ .

6 Physical context

It should be possible to observe the phenomena described in this study in a physical experiment. In a typical experimental set-up, a liquid film emerges from an inlet and flows down an inclined plane. We would expect to observe two-dimensional wave phenomena for not too large Reynolds number and sufficiently close to the inlet. However, three-dimensional behaviour is likely further away from the inlet, in particular for hanging films as was observed for example by Charogiannis & Markides (Reference Charogiannis and Markides2016).

Previous experimental studies have used films down to just a few microns in thickness. Assuming an aqueous film of thickness $h_{0}=100~\unicode[STIX]{x03BC}\text{m}$ at $20\,^{\circ }\text{C}$ , we find that $Bo=7.0\times 10^{-4}$ and $Re=4.9$ . Referring to figure 16(a), the threshold between convective and absolute instability is at $We=15.65$ and so in order to observe the regularisation of the dynamics discussed in § 5.1 we might take $We=20$ , say, for an experiment with $\cot \unicode[STIX]{x1D6FD}=0.95\unicode[STIX]{x03C0}$ (an inverted film). This requires an electric field of strength $E_{0}=1.487\times 10^{6}~\text{V}~\text{m}^{-1}$ . Assuming an electrode potential of 1000 V, a gap between the film and the electrode of approximately $672~\unicode[STIX]{x03BC}\text{m}$ would be required, and this is almost seven times the film height (our analysis assumes that the electrode is a long way from the film surface). For the same parameter values, but considering a slightly inclined plate with $\cot \unicode[STIX]{x1D6FD}=0.05\unicode[STIX]{x03C0}$ , the transition between convective and absolute instability occurs at $We=71.45$ (see figure 16 b). Taking $We=75$ would require an electric-field strength $E_{0}=2.880\times 10^{6}~\text{V}~\text{m}^{-1}$ , which is also below the dielectric breakdown limit (in this case with the same electrode potential the electrode–film gap is approximately 3.5 times the film height).