2.1. Governing equations

Two immiscible inviscid incompressible fluids are bound together in an infinite horizontal channel and separated by a sharp interface $y=\eta (x,t)$, where $x$ is the direction of wave propagation, and the $y$-axis points upwards, with $y=0$ at the undisturbed interface (see the sketch of the system in figure 1). We denote by $h^\pm$ and $\rho ^\pm$ the depth and density in each layer, where the superscripts ‘$+$’ and ‘$-$’ refer to fluid properties associated with the upper and lower fluid layers, respectively. The fluids, which are assumed to be perfect dielectrics with electrical permittivities $\varepsilon ^\pm$, are under the action of a uniform horizontal electric field of strength $E_0$ in the absence of perturbations. The flows are supposed to be irrotational, hence there exist potential functions $\phi ^\pm$ such that the velocity fields $\boldsymbol {u}^\pm$ satisfy $\boldsymbol {u}^\pm =\boldsymbol {\nabla } \phi ^\pm$, where $\boldsymbol {\nabla } =(\partial _x,\partial _y)$ is the gradient operator. If we denote the electric fields by $E^+$ and $E^-$ in the corresponding layers, then the electrostatic limit of Maxwell's equations yields $\boldsymbol {\nabla } \times E^\pm =0$; therefore voltage potentials $V^\pm$ can be introduced, such that $E^\pm =-\boldsymbol {\nabla } V^\pm$. In this scenario, both velocity potentials and voltage potentials satisfy the Laplace equation:

(2.1a,b)\begin{align} \varDelta \phi^+ = 0, \quad \varDelta V^+ = 0, \quad \eta< y< h^+, \end{align}

(2.2a,b)\begin{align} \varDelta \phi^- = 0, \quad \varDelta V^- = 0, \quad -h^-{<}y<\eta, \end{align}

where $\varDelta =\partial _{xx}+\partial _{yy}$ is the two-dimensional Laplace operator. The electric fields have to satisfy two boundary conditions at the interface $y=\eta (x,t)$, namely the continuity of electric potential and continuity of the electric displacement field given by

(2.3a,b)\begin{equation} V^+ = V^-, \quad\varepsilon^+\,\frac{\partial V^+}{\partial\boldsymbol{n}}=\varepsilon^-\,\frac{\partial V^-}{\partial\boldsymbol{n}}, \end{equation}

where $\boldsymbol {n}=(-\eta _x, 1)/\sqrt {1+\eta _x^2}$ is the unit normal vector pointing outwards from the lower layer. Far away from interfacial disturbances, we should impose the condition

(2.4)\begin{equation} V^\pm \rightarrow E_0x, \quad\text{as $x\rightarrow\pm\infty$}. \end{equation}

Following Barannyk et al. (Reference Barannyk, Papageorgiou and Petropoulos2012, Reference Barannyk, Papageorgiou, Petropoulos and Vanden-Broeck2015), the impermeability conditions should be satisfied for both fluids and electric fields on the channel walls, namely

(2.5)\begin{equation} \left. \begin{array}{c@{}} \dfrac{\partial\phi^+}{\partial y}=\dfrac{\partial V^+}{\partial y}=0, \quad\text{at $y=h^+$}, \\ \dfrac{\partial\phi^-}{\partial y}=\dfrac{\partial V^-}{\partial y}=0, \quad\text{at $y=h^-$}. \end{array}\right\}\end{equation}

On physical grounds, the no-current boundary condition for electric fields is used to model electrically insulating walls. The hydrodynamic boundary conditions at the interface are a kinematic condition

(2.6)\begin{equation} \eta_t=\phi^+_y-\eta_x\phi^+_x=\phi^-_y-\eta_x\phi^-_x, \end{equation}

which indicates continuity of normal velocity, and a modified dynamic boundary condition resulting from pressure continuity across the interface:

(2.7)\begin{align} &\rho^-\left(\phi^-_t+\tfrac12|\boldsymbol{\nabla} \phi^-|^2+g\eta\right)-\rho^+\left(\phi^+_t+\tfrac12|\boldsymbol{\nabla} \phi^+|^2+g\eta\right)\nonumber\\ &\quad +[\boldsymbol{n}\boldsymbol{\cdot} \varSigma\boldsymbol{\cdot} \boldsymbol{n}]^-_+{-}\frac{\sigma\eta_{xx}}{(1 +\eta_x^2)^{3/2}}=0, \end{align}

where $g$ is the acceleration due to gravity, $\sigma$ is the coefficient of surface tension, $\boldsymbol {n}\boldsymbol {\cdot } \varSigma \boldsymbol {\cdot } \boldsymbol {n}$ arises from the Maxwell stress tensor $(\varSigma _{ij})_{i,j=1,2}$ given by

(2.8a–c)\begin{equation} \varSigma_{11}=\frac{\varepsilon}{2}(V_x^2-V_y^2), \quad \varSigma_{12}=\varSigma_{21}=\varepsilon V_xV_y, \quad \varSigma_{22}=\frac{\varepsilon}{2}(V_y^2-V_x^2), \end{equation}

and $[\cdot ]^-_+$ represents the jump in this quantity across the interface. A straightforward calculation yields

(2.9)\begin{equation} [\boldsymbol{n}\boldsymbol{\cdot} \varSigma\boldsymbol{\cdot} \boldsymbol{n}]_+^- = \varepsilon^-\left[\left(\frac{\partial V^-}{\partial\boldsymbol{n}}\right)^2-\tfrac12|\boldsymbol{\nabla} V^-|^2\right]-\varepsilon^+\left[\left(\frac{\partial V^+}{\partial\boldsymbol{n}}\right)^2-\tfrac12|\boldsymbol{\nabla} V^+|^2\right]. \end{equation}

For convenience, we introduce new functions $W^\pm =V^\pm /E_0-x$, hence $W^\pm \rightarrow 0$ as $x\rightarrow \pm \infty$, and $W^\pm$ still satisfy Laplace's equation. Therefore, (2.3a,b) and (2.9) can be rewritten as

(2.10)\begin{gather} (\varepsilon^+{-}\varepsilon^-)\eta_x=\varepsilon^+ (W^+_y-\eta_xW^+_x)-\varepsilon^-(W^-_y-\eta_xW^-_x), \end{gather}

(2.11)\begin{align} [\boldsymbol{n}\boldsymbol{\cdot} \varSigma\boldsymbol{\cdot} \boldsymbol{n}]_+^-&= E_0^2\varepsilon^-\left[\left(\frac{\partial W^-}{\partial\boldsymbol{n}}\right)^2-\tfrac12|\boldsymbol{\nabla} W^-|^2\right] -E_0^2(\varepsilon^-W_x^--\varepsilon^+W_x^+)\nonumber\\ &\quad -E_0^2\varepsilon^+\left[\left(\frac{\partial W^+}{\partial \boldsymbol{n}}\right)^2-\tfrac12|\boldsymbol{\nabla} W^+|^2\right]+ \frac{E_0^2(\varepsilon^+{-}\varepsilon^-)\eta_x^2}{1+\eta_x^2}, \end{align}

