Nonlinear gravity electro-capillary waves in two-fluid systems: solitary and periodic waves and their stability

Starting from the Euler equations governing the flow of two immiscible incompressible fluids in a horizontal channel, allowing gravity and surface tension, and imposing an electric field across the channel, a nonlinear long-wave analysis is used to derive a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\times 2$$\end{document}2×2 system of evolution equations describing the interface position and a modified tangential velocity jump across it. Travelling waves of permanent form are shown to exist and are constructed in the periodic case producing wave trains and the infinite case yielding novel gravity electro-capillary solitary waves. Various regimes are analysed including a hydrodynamically passive but electrically active upper layer, pairs of perfect dielectric fluids and a perfectly conducting lower fluid. In all cases, the presence of the field produces both depression and elevation waves travelling at the same speed, for given sets of parameters. The stability of the non-uniform travelling waves is investigated by numerically solving appropriate linearised eigenvalue problems. It is found that depression waves are neutrally stable whereas elevation ones are unstable unless the surface tension is large. Stability or instability is shown to be linked mathematically to the type of local eigenvalues of the nonlinear flux matrix used to obtain travelling and solitary waves; if these are real (hyperbolic flux matrix), the system is stable, and if they are complex (elliptic), the system is unstable. The latter is a manifestation of Kelvin–Helmholtz instability in electrified flows.

by y = S(x, t). Taking the fluids to be perfect dielectrics, we introduce an electric potential V (x, y, t) within the channel by imposing constant voltages V (x, −h) = 0 and V (x, h) = V 0 at the lower and upper wall electrodes, respectively. In the absence of an electric field, the problem was studied by [16], and in fact, this is equivalent to having two fluids of equal permittivity. The superscripts + and − are used to denote variables in the upper and lower fluids, respectively, and the permittivities are denoted by 0 ± where 0 is the permittivity of free space-hence, ± are dimensionless and in what follows the ratio = + / − will become important. We also note that this setup was introduced and studied by Melcher [5] in his seminal monograph that established the field of interfacial electrohydrodynamics (see also Melcher and Taylor, Saville reviews). Our study extends that of [5] by introducing a physical dispersive regularisation by the inclusion of surface tension-in the absence of a regularisation, the 2 × 2 system of nonlinear model equations not only is at best hyperbolic, and hence terminates in shocks, but can also be elliptic (with catastrophic short-wave instabilities) if the electric field is sufficiently strong, a case that was not studied in [5]. Here, we present a complete study of the regularised system and in particular show the existence of solitary waves and construct families of them.
The fluids are irrotational and incompressible; hence, the velocity field can be written as u ± = ∇φ ± where φ ± denote the velocity potential in each region. Incompressibility then implies that φ ± are harmonic. The electrostatics is also governed by Laplace equations for V ± . To see this, start with the electrostatic limit of the Maxwell equations so that ∇ × E ± = 0, where E is the electric field, introduce voltage potentials so that E ± = −∇V ± , and use Gauss's law ∇ · ( 0 ± E ± ) = 0 in each region to obtain the result (the permittivities are constant). The boundary conditions at the walls are fluid impermeability and fixed voltage potentials. Those at the interface y = S(x, t) are the kinematic conditions, continuity of the voltage potential, continuity of the electric displacement deriving from Gauss's law, and continuity of normal stresses (tangential stresses cannot be prescribed, since the fluids are inviscid-slip is allowed at the interface). Hence, the governing equations and boundary conditions are

1d)
V + − V − = 0 at y = S, (2.1e) + n · ∇V + − − n · ∇V − = 0 at y = S, (2.1f) n T + n − n T − n = −σ κ at y = S. (2.1g) In (2.1f), (2.1g), is the transpose, n = (−S x , 1)/ 1 + S 2 x is the unit normal pointing into fluid +, and T is the combined fluid and Maxwell stress tensor which reads The parameter σ is the constant surface tension between the fluids, and κ = S x x /(1 + S 2 x ) 3/2 is the signed curvature of the interface. Imposition of initial conditions completes the mathematical statement of the problem.
We proceed to non-dimensionalise the problem taking a typical interfacial wavelength to be L, and a typical dimensional speed to be c. The interface is scaled by the channel height h and the pressure by its inertial scale. We write 3) where starred variables are dimensionless. It is convenient to introduce the following dimensionless parameters that arise from the scalings (2.3): These can be identified as a slenderness parameter β, a density ratio ρ, a Froude number F measuring the importance of gravity, a scaled inverse Weber number σ * that we retain in order to keep surface tension, and an electric Weber number measuring the ratio of electrostatic to inertial pressures. The scaled inverse Weber number σ * is taken to be an O 1 β parameter to allow the surface tension to compete with gravity and the electric field. In what follows we also use the Atwood ratio: Dropping the * from our dimensionless variables, we write the governing equations as follows: with the boundary conditions for ± at the interface taking the form: Equation (2.7b) is the Bernoulli equation at the interface and is derived in the familiar way of integrating the Euler equations to obtain (2.8) and evaluating this equation at each side of the interface-see [8] for its use in electrohydrodynamic problems.
The voltage potential boundary conditions at Y = S(X, T ) become while the conditions at the walls yield Using the expression (2.2) in the normal stress balance (2.1g) and non-dimensionalising as described above, yields the following expression for the pressure jump across the interface 11) to be used in the Bernoulli equation (2.7b), thus, eliminating the pressure. Our system of equations at this point is still exact, and a numerical procedure based on the full equations is necessary. In what follows we consider the long wave limit of β 1 and make progress asymptotically to derive reduced-dimension nonlinear evolution equations for the flow.

Asymptotic expansions and derivation of the evolution equations
In the limit β → 0, we seek a solution of the problem by the asymptotic expansions These expansions are valid as long as gradients, e.g. S 0x , remain bounded, and this is confirmed a posteriori by solving the resulting equations. Substituting expansions (3.1) into the dimensionless system (2.6), we obtain to leading order which can be readily integrated to give where the functions A ± 1,2 and B ± 1,2 must be determined. Using the no penetration conditions (2.10a) at the walls immediately implies that A ± 1 ≡ 0, and hence, ± 0 ≡ ± 0 (X, T ), i.e. they are independent of Y . Equation (2.6a) at order β 2 gives which when integrated once and on use of the no penetration wall conditions yields These expressions will be useful in the first-order kinematic conditions considered below. The wall voltage conditions (2.10b) yield, in turn, The unknown functions in (3.6) can be found from the leading order contributions to the voltage and displacement field continuity at the interface equations: .
In order to obtain equations for the remaining unknowns S 0 and ± 0 , we use an approach similar to [16]. From the O(β 2 ) kinematic conditions (2.7a) on either side of the interface, and using the solutions (3.5), we find We now define new variables, U and W , representing the average horizontal velocity and half its jump across the interface, Subtracting either of the equations in (3.8) from the other and integrating with respect to X gives where χ(T ) is an arbitrary function of T . Adding equations (3.8a) and (3.8b) and using (3.10) to eliminate U yield A second-evolution equation follows from the Bernoulli equation (2.7b) after its differentiation with respect to X and use of (3.9) and (3.10). The result is (recall that α is the Atwood number defined in (2.5)) Defining a new variable allows us to write (3.13) as the following dispersively modified 2 × 2 system of conservation laws (3.15) where the matrix entries are (3.16d) As we will see below, the electric field has a considerable effect on the flow by introducing complex eigenvalues for the matrix Q. In the absence of an electric field, most easily recovered by setting = 1 in (3.15), we recover the equations of [16] where the eigenvalues of the analogous matrix Q are now real, and hence, in the absence of dispersion, the system is hyperbolic. It is instructive to consider the linear limit of this problem for expanding S, like where ζ 1. The system (3.15) can be reduced to a single equation forS 0 Looking for a wave-like solution of the formS 0 = e ik X+ωtŜ (X, T ) gives the dispersion relation we have growth, and if the problem is dispersive, so under no circumstances is the linear system dissipative.
The full problem stated in (3.15) is a nonlinear free boundary one with two fluid phases and a large number of physical parameters. The asymptotic reduction supports analytical progress and in particular allows for the exploration of nonlinear coherent structures such as travelling wave trains and solitary waves. We explore such possibilities next.