respectively, where the constant $E_0^2(\varepsilon ^+-\varepsilon ^-)$ in (2.11) has been neglected since it can always be absorbed by redefining $\phi ^\pm$ in the dynamic boundary condition. Finally, (2.1a,b), (2.2a,b), (2.5), (2.6), (2.7), (2.10) and the continuity condition $W^+=W^-$ at the interface complete the whole problem.

Figure 1. The schematic of the problem. The dashed line represents the undisturbed interface.

2.2. Zakharov–Craig–Sulem formulation

The Dirichlet–Neumann operator (DNO), which sends boundary value data to normal derivative data via solving the Laplace equation, is essential for investigating free boundary problems in potential theory. We define the velocity potentials and voltage potential at the interface as

(2.12)\begin{equation} \left. \begin{array}{c@{}} \xi^+(x,t) = \phi^+(x,\eta(x,t),t), \quad \xi^-(x,t)=\phi^-(x,\eta(x,t),t), \\ W(x,t) = W^+(x,\eta(x,t),t)=W^-(x,\eta(x,t),t). \end{array}\right\} \end{equation}

First, it follows from the kinematic boundary condition and the definitions of $\xi ^\pm$ that at the interface,

(2.13a,b)\begin{equation} \phi^\pm_x=\frac{\xi^\pm_x-\eta_x\eta_t}{1+\eta_x^2} \quad\text{and}\quad \phi^\pm_y= \frac{\eta_t+\eta_x\xi^\pm_x}{1+\eta_x^2}. \end{equation}

Following Craig & Sulem (Reference Craig and Sulem1993), we introduce the DNOs

(2.14)\begin{equation} \left. \begin{array}{c@{}} G^+(\eta,h^+)\,\xi^+ = [\eta_x\phi^+_x-\phi^+_y]_{y=\eta}=[\boldsymbol{\nabla} \phi^+\boldsymbol{\cdot} (-\boldsymbol{n})\sqrt{1+\eta_x^2}]_{y=\eta}, \\ G^-(\eta,h^-)\,\xi^- = [\phi^-_y-\eta_x\phi_x^-]_{y=\eta}=[\boldsymbol{\nabla} \phi^-\boldsymbol{\cdot} \boldsymbol{n}\sqrt{1+\eta_x^2}]_{y=\eta}, \\ G^+(\eta,h^+)\,W = [\eta_xW_x^+{-}W_y^+]_{y=\eta}=[\boldsymbol{\nabla} W^+\boldsymbol{\cdot} (-\boldsymbol{n})\sqrt{1+\eta_x^2}]_{y=\eta}, \\ G^-(\eta,h^-)\,W = [W^-_y-\eta_xW_x^-]_{y=\eta}=[\boldsymbol{\nabla} W^-\boldsymbol{\cdot} \boldsymbol{n}\sqrt{1+\eta_x^2}]_{y=\eta}. \end{array} \right\} \end{equation}

We suppress the dependency of DNOs on $h^\pm$ and $\eta$ in subsequent analyses for simplicity of notation. Then the derivatives of $W^\pm$ can be expressed as

(2.15a,b)\begin{equation} W^\pm_x=\frac{W_x\pm\eta_xG^\pm W}{1+\eta_x^2}\quad\text{and}\quad W^\pm_y=\frac{\mp G^\pm W+\eta_xW_x}{1+\eta_x^2}, \end{equation}

the kinematic boundary condition can be rewritten as a compact form

(2.16)\begin{equation} \eta_t=G^-\xi^- = -G^+\xi^+, \end{equation}

and the continuity equation for the electric displacement field at the interface now reads

(2.17)\begin{equation} \varepsilon^-G^-W+\varepsilon^+G^+W= (\varepsilon^--\varepsilon^+)\eta_x. \end{equation}

Following Zakharov (Reference Zakharov1968), who formulated the free-surface water wave problem in terms of surface quantities, we can rewrite the pressure equation as

(2.18)\begin{align} &\rho^-\left[\xi^-_t+\tfrac12(\xi^-_x)^2- \frac{(\eta_t+\eta_x\xi^-_x)^2}{2(1+\eta_x^2)}\right]- \rho^+\left[\xi^+_t+\tfrac12(\xi^+_x)^2- \frac{(\eta_t+\eta_x\xi^+_x)^2}{2(1+\eta_x^2)}\right] \nonumber\\ &\quad +\frac{E_0^2}{2(1+\eta_x^2)}[\varepsilon^-(G^-W)^2- \varepsilon^+(G^+W)^2-(\varepsilon^--\varepsilon^+) (W_x^2+2W_x)] \nonumber\\ &\quad +(\rho^--\rho^+)g\eta-\frac{\sigma\eta_{xx}}{(1+\eta_x^2)^{3/2}}=0, \end{align}

upon substituting (2.11), (2.13a,b), and (2.15a,b) into (2.7), and noticing that $\xi ^\pm _t=\phi ^\pm _t+\eta _t\phi ^\pm _y$. We denote by $\xi =\rho ^-\xi ^--\rho ^+\xi ^+$ the potential difference at the interface, and it is easy to see that (2.18) is the evolution equation for $\xi$. Recalling the kinematic boundary condition (2.16), one obtains the relations between $\xi$ and $\xi ^\pm$ as

(2.19)\begin{gather} G^+\xi =(\rho^-G^+{+}\rho^+G^-)\xi^-\quad\Longrightarrow\quad\xi^- = (\rho^-G^+{+}\rho^+G^-)^{{-}1}G^+\xi, \end{gather}

(2.20)\begin{gather}G^-\xi ={-}(\rho^-G^+{+}\rho^+G^-)\xi^+\quad\Longrightarrow\quad\xi^+ = -(\rho^-G^+{+}\rho^+G^-)^{{-}1}G^-\xi. \end{gather}

On the other hand, due to (2.17), the voltage potential $W$ and the interface displacement are related through

(2.21)\begin{equation} W=(\varepsilon^--\varepsilon^+) (\varepsilon^-G^-{+}\varepsilon^+G^+)^{{-}1}\eta_x. \end{equation}

Finally, by replacing $\xi ^\pm$ and $W$ with (2.19)–(2.21), equations (2.16) and (2.18) form a closed system with two unknowns: $\eta$ and $\xi$.

2.3. Linear theory

Coifman & Meyer (Reference Coifman and Meyer1985) proved that if the $L^1$-norm and Lipschitz-norm of $\eta$ are smaller than a certain constant, then $G^-$ is an analytic function of $\eta$. It then follows that $G^-(\eta,h^-)$ can be written naturally in the form of a Taylor expansion in $\eta$, namely $G^-(\eta,h^-)=\sum _{j=0}^\infty G^-_j(\eta,h^-)$, and Craig & Sulem (Reference Craig and Sulem1993) initially obtained a recursive formula for the expansion. The first three terms of the Taylor series are given by

(2.22)\begin{equation} \left. \begin{array}{c@{}} G_0^- = (-\partial_{xx})^{1/2}\tanh((-\partial_{xx})^{1/2}h^-), \\ G_1^- = -\partial_x\eta\partial_x-G^-_0\eta G^-_0, \\ G_2^- = \dfrac12G^-_0\eta^2\partial_{xx}+\dfrac12\partial_{xx}\eta^2G^-_0+G^-_0\eta G^-_0\eta G^-_0, \end{array} \right\} \end{equation}

where the dependency of $G_j^-$ on $\eta$ and $h^-$ has been suppressed for simplicity as usual. In the same vein, the DNO for the upper-layer fluid, $G^+(\eta,h^+)$, can be expanded as

(2.23)\begin{equation} G^+(\eta,h^+)=\sum_{j=0}^\infty G_j^+(\eta,h^+)=\sum_{j=0}^\infty({-}1)^j\,G^-_j(\eta,h^+), \end{equation}

which is obtained by replacing $\eta$ and $h^-$ with $-\eta$ and $h^+$, respectively, in $G_j^-(\eta,h^-)$.

In the subsequent analyses, we derive the linear dispersion relation of the problem. We first linearise the whole system around the trivial uniform stream solution $\eta =\phi\! =W\!=0$. Dropping nonlinear terms in (2.16), (2.18) and (2.21) gives

(2.24)\begin{equation} \left. \begin{array}{c@{}} \eta_t=G^-_0(\rho^-G_0^+{+}\rho^+G_0^-)^{{-}1}G^+_0\xi, \\ \xi_t-E_0^2(\varepsilon^--\varepsilon^+)W_x+ (\rho^--\rho^+)g\eta-\sigma\eta_{xx}=0, \\ W=(\varepsilon^--\varepsilon^+)(\varepsilon^-G^-_0+\varepsilon^+G^+_0)^{{-}1}\eta_x. \end{array}\right\}\end{equation}

For wavenumber $k$ and frequency $\omega$, substituting the ansatz

(2.25)\begin{equation} (\eta, \xi, W)=(\hat{\eta}, \hat{\xi}, \hat{W})\, {\rm e}^{\mathrm{i}(kx-\omega t)},\end{equation}

into the system (2.24) yields

(2.26)\begin{equation} \left. \begin{array}{c@{}} -\mathrm{i}\omega \hat{\eta}= \dfrac{|k|\tanh(|k|h^-)\tanh(|k|h^+)}{\rho^-\tanh(|k|h^+)+ \rho^+\tanh(|k|h^-)}\,\hat{\xi}, \\ {[}\varepsilon^+|k|\tanh(|k|h^+)+\varepsilon^-|k|\tanh(|k|h^-)] \hat{W}=(\varepsilon^--\varepsilon^+)\mathrm{i}k \hat{\eta}, \\ -\mathrm{i}\omega \hat{\xi}+(\rho^--\rho^+)g \hat{\eta}-E_0^2(\varepsilon^--\varepsilon^+)\mathrm{i}k \hat{W}+\sigma k^2 \hat{\eta}=0. \end{array}\right\} \end{equation}

This is a homogeneous linear algebraic system for $\hat {\eta }$, $\hat {\xi }$ and $\hat {W}$. By requiring the solution for this system to be non-trivial, we obtain, after some algebra, the linear dispersion relation in the form

(2.27)\begin{equation} \omega^2=\frac{|k|\tanh(|k|h^-)\tanh(|k|h^+)}{\mathcal{A}(|k|)} \left[(\rho^--\rho^+)g+\frac{E_0^2(\varepsilon^-- \varepsilon^+)^2|k|}{\mathcal{B}(|k|)}+\sigma|k|^2\right], \end{equation}

with

(2.28)\begin{gather} \mathcal{A}(|k|) = \rho^-\tanh(|k|h^+)+\rho^+\tanh(|k|h^-), \end{gather}

(2.29)\begin{gather}\mathcal{B}(|k|) = \varepsilon^+\tanh(|k|h^+)+\varepsilon^-\tanh(|k|h^-). \end{gather}

In the long-wave regime, i.e. the wavenumber $k$ is close to zero, the dispersion relation (2.27) can be approximated by

(2.30)\begin{equation} \omega^2\approx\frac{h^-h^+}{\rho^-h^+{+}\rho^+h^-}\left[\rho^--\rho^+{+}\frac{E_0^2(\varepsilon^--\varepsilon^+)^2}{\varepsilon^+h^+{+}\varepsilon^-h^-}\right]|k|^2, \end{equation}

which shows clearly the stabilising effect of the electric field. However, if one of the fluid layers is of infinite depth, then the electric field can provide dispersive regularisation only to short waves. For example, we let $h^-\rightarrow \infty$, and then the leading-order dispersion relation near $k=0$ reads $\omega ^2\approx |k|^2(\rho ^--\rho ^+)h^+/\rho ^+$, indicating that the system is linearly ill-posed when $\rho ^+>\rho ^-$ (namely the Rayleigh–Taylor unstable regime).

3.1. Equations for the shallow–shallow configuration

In what follows, we analyse the problem when the thickness of the upper fluid layer is comparable to that of the lower layer, namely $h=h^+/h^-=O(1)$. This means that both layers are shallow in comparison with the typical length of the interfacial wave. We derive weakly nonlinear models in various asymptotic regimes, each of which features a nonlinearity-dispersion balance.

We first consider the classic Boussinesq scaling: the interfacial wave varies on temporal and horizontal spatial scales of $1/\mu$, and the amplitude of the interface displacement scales by $\mu ^2$. All the parameters, including the Bond number, electric Bond number, depth ratio and permittivity ratio, are assumed to be of order $\mu ^0$ to enable the forces exerted by the electric field, surface tension and gravity to compete with each other. In summary, the following scales are chosen:

(3.3)\begin{equation} \left. \begin{array}{c@{}} \partial_x=O(\mu), \quad\partial_t=O(\mu), \quad\eta=O(\mu^2), \quad\xi=O(\mu), \\ R=O(1), \quad E_b=O(1),\quad\tau=O(1), \quad\varepsilon=O(1). \end{array}\right\} \end{equation}

We use order-of-magnitude arguments to derive reduced nonlinear models. First, due to the smallness of the spatial derivative, the DNOs can be expanded as

(3.4)\begin{equation} \left. \begin{array}{c@{}} G^- = -\partial_{xx}-\dfrac13\partial_{xxxx}-\partial_x\eta\partial_x+O(\mu^6), \\ G^+ = -h\partial_{xx}-\dfrac{h^3}{3}\partial_{xxxx}+\partial_x\eta\partial_x+O(\mu^6), \end{array}\right\}\end{equation}

and then a direct symbolic calculation yields