Solitary waves at arbitrary density ratios
We will construct solitary wave solutions of the system (3.13) and find the range of parameters for which they exist. Looking for travelling wave solutions of (3.13) of the form where c is the positive wavespeed, yields the coupled ordinary differential equations where c is the wave-speed. It is trivial to solve the first equation for W in terms of S, thereby eliminating W from the problem; we find where D is a constant of integration. In order to integrate (4.2b), we will need where C 1 is a constant-this follows by use of the expressions (4.3). Next, we integrate (4.2b) with respect to ξ , multiply the result by S ξ , carry out an additional ξ integration using (4.4), and yield the following dynamical system for S: , (4.5) where G and H are constants and the parameters measure surface tension and electric field strength, noting that γ > 0 and δ ≥ 0. By requiring that as |ξ | → ∞, S and all its derivatives tend to 0 (thus restricting our solutions to solitary waves), we see from (4.3) that D is the jump in horizontal velocity across the interface, i.e. it is the undisturbed vortex sheet strength. Equation (4.5) becomes, after appropriate selection of G and H to achieve the desired decay at infinity, Note that for solitary waves, γ can be scaled to unity by a redefinition of ξ . We can see that this equation implies that α > 0 is a necessary condition for solitary waves to exist, hence excluding Rayleigh-Taylor unstable configurations as would be expected. Note that this condition is the same as in the non-electrified case considered by [16], consistent from the fact that a vertical electric field is destabilising. In order to analyse the properties of the solitary waves that may be present in the channel, it is useful to rewrite Eq. (4.7) in the form where p 3 (S) is a cubic polynomial given by The constants j, k, l and m are given by Inspection of (4.8) and noting the facts |S| < 1 and + 1 + S( − 1) > 0, implies that the sign of p 3 determines the existence of solitary waves. First, we must have p 3 (S) > 0 for small |S| so that S 2 ξ is positive. By continuity, we require p 3 (0) = m > 0 which gives the following necessary condition for the existence of solitary waves: Solitary waves emerge when we can connect S = 0 to another root(s) S r , say, as long as |S r | < 1. For this to happen, we need at least one real root of p 3 (S) to lie in −1 < S < 1. Now, with p 3 (1) = 0 when α = 1, i.e. ρ = 0. This is possible when the upper fluid is passive -this is considered in more detail later-see also [16]. Further special cases arise when D = ±c, in which case p 3 has roots at either S = −1 or S = 1, and in both cases, it is possible to have zero, one, or two solitary waves produced depending on the other parameters-this has been confirmed numerically, but details are not included for brevity. As these cases do not correspond to a meaningful physical limit since the interface touches the wall, they have not been investigated further. From (4.12), we conclude that if it is assumed that p 3 (S) has real roots in −1 < S < 1 (something that is necessary for solitary waves to exist), then it must have exactly two roots or a repeated double root; otherwise, it is impossible to satisfy all the inequalities of (4.11), (4.12). The latter scenario is inadmissible because it forces S 2 ξ < 0 in the vicinity of S = 0 even if the repeated root is at S = 0. The former scenario is again inadmissible unless the roots have opposite signs, and this is the generic situation for the existence of solitary waves in this rather general case. Interestingly, we establish that solitary waves come in pairs, an elevation wave corresponding to 1 > S + r > 0 and a depression (or dark) solitary wave for the root −1 < S − r < 0. An illustration of this is provided in Fig. 2, in which regions where (S ξ ) 2 < 0 are shown for mathematical completeness despite this not being physical. Integration of (4.8) provides integrals for the solutions analogous to (4.19) constructed below.

Solitary waves with upper fluid of zero density (α = 1)
When the upper fluid has negligible density compared to the lower fluid the Atwood ratio α = 1. The equations simplify in this case, and in particular are not prone to a Kelvin-Helmholtz instability since the system is one-sided for the hydrodynamics. We begin with a general dielectric fluid in the lower region and then consider the additional limit of a perfectly conducting lower liquid first studied by Melcher [5] (Chapter 7).