(3.5)\begin{align} (G^+{+}RG^-)^{{-}1}&=\left\{-(h+R)\partial_{xx}\left[1+ \frac{h^3+R}{3(h+R)}\,\partial_{xx}-\frac{1-R}{h+R}\, \partial_x^{{-}1}\eta\partial_x+O(\mu^4)\right]\right\}^{{-}1}\nonumber\\ &={-}\frac{1}{h+R}\left[1+\frac{h^3+R}{3(h+R)}\,\partial_{xx}- \frac{1-R}{h+R}\,\partial_x^{{-}1}\eta\partial_x+O(\mu^4)\right]^{{-}1}\partial_{xx}^{{-}1}\nonumber\\ &={-}\frac{\partial_{xx}^{{-}1}}{h+R}+\frac{h^3+R}{3(h+R)^2}- \frac{1-R}{(h+R)^2}\,\partial_x^{{-}1}\eta\partial_x^{{-}1}+O(\mu^2). \end{align}

In the same vein, one can obtain

(3.6)\begin{equation} (\varepsilon G^+{+}G^-)^{{-}1}={-}\frac{\partial_{xx}^{{-}1}}{\varepsilon h+1}+ \frac{\varepsilon h^3+1}{3(\varepsilon h+1)^2}+\frac{1- \varepsilon}{(\varepsilon h+1)^2}\,\partial_x^{{-}1}\eta\partial_x^{{-}1}+O(\mu^2). \end{equation}

With the aid of these formulae, we have the following expansions, after some calculations of (2.19)–(2.21):

(3.7a,b)\begin{equation} \xi^- = \frac{h\xi}{h+R}+O(\mu^2), \quad\xi^+ = \frac{\xi}{h+R}+O(\mu^2), \end{equation}

and

(3.8)\begin{equation} W_x={-}\frac{1-\varepsilon}{\varepsilon h+1}\,\eta+\frac{(1-\varepsilon)(\varepsilon h^3+1)}{3(\varepsilon+1)^2}\,\eta_{xx}+\frac{(1-\varepsilon)^2}{(\varepsilon h+1)^2}\,\eta^2+O(\mu^4). \end{equation}

Substituting the above formulae into boundary conditions at the interface, and retaining terms valid up to $O(\mu ^5)$ in (2.16) and $O(\mu ^4)$ in (3.2), one then obtains

(3.9)\begin{gather} \eta_t={-}\frac{h}{h+R}\,\xi_{xx}-\frac{h^2(1+hR)}{3(h+R)^2}\,\xi_{xxxx}- \frac{h^2-R}{(h+R)^2}\,(\eta\xi_x)_x,\end{gather}

(3.10)\begin{align} \xi_t&={-}\left[1-R+\frac{E_b(1-\varepsilon)^2}{\varepsilon h+1}\right] \eta-\frac{h^2-R}{2(h+R)^2}\,\xi_x^2+\frac{3E_b(1- \varepsilon)^3}{2(\varepsilon h+1)^2}\,\eta^2\nonumber\\ &\quad +\left[\frac{E_b(1-\varepsilon)^2(\varepsilon h^3+1)}{3(\varepsilon h+1)^2}+\tau\right] \eta_{xx}, \end{align}

a Boussinesq-type system in the presence of a horizontal electric field. From the system (3.9) and (3.10), we can infer that the electric field can provide dispersive regularisation to both gravity and surface tension, and introduce new nonlinearity.

The KdV equation can be derived from (3.9) and (3.10) upon further simplifications. Taking the derivative of (3.10) with respect to $x$ yields

(3.11)\begin{gather} \eta_t+u_x={-}\frac{h(1+hR)}{3(h+R)}\,u_{xxx}- \frac{h^2-R}{h(h+R)}\,(\eta u)_x,\end{gather}

(3.12)\begin{align} u_t+c^2\eta_x&={-}\frac{h^2-R}{2h(h+R)}\,(u^2)_x+ \frac{3hE_b(1-\varepsilon)^3}{2(\varepsilon h+1)^2(h+R)}\,(\eta^2)_x\nonumber\\ &\quad +\frac{h}{h+R}\left[\frac{E_b(1-\varepsilon)^2(\varepsilon h^3 +1)}{3(\varepsilon h+1)^2}+\tau\right] \eta_{xxx}, \end{align}

where

(3.13a,b)\begin{equation} u=\frac{h\xi_x}{h+R}\quad\text{and}\quad c^2=\frac{h}{h+R}\left[1-R+\frac{E_b(1-\varepsilon)^2}{\varepsilon h+1}\right]. \end{equation}

Note that the leading order of the system is a classic wave equation $\eta _{tt}-c^2\eta _{xx}=0$. Assuming that the wave travels mainly towards one direction (say, the right direction), the solution must take the form $\eta (x-ct)+\cdots$, and it then follows that $u=c\eta +\tilde {u}(x-ct)$, where $\tilde {u}$ is of order $\mu ^4$. Upon substituting $u$ with $c\eta +\tilde {u}$ in (3.11) and (3.12), and eliminating the leading-order terms via subtracting one equation from the other, we find

(3.14)\begin{align} \tilde{u}&={-}\frac{h}{2c(h+R)}\left[\frac{c^2(1+hR)}{3}+ \frac{E_b(1-\varepsilon)^2(\varepsilon h^3+1)}{3(\varepsilon h+1)^2}+\tau\right]\eta_{xx}\nonumber\\ &\quad -\frac{1}{2c}\left[\frac{c^2(h^2-R)}{2h(h+R)}+ \frac{3hE_b(1-\varepsilon)^3}{2(\varepsilon h+1)^2(h+R)}\right]\eta^2. \end{align}

Finally, substituting (3.14) into (3.11) or (3.12) yields the KdV equation

(3.15)\begin{equation} \eta_t+c\eta_x+\alpha\eta_{xxx}+\beta\eta\eta_x=0, \end{equation}

where

(3.16)\begin{gather} \alpha ={-}\frac{h}{2c(h+R)}\left[\frac{E_b(1-\varepsilon)^2(\varepsilon h^3+1)}{3(\varepsilon h+1)^2}+\tau-\frac{c^2(1+hR)}{3}\right], \end{gather}

(3.17)\begin{gather}\beta = \frac{3}{2c}\left[\frac{c^2(h^2-R)}{h(h+R)}-\frac{hE_b(1-\varepsilon)^3}{(\varepsilon h+1)^2(h+R)}\right], \end{gather}

are the dispersive and nonlinear coefficients, respectively. The famous soliton solution admitted by (3.15) is given by

(3.18)\begin{equation} \eta=\frac{12\alpha}{\beta}\,\delta^2\, {\rm sech}^2[\delta(x-ct-4\delta^2\alpha t-x_0)],\end{equation}

where $\delta$ and $x_0$ are arbitrary constants associated with the wave amplitude and initial phase, respectively.

It is well known that the soliton formation arises from a dynamic balance between dispersive and nonlinear effects. However, the dispersive and nonlinear coefficients in (3.15), $\alpha$ and $\beta$, may be close to zero for certain parameter sets, hence rescaling is required in the asymptotic analysis. For the first case, $\beta \ll 1$, a higher-order nonlinearity needs to be introduced to balance the dispersion, resulting in the modified KdV equation. For this purpose, we choose the scales

(3.19a–e)\begin{equation} \partial_x=O(\mu), \quad\partial_t=O(\mu), \quad\eta=O(\mu), \quad\xi=O(1), \quad\beta=O(\mu), \end{equation}

and the other parameters ($h$, $R$, $\varepsilon$ and $E_b$) keep the original order of magnitude. The new scaling indicates that the balance between the two effects can be achieved for larger-amplitude waves, though the temporal and spatial scales for solitons remain the same. Following the same procedure, symbolic calculations yield

(3.20)\begin{align} \left. \begin{array}{c@{}} \xi^-_x = \dfrac{h}{h+R}\,\xi_x-\dfrac{R(1+h)}{(R+h)^2}\,\eta\xi_x+O(\mu^3), \\ \xi^+_x ={-}\dfrac{1}{h+R}\,\xi_x-\dfrac{1+h}{(R+h)^2}\,\eta\xi_x+O(\mu^3), \\ W_x ={-}\dfrac{1-\varepsilon}{1+\varepsilon h}\,\eta+\dfrac{(1-\varepsilon)^2}{(1+\varepsilon h)^2}\,\eta^2+\dfrac{(1-\varepsilon)(1+\varepsilon h^3)}{3(1+\varepsilon h)^2}\,\eta_{xx}-\dfrac{(1-\varepsilon)^3}{(1+\varepsilon h)^3}\,\eta^3+O(\mu^4). \end{array}\right\} \end{align}

Substituting the above formulae into the kinematic and dynamic boundary conditions, namely (2.16) and (3.2), and then truncating the expansions at appropriate orders in $\mu$, one obtains

(3.21)\begin{gather} \eta_t={-}\frac{h}{h+R}\,\xi_{xx}+\frac{R-h^2}{(h+R)^2}\, (\eta\xi_x)_x-\frac{h^2(1+hR)}{3(h+R)^2}\,\xi_{xxxx}+\frac{R(1+h)^2}{(h+R)^3}\, (\eta^2\xi_x)_x, \end{gather}

(3.22)\begin{align} \xi_t&={-}\left[1-R+\frac{E_b(1-\varepsilon)^2}{1+\varepsilon h}\right] \eta+\frac{R-h^2}{2(h+R)^2}\,\xi_x^2+\frac{3E_b(1-\varepsilon)^3}{2(1+ \varepsilon h)^2}\,\eta^2\nonumber\\ &\quad +\left[\frac{E_b(1-\varepsilon)^2(1+\varepsilon h^3)}{3(1+\varepsilon h)^2}+\tau\right]\eta_{xx}+\frac{R(1+h)^2}{(h+R)^3}\,\eta\xi_x^2- \frac{2E_b(1-\varepsilon)^4}{(1+\varepsilon h)^3}\,\eta^3. \end{align}

By assuming one-way wave propagation along a major direction (right-going waves, say) and following a standard argument, one finally, after some tedious calculation, obtains a modified KdV equation

(3.23)\begin{equation} \eta_t+c\eta_x+\alpha\eta_{xxx}+\beta\eta\eta_x+\tilde{\beta}\eta^2\eta_x=0,\end{equation}

where $c$ is the long-wave speed defined in (3.13a,b), $\alpha$ and $\beta$ are the same as given in (3.16) and (3.17), and

(3.24)\begin{equation} \tilde{\beta}=\frac{3hE_b(1-\varepsilon)^4}{c(1+\varepsilon h)^3(h+R)}+ \frac{3E_b(1-\varepsilon)^3(R-h^2)}{c(1+\varepsilon h)^2(h+R)^2}- \frac{3cR(1+h)^2}{h(h+R)^2}. \end{equation}

For the second case, $\alpha \ll 1$, a fifth-order derivative term is commonly introduced to balance the quadratic nonlinearity, implying the scaling

(3.25a–e)\begin{equation} \partial_x = O(\mu), \quad\partial_t = O(\mu), \quad\eta=O(\mu^4), \quad\xi=O(\mu^3), \quad\alpha=O(\mu^2). \end{equation}

Following the same argument, after some algebra, a fifth-order KdV equation can be obtained to govern the interface dynamics:

(3.26)\begin{equation} \eta_t+c\eta_x+\alpha\eta_{xxx}+\beta\eta\eta_x+\tilde{\alpha}\eta_{xxxxx}=0, \end{equation}

where

(3.27)\begin{align} \tilde{\alpha}&= \frac{hc}{h+R}\left[\frac{1+h^3R}{15}- \frac{R(1-h^2)^2}{18(h+R)}\right]-\frac{h^2(1+hR)\tau}{6c(h+R)^2}\nonumber\\ &\quad +\frac{hE_b(1-\varepsilon)^2}{2c(h+R)}\left[\frac{(1+\varepsilon h^3)^2}{9(1+\varepsilon h)^3}-\frac{2(1+\varepsilon h^5)}{15(1+\varepsilon h)^2}-\frac{h(1+hR)(1+\varepsilon h^3)}{9(h+R)(1+\varepsilon h)^2}\right]. \end{align}

There is a third possibility: $\alpha$ and $\beta$ are very close to zero for the same set of parameters. This is possible even when the electric field is absent, as long as one chooses $R\approx h^2$ and $\tau \approx (1-h)(1+h^3)/3$. In this situation, one can obtain a higher-order modified KdV equation, including fifth-order dispersion and cubic nonlinearity (see Laget & Dias (Reference Laget and Dias1997) for the case when the electric field is absent). We omit the derivation of the model equation since it is tedious but straightforward if one follows the previous procedure.

3.2. Equation for the case of great upper-layer depth

In this part, we consider the case that the typical wavelength $\lambda$ is much larger than the mean depth of the lower layer but comparable to the thickness of the upper layer, which is termed the Benjamin scaling. The non-dimensionalisation is achieved in the same way as in the shallow–shallow case, where $h^-$ is taken as the length scale. Based on the linear analysis, $R$ is always set to be less than 1 in the current subsection, such that the problem is linearly well-posed. The Benjamin regime includes the additional expectation of small-amplitude motions and large interfacial tension (Benjamin Reference Benjamin1992). One chooses the scaling

(3.28)\begin{equation} \left.\begin{array}{c@{}} \partial_x=O(\mu), \quad\partial_t=O(\mu), \quad\eta=O(\mu), \quad\xi=O(1), \quad h\geqslant O(1/\mu), \\ R=O(1), \quad\varepsilon=O(1), \quad E_b=O(1), \quad\tau=O(1/\mu). \end{array}\right\} \end{equation}

In this case, the expansions of the DNOs can be expressed as

(3.29)\begin{gather} G^- = -\partial_{xx}-\tfrac{1}{3}\partial_{xxxx}-\partial_x\eta\partial_x+O(\mu^6), \end{gather}

(3.30)\begin{gather}G^+ = G^+_0+\partial_{x}\eta\partial_{x}+G_0^+\eta G_0^+{+}O(\mu^5), \end{gather}