Perfect dielectric lower fluid
In order for this model to be physical we will assume that the ratio of permittivities is < 1. In this special case, S = 1 is a root of p 3 so we can simplify (4.8) to where (4.14) Evaluating p 2 at S = 0 and requiring p 2 (0) > 0, it can be seen that in this case, the necessary condition for solitary waves to be produced is As in the previous section, we have that p 2 (−1) = −(D − c) 2 < 0, so for the existence of solitary waves, we require one root of p 2 is −1 < S < 0. However, the right-hand side of (4.13) is now finite as S → 1 and can take either a positive or negative value. Assuming that there is one root of p 2 in the range −1 < S < 0, if the second root is in 0 < S < 1, then there will be a second solitary wave produced. This extra condition is equivalent to requiring is positive we have only one depression solitary wave, if it is negative, we have an additional elevation solitary wave. An example of how the value of p 2 (1) depends on is given in Fig. 3b, with = 0.15 0.2 giving p 2 (1) < 1 so producing two solitary waves. From (4.15) and the expression for p 2 (1), we have the explicit condition for the existence of two solitary waves: For one solitary wave, we must have The dependence on the number of roots on the ratio of permitivitties can be seen in Fig. 3a for a typical set of parameter values given in the caption. These can be found explicitly from (4.14), and they are plotted when they are real as functions of . In region 1, there are no admissible solitary waves since the roots are both negative. For > 0.07 we have one positive and one negative root, admitting two solitary waves. Region 2 terminates at ≈ 0.23 since the positive root now equals unity, i.e. touches the wall and must be excluded. Beyond this value we have region 3 which supports a single depression solitary wave. The region 3 solutions are consistent with the results of [16] who found only depression solitary waves in the case of no electric field equivalent to the = 1 case here. Example phase planes for certain values of are given in Fig. 3b.
We can also numerically solve Eq. (4.13) for S and construct solitary waves. In the present case, we have either one or two roots of p 2 (S) in the physical range −1 < S < 1. As discussed above, we require exactly one root, denoted by S − , to be in the range −1 < S − < 0. Denoting the other root of p 2 (S) to be S 0 (note that we do not assume 0 < S 0 < 1), we can write (4.13) in the form . (4.18) Separating variables and integrating from the trough S = S − to an elevation height S yields The positive square root was taken in (4.19) to produce half the wave-the negative root gives the other symmetric half. Numerically solving this equation for varying parameters produces the lower (red) solitary wave depicted in Fig. 4. A depression solitary wave is always produced if the necessary condition (4.15) for solitary waves is met. If in addition (4.16) is met so that 0 < S 0 < 1, then a second elevation solitary wave is supported having the same speed as its depression counterpart. In the example of Fig. 4 we also show the elevation wave in blue.  (4.20) where in this case, we have that for solitary waves, If the eigenvalues of the nonlinear flux matrix in (4.20) are denoted by μ(S, ), then we find that they are determined through the equation: We find that complex eigenvalues of the flux matrix in (4.20) arise when When we have μ is complex for a range of S we call such a region elliptic, and when μ ∈ Re, we refer to the region as hyperbolic, so S tr is the transition point where the system changes from hyperbolic to elliptic or vice versa. An immediate consequence of (4.23) is that when is near 1 we can make S tr > 1, i.e. elevation waves can be found which are what we term hyperbolic; this is in contrast to the = 0 case discussed in Sect. 4.2.2. We also note that the bound (4.23) only makes physical sense when the parameters are such that additionally |S| < 1. In the case of two solitary waves, we have from the inequality: (4.16) Using (4.24) in (4.23) provides the following range of values for the transition boundary S = S tr : This inequality in turn predicts that the solitary waves produced when α = 1, 0 < < 1, can potentially have parts where the equations are locally elliptic and other parts where they are locally hyperbolic, in addition to wholly elliptic or hyperbolic. It is quite easy to check numerically the upper and lower bounds of S tr in (4.25) for 0 < < 1, and to confirm that transitions occur for both depression and elevation waves when 0 < 0.587. Such analytical estimates were used prior to searching for transitional solitary waves like those in Fig. 6.
We note that the classification of solitary waves into elliptic and hyperbolic regions presented here does not take into account the dispersive term on the right-hand side of (4.20) so strictly speaking only concerns the conservation laws when σ = 0. However, this diagnostic tool yields useful information because a change from real to complex eigenvalues of the conservation laws predicts the presence of instabilities that destroy the solitary wave structures, which enables us to predict instability of nonlinear solitary waves without doing any spectral analysis. Analogous analyses and classifications have been used in related viscous multifluid interfacial problems where the regularising terms are diffusive -see for example [17,18].