where

(3.31)\begin{equation} G_0^+ = (-\partial_{xx})^{1/2}\tanh(h(-\partial_{xx})^{1/2})=O(\mu). \end{equation}

It then follows that

(3.32)\begin{gather} (G^+{+}RG^-)^{{-}1} = (G_0^+)^{{-}1}+R(G_0^+)^{{-}1}\partial_{xx}(G_0^+)^{{-}1}+O(\mu), \end{gather}

(3.33)\begin{gather}(\varepsilon G^+{+}G^-)^{{-}1} = \frac{1}{\varepsilon}(G_0^+)^{{-}1}+O(1). \end{gather}

Substituting these expansions into boundary conditions at the interface gives

(3.34)\begin{gather} \eta_t ={-}\xi_{xx}-(\eta\xi_x)_x-R\mathcal{P}\xi_{xxx}, \end{gather}

(3.35)\begin{gather}\xi_t ={-}\frac12\xi^2_x+\frac{E_b(1-\varepsilon)^2}{\varepsilon}\,\mathcal{P} \eta_x-(1-R)\eta+\tau\eta_{xx}, \end{gather}

where $\mathcal {P}=\partial _x(G_0^+)^{-1}$ is a pseudo-differential operator. In the limiting case $h\rightarrow \infty$, this operator becomes the Hilbert transform $\mathcal {H}$ defined by

(3.36)\begin{equation} \mathcal{H}[f](x) = \frac{1}{\rm \pi}\,\mathrm{P.V.} \int\frac{f(x')}{x'-x}\,{{\rm d} x}', \end{equation}

where ‘$\mathrm {P.V.}$’ means that the integral is in the Cauchy principal value sense. Finally, by assuming one-way wave propagation and following the same procedure as the previous subsection, one can reduce the system to a single evolution equation for $\eta$:

(3.37)\begin{equation} \eta_t+c\eta_x+\frac{3c}{2}\eta\eta_x+\left[\frac{cR}{2}- \frac{E_b(1-\varepsilon)^2}{2\varepsilon c}\right]\mathcal{P}\eta_{xx}- \frac{\tau}{2c}\,\eta_{xxx}=0, \end{equation}

where $c=\sqrt {1-R}$ is the long-wave speed. For the limiting case $h\rightarrow \infty$, (3.37) reduces to the famous Benjamin equation for long interfacial waves with strong surface tension (Benjamin Reference Benjamin1992). In the context of electrohydrodynamic waves, the Benjamin equation was first derived by Gleeson et al. (Reference Gleeson, Hammerton, Papageorgiou and Vanden-Broeck2007), who investigated the free-surface evolution of a conducting fluid under a normal electric field.

While weakly nonlinear dispersive models, such as (3.37), have been successful in many aspects of the free-surface/interfacial wave problems, strongly nonlinear models still need to be developed to describe moderate-amplitude wave phenomena as well as to afford a significant simplification over the primitive equations. Next, we derive strongly nonlinear models without the smallness assumption on the wave amplitude. The following scaling is made:

(3.38)\begin{equation} \left. \begin{array}{c@{}} \partial_x=O(\mu), \quad\partial_t=O(\mu), \quad\eta=O(1), \quad\xi=O(1/\mu), \\ R=O(1), \quad\varepsilon=O(1), \quad E_b=O(1), \quad\tau=O(1/\mu). \end{array} \right\} \end{equation}

Under this scaling assumption, one proceeds by expanding the DNOs in powers of $\mu$:

(3.39)\begin{gather} G^- = -\partial_{xx}-\partial_x\eta\partial_x+O(\mu^4), \end{gather}

(3.40)\begin{gather}G^+ = G^+_0+\partial_x\eta\partial_x+G^+_0\eta G^+_0+O(\mu^3). \end{gather}

We can then expand the crucial pseudo-differential operator $(G^++RG^-)^{-1}$ symbolically as

(3.41)\begin{align} &(G^+{+}RG^-)^{{-}1} = [G^+_0+\partial_x\eta\partial_x+G^+_0\eta G^+_0- R\partial_{xx}-R\partial_x\eta\partial_x+O(\mu^3)]^{{-}1}\nonumber\\ & = (G^+_0)^{{-}1}-\eta-(1-R)(G^+_0)^{{-}1}\partial_x\eta\partial_x (G^+_0)^{{-}1}+R(G^+_0)^{{-}2}\partial_{xx}+O(\mu). \end{align}

Similarly, we have

(3.42)\begin{equation} (\varepsilon G^+{+}G^-)^{{-}1}=\frac{1}{\varepsilon}(G^+_0)^{{-}1}+O(1)\quad\Longrightarrow\quad W_x=\frac{1-\varepsilon}{\varepsilon}(G^+_0)^{{-}1}\partial_{xx}\eta+O(\mu^2), \end{equation}

and

(3.43)\begin{gather} \xi^- = \xi+R(G^+_0)^{{-}1}\partial_{xx}\xi+R(G^+_0)^{{-}1} \partial_x\eta\partial_x\xi+O(\mu), \end{gather}

(3.44)\begin{gather}\xi^+ = (G^+_0)^{{-}1}\partial_{xx}\xi+(G^+_0)^{{-}1} \partial_x\eta\partial_x\xi+O(\mu). \end{gather}

Finally, we obtain the governing equations for $\eta$ and $\xi$ upon substituting these expansions into the kinematic and dynamic boundary conditions

(3.45)\begin{gather} \eta_t+u_{x}+(\eta u)_x-R[\mathcal{T}u+\mathcal{T}\eta u+\eta\mathcal{T}u+ \eta\mathcal{T}(\eta u)]_x=0, \end{gather}

(3.46)\begin{gather}u_t+(1-R)\eta_x+uu_x-\tau\eta_{xxx}+\frac{E_b(1- \varepsilon)^2}{\varepsilon}\,\mathcal{T}\eta_x-R [u\mathcal{T}u+u\mathcal{T}(\eta u)]_x=0, \end{gather}

where the pseudo-differential operator $\mathcal {T}$ is defined as

(3.47)\begin{equation} \mathcal{T}=\frac{(-\partial_{xx})^{1/2}}{\tanh (h(-\partial_{xx})^{1/2})}. \end{equation}

Note that the governing system (3.45) and (3.46) has been rewritten in a more revealing form via differentiating the dynamic boundary condition with respect to $x$ and letting $u=\xi _x$.

A slightly different nonlinear model can be derived by using $\xi ^-$ instead of $\xi$ as the unknown. Recalling the kinematic boundary condition (2.16) and the related operator expansions, one can express $\xi ^+$ as

(3.48)\begin{equation} \xi^+ = -(G^+)^{{-}1}\eta_t={-}(G^+_0)^{{-}1}\eta_t+O(\mu). \end{equation}

Substituting (3.48) into the dynamic boundary condition (3.2) and retaining terms valid up to $O(\mu )$, one obtains