Lower fluid a perfect conductor
As mentioned earlier this limit was considered in [5] in the absence of surface tension, and hence the solutions there cannot produce solitary waves. The perfect conductor lower fluid limit is found by sending its permittivity to infinity, − → ∞. Hence, we have = + / − = 0 and Eq. (4.13) reduces to Using similar reasoning as before, the necessary condition for solitary waves for this sub-case is p * and consequently solitary waves can only exist if F < 1/δ, i.e. for sufficiently small Froude numbers, given an electric field strength. We also note that p * 2 (−1) = −(D − c) 2 < 0, p * 2 (1) = −2δ < 0; hence, a solitary wave of depression and one of elevation coexist in general. Typical pairs of solitary waves are given in Fig. 5 as F varies.
For completeness, we consider the mathematical structure of the waves constructed in this section. The flux matrix in this case now reads It can be shown that S a defined above is always smaller than the amplitude of elevation solitary waves given by the positive root of (4.27). From (4.28), we have δ F < 1; combining this with (4.30) and the physical fact that δ and F are non-negative, we arrive at the following general condition that is necessary for the evolution equations to become locally elliptic (we exclude the boundaries where (4.22) gives two equal and real eigenvalues) The inequality (4.30) in turn predicts that any solitary wave satisfying 1 > S > 1 − 2 1/3 ≈ −0.26, can potentially have parts where the equations are locally elliptic and other parts where they are locally hyperbolic, in addition to wholly elliptic or hyperbolic. These properties hold for depression as well as elevation waves, and examples of such mixed behaviour are given in Fig. 6. Figure 6a has Froude number F = 0.65, while Fig. 6b has F = 0.1, the other parameters being δ = 1, σ = 0.1 (note that α = 1 and = 0 here). For F = 0.65, the depression wave exhibits mixed behaviour as seen in the figure that depicts the elliptic parts in red and the hyperbolic ones in blue. When F = 0.1, the elevation wave now supports ellipticity where its amplitude is sufficiently large, as seen in Fig. 6b. We note that if transition takes place in the elevation wave then the depression wave is wholly hyperbolic, and if it takes place in the depression wave then the positive wave is wholly elliptic. This can be inferred from the monotonicity of the sufficient condition for ellipticity 1 > S > 1 − 2 1/3 ≈ −0.26. Of course, exact diagnostics are calculated directly from the eigenvalues (4.22).