(3.49)\begin{equation} \xi^-_t+R(G^+_0)^{{-}1}\eta_{tt}+\frac12(\xi^-_x)^2+(1-R)\eta+ \frac{E_b(1-\varepsilon)^2}{\varepsilon}\,\mathcal{T}\eta-\tau\eta_{xx}=0. \end{equation}

The approximate kinematic boundary condition is constructed by the retention of terms up to $O(\mu ^2)$ in the expansion of $G^-$. One finally obtains

(3.50)\begin{gather} \eta_t+u_x+(\eta u)_x=0, \end{gather}

(3.51)\begin{gather}u_t+(1-R)\eta_x+uu_x+\frac{E_b(1-\varepsilon)^2}{\varepsilon}\, \mathcal{T}\eta_x-\tau\eta_{xxx}+R\mathcal{T}[(1+\eta)u]_t=0, \end{gather}

where we have written $u = \xi ^-_x$ as before. The system (3.50) and (3.51) has an inbuilt advantage over the system (3.45) and (3.46) in computing solitary waves. To illustrate this, we assume that a solitary wave propagates with a constant speed $c$, and the wave becomes static in the moving frame $X=x-ct$. By integrating (3.50), one can show easily that $u={c\eta }/({1+\eta })$. After substituting this relation into (3.51), an integrodifferential equation for $\eta$ is obtained:

(3.52)\begin{equation} -\frac{c^2\eta(2+\eta)}{2(1+\eta)^2}+(1-R)\eta+ \left[\frac{E_b(1-\varepsilon)^2}{\varepsilon}-c^2R\right] \mathcal{T}\eta-\tau\eta_{XX}=0,\end{equation}

which can be solved by some numerical iteration method. However, (3.45) cannot relate $\eta$ and $u$ explicitly, complicating the problem.

3.3. Barannyk–Papageorgiou–Petropoulos scaling

Barannyk et al. (Reference Barannyk, Papageorgiou and Petropoulos2012, Reference Barannyk, Papageorgiou, Petropoulos and Vanden-Broeck2015) considered the problem of interfacial electrocapillary-gravity waves between two perfect dielectrics under a horizontal electric field. To understand the effects of the electric field on the Rayleigh–Taylor instability, they proposed a nonlinear long-wave model based on novel scaling. Their scaling assumes that interfacial waves are long relative to the upper-layer thickness but of comparable length with the lower-layer thickness. More precisely, the scales are chosen as

(3.53)\begin{equation} \left. \begin{array}{c@{}} \partial_x=O(1), \quad\partial_t=O(\mu), \quad \eta=O(\mu^2), \quad\xi^+ = O(\mu), \quad h=O(\mu^2), \\ R=O(1), \quad\varepsilon=O(1), \quad E_b=O(1), \quad\tau=O(1). \end{array} \right\} \end{equation}

Under the Barannyk–Papageorgiou–Petropoulos scaling, the kinematic boundary condition (2.16) indicates that

(3.54)\begin{equation} \xi^- = -(G^-)^{{-}1}G^+\xi^+ = -(G^-_0)^{{-}1}(h|D|^2+\partial_x\eta\partial_x)\xi^+{+}O(\mu^5), \end{equation}

where $G^-_0=|D|\tanh (|D|)$. It follows directly that $\xi ^-=O(\mu ^3)$. One proceeds by expanding the pseudo-differential operators in powers of $\mu$ in the boundary conditions at the interface, as before. Retaining terms valid up to $O(\mu ^3)$ for the kinematic boundary condition, and up to $O(\mu ^2)$ for the dynamic boundary condition, yields

(3.55)\begin{gather} \eta_t-h\xi^+_{xx}+(\eta\xi^+_x)_x=0, \end{gather}

(3.56)\begin{gather}\xi^+_t+\frac12(\xi^+_x)^2+\frac{E_b(1-\varepsilon)^2}{R}\,(G^-_0)^{{-}1} \eta_{xx}-\frac{1-R}{R}\,\eta+\frac{\tau}{R}\,\eta_{xx}=0. \end{gather}

Changing variables as

(3.57a–c)\begin{equation} \eta=h(1-\tilde{\eta}), \quad\xi^+_x=\sqrt{h}u, \quad \partial_t=\sqrt{h}\partial_\tau, \end{equation}

one arrives at the system

(3.58)\begin{gather} \tilde{\eta}_{\tau}+(u\tilde{\eta})_x = 0, \end{gather}

(3.59)\begin{gather}R(u_{\tau}+uu_x)-E_b(1-\varepsilon)^2(G_0^-)^{{-}1}\tilde{\eta}_{xxx} +(1-R)\tilde{\eta}_x-\tau\tilde{\eta}_{xxx}=0, \end{gather}

which is the same as obtained by Barannyk et al. (Reference Barannyk, Papageorgiou and Petropoulos2012). It is noted that, unlike the Benjamin equation, the dispersive effects in the system (3.58) and (3.59) come from only the electric field and surface tension. That is because the Barannyk–Papageorgiou–Petropoulos scaling gives rise to a very small $\xi ^-$, which excludes the motion of the deep layer when truncating the asymptotic expansion at an appropriate order. Though there is a smallness assumption on wave amplitude, the system can be viewed as a strongly nonlinear model since the amplitude is of the same order as the upper-layer thickness. The system (3.58) and (3.59) demonstrates interesting dynamics, including suppressing the Rayleigh–Taylor instability, self-similar solutions, and touch-up singularities, and interested readers are referred to Barannyk et al. (Reference Barannyk, Papageorgiou, Petropoulos and Vanden-Broeck2015) for more details.

4.1. Boundary integral equations

Define the complex velocity potentials as

(4.1)\begin{equation} f^{{\pm}}(z)=\phi^{{\pm}}(z)+\mathrm{i}\,\psi^{{\pm}}(z), \end{equation}

where $\psi ^{\pm }$ are stream functions as well as the complex conjugates of the corresponding velocity potentials, and $z=x+\mathrm {i}y$ stands for a point on the complex plane. Without loss of generality, we assume $\phi ^{\pm }=0$ at $x=0$, and $\psi ^{\pm }=0$ at the interface. It follows from the non-dimensionalisation scheme that $\psi ^-=-1$ at the bottom, and $\psi ^+=h$ at the top. Using a hodograph transformation that exchanges the independent and dependent variables, the lower (upper) layer is mapped conformally onto the infinite strip $-1<\psi ^-<0$ ($0<\psi ^+< h$) in the $\phi ^\pm +\mathrm {i}\psi ^\pm$ plane. The impermeability boundary conditions at two rigid walls can be satisfied automatically by using the method of images. Specifically, the function ${\mathrm {d}z}/{\mathrm {d}f^+}-1$ can be reflected in the line $\psi ^+=h$ to produce a function analytic in $0<\psi ^+<2h$, and the same procedure can be applied to the lower-layer fluid and the analytic function ${\mathrm {d}z}/{\mathrm {d}f^-}-1$. Since the velocity potential is discontinuous across the interface, an auxiliary function $\chi$ can be introduced such that $\chi (\phi ^-) = \phi ^+$ at the interface. Hereafter, we use $\phi ^-$ as the primary independent variable and suppress its superscript for ease of notation. Applying the Cauchy integral formula to both lower- and upper-layer fluids results in the two integral equations