Periodic travelling waves
Section 4 is concerned with solitary waves that decay far away. Here, we construct travelling waves of finite spatial period governed by equations (3.15). We will concentrate on the physical regime of Sect. 4.2.2, i.e. we take α = 1 and − = ∞ (i.e. = 0) which restricts our analysis to the case of a hydrodynamically passive upper dielectric region and a perfectly conducting lower fluid (the general case can be analysed in analogous ways and is excluded for brevity). Looking for travelling waves in a frame ξ = X − cT as before, casts the system (3.15) into where 2L is the wavelength, and note that due to the assumed symmetry, we only need to solve the system on the half domain 0 ≤ ξ ≤ L.
For a given set of parameters F, σ, δ, to solve (5.1), we require 4 boundary conditions, two of which we impose on S and two on , both at the points ξ = 0, L. Integrating each of (5.1a) and (5.1b) between 0 and L yields the conditions which can be thought of as two equations for the unknowns c and L for given physical parameters F, σ, δ, and given wave and velocity amplitudes, i.e. the four values S(0), S(L), (0), (L). Note that not all sets of parameters and end conditions give rise to physically admissible solutions. The integrals (5.2a)-(5.2b) allow elimination of c and derivation of a nonlinear relation between the end conditions which we conveniently state as the following quadratic equation for (0): For numerical purposes, it is simpler to specify the spatial period L and think of (0) as the unknown quantity to be found as part of the nonlinear two-point boundary value problem. The present system allows for an exact evaluation of (0) from (5.3) in contrast to iterations that would be typically required. For real admissible solutions, we require that a 2 2 − 4a 1 a 3 > 0 where equation (  speeds c 1 > c 2 for a given amplitude. An example when the two waves have equal and opposite speeds is given in Fig. 7 for parameter values F = 0.1, δ = σ = 0.01. The two branches in this case are symmetric due to the fact that the choice (L) = 0 has been made, which implies from (5.3) that the (real) roots for (0) have equal and opposite signs. It follows from expression (5.2a) that c has two equal and opposite values as seen in Fig. 7, where in addition we have fixed the length to L = 1/2 and also took S(L) = 0. Note that the choices S(L) = (L) = 0 are natural for the recovery of solitary waves from periodic ones as we see later. The graph in Fig. 7 has been extended to the y-axis by using the linear dispersion relation: if |S(0) − S(L)| 0, and to the x-axis by noting that as c → 0 we have S(0) → −1, i.e the wave touches the bottom channel wall for waves with zero speed.
While it is necessary to solve the full nonlinear system (5.1) numerically, it is useful to first obtain the linear solutions analytically and use them in numerical continuation to larger amplitudes. We look for linear travelling wave solutions in the form: where the amplitude |S(0) − S(L)| 1. Linearising the equations and looking for solutions S, proportional to cos(π ξ/L) yields (note that having fixed the end conditions, the period 2L is an eigenvalue to be determined): where the frequency is given by In the results that follow, we fix S(L) = (L) = 0 as discussed earlier and then construct periodic waves numerically for given L and amplitude S(0). Typical results for parameter values F = 0.01, δ = 1 and σ = 1, with The smallest amplitude waves are in agreement with the linear results described above, but as the amplitude increases, we observe convergence to solitary waves. Such convergence is corroborated further in Fig. 9 for a slightly different set of parameters F = δ = σ = 0.1 that superimposes the depression solitary wave found using the analysis of Sect. 4. Agreement is seen to be excellent.

Stability of periodic and solitary travelling waves
In this section, we consider a linear stability analysis of both the solitary wave solutions derived in Sect. 4 as well as the periodic travelling waves of Sect. 5. We shall perform this analysis in the case of a hydrodynamically passive upper fluid (i.e. α = 0) and the lower fluid a perfect conductor (i.e = 0) -the analysis can easily be extended to arbitrary parameters. The non-uniform travelling wave solutions in this case are governed by the coupled equations (5.1) with α = 1, = 0, and we denote them by S(ξ ) and (ξ ), recalling that they travel with speed c and ξ = X − cT . We perturb these solutions so that where |S (ξ, T )|, | (ξ, T )| 1, substitute into the governing equations (3.15) and linearise to find We assume that the perturbed solutions are of the form S = e λT e ikξ S(ξ ), = e λT e ikξ (ξ ), (6.3) where S(ξ ) and (ξ ) are 2L periodic and the wavenumber 0 ≤ k ≤ π/L allows perturbations that are longer than 2L, e.g. subharmonic. In what follows we take k = 0, i.e. we consider waves with the same wavelength as the basic period, and use these results with comparisons with corresponding solitary wave stability later. Substituting (6.3) into (6.2a)-(6.2b) (with k = 0), yields the following eigenvalue problem to determine λ: The eigenvalue problem (6.4) is solved numerically after specifying periodic boundary conditions. Finite difference methods are used, and the system is cast into a matrix eigenvalue problem of the form Ar = λr where r is the discretisation of S, on the grid. We begin with the stability of the depression waves shown in Fig. 9, recalling that we have superimposed a calculated solitary wave along with its large wavelength analogue. We have carried out the stability by using null conditions far away for the solitary wave profiles, as well as spatially periodic boundary conditions for large enough periods. The resulting spectra are the same and shown in Fig. 10. It can be seen that all eigenvalues lie on the imaginary axis, and hence, the waves are neutrally stable. The vertical extent of the spectrum increases as resolution increases and higher wave numbers enter into the calculation-this is expected due to the dispersive regularisation provided by surface tension. We also note that according to the criterion (4.30), the waves in Fig. 10 are what we termed hyperbolic, i.e. the eigenvalues of the accompanying nonlinear flux matrix are real, and this provides an explanation for the stability found in our computations. In fact, all wholly hyperbolic depression solitary waves studied were found to be neutrally stable, and hence, criterion (4.30) can be used as a simple rule of thumb to determine linear stability (albeit to waves with wavelengths at most 2L-modulational stability will be considered elsewhere). The converse of this rule of thumb also holds: namely, if the wave is locally elliptic, then the spectrum contains at least one eigenvalue with positive real part and the system is unstable. As all elevation solitary waves are proven to be at least locally elliptic, these are found to be unstable. The effect of increasing the amplitude of periodic depression waves of permanent form was investigated in Fig. 8 for F = 0.01, δ = 1 and σ = 1. For the stability of those waves, in Fig. 11 we give results for the smallest amplitude wave having S min = −0.1, and largest amplitude one having S min = −0.8. The results show that the smaller wave is neutrally stable, whereas the larger one has one unstable mode with real part ≈ 0.259. We note that the number of unstable modes increases as the surface tension σ is decreased (not shown) due to the reduced dispersive regularisation when surface tension is weaker.
Turning to elevation waves next, in general, we find that these are typically more unstable than depression waves for a given set of parameters and amplitude. The instability behaviour is analogous to that of depression waves, in the sense that longer waves appear to be more stable that shorter ones having the same amplitude, and larger amplitude waves are also more unstable. Typical results are included in Fig. 12 for parameters chosen to give a long waves of period 14 units (essentially solitary), and a shorter wave of period 2 units, both of relatively large amplitude equal to 0.4. Panels (a) and (b) show the solitary wave and its spectrum, with panels (c) and (d) depicting the finite-period wave and its instability characteristics. Both of these waves are wholly elliptic according to the criterion (4.30), and hence, it is not surprising to observe unstable eigenvalues. We note that as the surface tension parameter σ is increased, the waves become more sinusoidal and also linearly stable (this was confirmed but not shown for brevity).
Finally, we consider the stability of the non-electrified solitary waves that were reported in [16]. These were constructed for different parameters but their stability was not investigated. Our computations indicate that all the solitary waves in [16] are neutrally stable, and we illustrate this for the case α = 1 (upper layer is hydrodynamically passive). To switch off the field in our model, it is sufficient to take = 1 and δ = 0, construct solitary waves as in Sect. 4.2.1, and study their stability according to the eigenvalue problem (6.4). As noted in [16], depression solitary waves alone exist, and calculation of the eigenvalues of the flux matrix gives real values only, hence, the waves are what we termed of hyperbolic type. Their stability is analogous to the = 0 case considered in Sect. 4.2.2. The profile for one of the exact solutions found in [16] and its spectrum, showing neutral stability, are given in Fig. 13 for the set of parameters F = c = σ = 1. The exact explicit solution used here is given by Eqs. (54)-(55) in [16].

Conclusions
In this paper, we considered a two-fluid flow in a horizontal channel and allowing for the effects of gravity, electrical fields, and surface tension. We derived a 2 × 2 system of nonlinear 1-D evolution equations by carrying out a fully nonlinear long-wave asymptotic analysis that describes wave amplitudes scaling with the channel thickness. The conditions under which solitary wave solutions exist were derived in the case of a lighter upper fluid (i.e. a Rayleigh-Taylor stable regime), and the role of the electric field was elucidated. In the case when the upper fluid has negligible density (i.e. is hydrodynamically passive), and the lower fluid is a perfect conductor, it was shown that two solitary waves are always produced, an elevation and a depression wave, having the same speeds. Using the eigenvalues of the nonlinear flux matrix of the system, all travelling waves were classified to be locally elliptic or hyperbolic depending on whether the local eigenvalues are complex or real, respectively. The zero surface tension limit recovers the equations of [5], and the elliptic/hyperbolic diagnosis was extended to more general dielectric fluids where − = ∞. Periodic travelling waves were also constructed and the relationship between amplitude and wave speed was investigated, prior to carrying out a linear stability analysis of both classes of non-uniform solutions. In general, depression solitary and periodic waves were found to be stable whereas elevation ones can be unstable if the surface tension parameter is not too large. Our results show, as expected that the presence of locally elliptic regions in the travelling waves produces instability, and a necessary condition for neutral stability is that the profiles are locally hyperbolic over the whole spatial period. Future directions include the investigation of disturbances that are longer than the basic period of travelling waves (modulational instabilities), and a more complete study of shock type solutions involving heteroclinic connections and their stability.