(4.2)\begin{gather} x_\phi(\phi_m)-1=\frac{1}{{\rm \pi} }\int_{-\infty}^{\infty} \left[\frac{2(x_\phi-1) - y_\phi(\phi-\phi_m)}{(\phi-\phi_m)^2+4} -\frac{y_\phi}{\phi-\phi_m}\right] \mathrm{d}\phi, \end{gather}

(4.3)\begin{gather}\frac{x_\phi}{\chi_\phi}\,(\phi_m)-1=\frac{1}{{\rm \pi} }\int^{\infty}_{-\infty}\left[ \frac{2h(x_\phi-\chi_\phi) +y_\phi(\chi(\phi)-\chi(\phi_m))}{(\chi(\phi)- \chi(\phi_m))^2+4h^2}+\frac{y_\phi}{\chi(\phi)-\chi(\phi_m)}\right] \mathrm{d}\phi, \end{gather}

where the functions $x_\phi$ and $y_\phi$ are evaluated at the interface, and both $\phi$ and $\phi _m$ represent the horizontal coordinate of the $\phi ^-+\mathrm {i}\psi ^-$ plane. For the electric field, we introduce the complex voltage potentials

(4.4)\begin{equation} q^{{\pm}}(z) = W^{{\pm}}(z)+\mathrm{i}\,U^{{\pm}}(z), \end{equation}

where $W^{\pm }$ are defined in § 2.1, and $U^\pm$ are newly introduced complex conjugate functions. Based on the Cauchy integral formula, two integral equations are obtained:

(4.5)\begin{gather} W_\phi(\phi_m) = \frac{1}{{\rm \pi} }\int_{-\infty}^{\infty}\left[\frac{2W_\phi + U^-_\phi(\phi-\phi_m)}{(\phi-\phi_m)^2+4}+\frac{U^-_\phi}{\phi- \phi_m}\right]\mathrm{d}\phi, \end{gather}

(4.6)\begin{gather}\frac{W_\phi}{\chi_\phi}\,(\phi_m)=\frac{1}{{\rm \pi} } \int^{\infty}_{-\infty} \left[\frac{2hW_\phi - U^+_\phi(\chi(\phi)-\chi(\phi_m))}{(\chi(\phi)- \chi(\phi_m))^2+4h^2}-\frac{U^+_\phi}{\chi(\phi)-\chi(\phi_m)} \right]\mathrm{d}\phi, \end{gather}

where $W$ stands for the perturbation of the voltage potential at the interface, as before. In the transformed plane, the continuity condition of the electric displacement field and the dynamic boundary condition are recast to

(4.7)\begin{gather} (\varepsilon-1)y_\phi=\varepsilon U^+_\phi - U^-_\phi, \end{gather}

(4.8)\begin{align} &\frac{F^2}{2J}\,(1-R\chi_\phi^2)+(1-R) \left(y-\frac{F^2}{2}\right)+\frac{E_b}{2J}\, (U_\phi^{{-}2}-W_\phi^2)-\frac{\varepsilon E_b}{2J}\, (U_\phi^{{+}2}-W_\phi^2)\nonumber\\ &\quad -\frac{E_b(1-\varepsilon)}{J}\,x_\phi W_\phi- \frac{\tau}{J^{3/2}}\,(y_{\phi \phi} x_\phi - x_{\phi \phi} y_\phi)=0, \end{align}

respectively, where $J=x_\phi ^2+y_\phi ^2$ is the Jacobian of the transformation, and $F={c}/{\sqrt {gh^-}}$ is the Froude number measuring the wave speed.

The problem being considered is as follows: given one of the two parameters $H=y(0)$ and $F$, we seek solitary-wave solutions to the system (4.2), (4.3), (4.5)–(4.8). Since we confine our attention to symmetric solutions in the present paper, it is sufficient to perform numerical computations on a half real line (say, $0\le \phi <\infty$). Due to the decaying nature of solitary waves, a truncated domain $[0,L)$ with large $L$ is used in numerical computations. We introduce a set of mesh grids

(4.9)\begin{equation} \phi_i=\frac{(i-1)L}{N}, \quad i = 1,2,\ldots,N, \end{equation}

and the corresponding unknowns

(4.10)\begin{equation} \frac{\mathrm{d} y}{\mathrm{d}\phi}(\phi_i), \quad\frac{\mathrm{d} x}{\mathrm{d}\phi}(\phi_i), \quad\frac{\mathrm{d}\chi}{\mathrm{d}\phi}(\phi_i), \quad \frac{\mathrm{d}W}{\mathrm{d}\phi}(\phi_i), \quad\frac{\mathrm{d}U^-}{\mathrm{d}\phi}(\phi_i), \quad F; \end{equation}

$U^+_\phi$ can be obtained through (4.7). The symmetry of solutions gives

(4.11a,b)\begin{equation} y_\phi(0)=0, \quad U^-_\phi(0)=0, \end{equation}

and additionally, the decaying nature of solitary waves indicates that

(4.12a–c)\begin{equation} x_\phi(\phi_N)=1, \quad\chi_\phi(\phi_N)=1, \quad W_\phi(\phi_N)=0. \end{equation}

Overall, there are $5N-4$ unknowns, and we need the same number of equations. To avoid the singularities in computations of the Cauchy integrals, we introduce another set of mesh grids

(4.13)\begin{equation} \phi_{i+1/2}=\frac{\phi_i+\phi_{i+1}}{2}, \quad i=1,2,\ldots,N-1. \end{equation}

Evaluating (4.2)–(4.6) and (4.8) at the midpoints $\phi _{i+1/2}$ results in $5N-5$ algebraic equations. We can specify the wave displacement at the centre to close the system, namely $\eta (0)=H$. It is noted that this condition should be omitted while computing a solitary wave with a prescribed Froude number. Here, $x$, $y$ and $\chi$ can be obtained by integrating their derivatives using the trapezoidal rule. The derivatives $x_{\phi \phi }$ and $y_{\phi \phi }$ are computed via the five-point centred difference scheme. And the unknowns at the midpoints are calculated by a fourth-order interpolation formula.

The fully nonlinear equations are solved via a Newton iteration method. The program is considered to have converged to a solution when the $l^{\infty }-$norm of the residual error is less than $10^{-10}$. The initial guess of the iteration is obtained by applying an artificial pressure to the interface and then eliminating it gradually via a numerical continuation. Once a solitary-wave solution is found, other solutions on the same branch are computed via a straightforward continuation method by changing the Froude number, the wave amplitude or other parameters. Finally, the number of grid points is chosen sufficiently large ($N\geqslant 2000$), and the grid spacing is chosen sufficiently small ($L/N\leqslant 0.1$), to ensure accurate enough solutions.