2.1. Governing equations

In the interior regions of the fluids, the hydrodynamics of the system is governed by the Euler equations

(2.1)\begin{equation} \hat{\rho}_i \left(\frac{\partial \boldsymbol{\hat{u}_i}}{\partial \hat{t}} + \boldsymbol{\hat{u}_i} \boldsymbol{\cdot} \hat{\boldsymbol{\nabla}} \boldsymbol{\hat{u}_i} \right) ={-} \hat{\boldsymbol{\nabla}} \hat{p}_i - \hat{\rho}_i \hat{g} {\boldsymbol{j}}, \quad i=1,2 , \end{equation}

where $\hat {\boldsymbol {\nabla }} \equiv ({\partial }/{\partial \hat {x}},{\partial }/{\partial \hat {z}})$ is the dimensional spatial differential operator, and $\boldsymbol {\hat {u}_i}$ and $\hat {p}_i$ are the velocity and pressure of the fluid in phase $i$ ($i=1,2$), respectively. Since both fluids are incompressible and irrotational, we have

(2.2)\begin{gather} \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} \boldsymbol{\hat{u}_i} = 0, \quad i=1,2, \end{gather}

(2.3)\begin{gather}\hat{\boldsymbol{\nabla}} \times \boldsymbol{\hat{u}_i} = 0, \quad i=1,2. \end{gather}

For a perfectly dielectric fluid, the electrical conductivity is zero. This means that no free charges exist in the bulk of fluids and only polarized charges need to be counted. In this case, the dynamic current is very small and the induced magnetic fields are negligible. Consequently, Faraday's law reduces to

(2.4)\begin{equation} \hat{\boldsymbol{\nabla}} \times \boldsymbol{\hat{E}_i} = 0, \quad i=1,2, \end{equation}

where $\boldsymbol {\hat {E}_i}$ is the electric field in fluid $i$ ($i=1,2$). Due to the absence of free charges, Gauss's law can be written as

(2.5)\begin{equation} \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} (\hat{\epsilon}_i \boldsymbol{\hat{E}_i}) = 0, \quad i=1,2. \end{equation}

At the vertical far fields, the velocities of the fluids are zero, i.e.

(2.6)\begin{equation} \boldsymbol{\hat{u}_1}(\hat{z} = \infty) = \boldsymbol{\hat{u}_2}(\hat{z} ={-}\infty) = 0, \end{equation}

and the electric fields are constant in each fluid phase, i.e.

(2.7)\begin{equation} \boldsymbol{\hat{E}_1}(\hat{z}=\infty) = \boldsymbol{\hat{E}_2}(\hat{z}={-}\infty) = \hat{E}_h \boldsymbol{i}. \end{equation}

At the material interface, the fluid velocity is continuous in the normal direction, i.e.

(2.8)\begin{equation} \boldsymbol{\hat{u}_1} \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{\hat{u}_2} \boldsymbol{\cdot} \boldsymbol{n} = \frac{\partial \boldsymbol{\hat{X}}}{\partial \hat{t}} \boldsymbol{\cdot} \boldsymbol{n} \quad \mathrm{at}\ \hat{z}={\hat{\eta}}, \end{equation}

where $\boldsymbol {\hat {X}} \equiv (\hat {x}, \hat {\eta }(\hat {t},\hat {x}))$ denotes the location of the material interface at time $\hat {t}$ and $\boldsymbol {n} = (-\partial \hat {\eta } / \partial \hat {x}, 1) / \sqrt {1+(\partial \hat {\eta } / \partial \hat {x})^2}$ denotes the normal unit vector at the material interface (see figure 1). Other than the bound charges due to polarization, there are no unbounded surface charges at the material interface. Therefore, the normal component of the electric displacement field $\hat {\epsilon } \boldsymbol {\hat {E}}$ is continuous across the material interface, and so is the tangential component of the electric field $\boldsymbol {\hat {E}}$. These conditions lead to

(2.9)\begin{gather} \left[ \hat{\epsilon} \boldsymbol{\hat{E}} \boldsymbol{\cdot} \boldsymbol{n} \right]_2^1 = 0 \quad \mathrm{at} \ \hat{z}=\hat{\eta}, \end{gather}

(2.10)\begin{gather}\left[ \boldsymbol{\hat{E}} \boldsymbol{\cdot} \boldsymbol{s} \right]_2^1 = 0 \quad \mathrm{at} \ \hat{z}=\hat{\eta}, \end{gather}

where $\boldsymbol {s} = (1, \partial \hat {\eta } / \partial \hat {x}) / \sqrt {1+(\partial \hat {\eta } / \partial \hat {x})^2}$ is the tangential unit vectors at the material interface (see figure 1). The normal stress at the material interface is also continuous, i.e.

(2.11)\begin{equation} \boldsymbol{n} \boldsymbol{\cdot} \hat{\boldsymbol{\mathsf{T}}}_{\boldsymbol{\mathsf{1}}} \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{n} \boldsymbol{\cdot} \hat{\boldsymbol{\mathsf{T}}}_{\boldsymbol{\mathsf{2}}} \boldsymbol{\cdot} \boldsymbol{n} \quad \mathrm{at} \ \hat{z}=\hat{\eta}, \end{equation}

where

(2.12)\begin{equation} \hat{\boldsymbol{\mathsf{T}}}_{\boldsymbol{\mathsf{i}}} ={-}\hat{p}_i \boldsymbol{\mathsf{I}} + \hat{\epsilon} \boldsymbol{\hat{E}_i} \otimes \boldsymbol{\hat{E}_i} - \tfrac{1}{2} \hat{\epsilon}_i |\boldsymbol{\hat{E}_i}|^2 \boldsymbol{\mathsf{I}} \end{equation}

is the stress tensor for fluid $i$ and $\boldsymbol{\mathsf{I}}$ is the identity matrix. We comment that, on the right-hand side of (2.12), the first term is the hydrodynamic component, and the second and third terms are the electrodynamic components. The expressions of the electrodynamic components come from the Maxwell stress tensor. Here we have used the condition of constant permittivity.

For our system, (2.1)–(2.5) are the governing equations in the interior regions of the fluids, (2.6) and (2.7) are the boundary conditions at the far fields and (2.8)–(2.11) are the boundary conditions at the material interface. We comment that (2.1)–(2.3), (2.6) and (2.8) do not involve electric fields and are identical to the equations for the classical RTI. Since our system contains electric fields, it requires additional equations, namely (2.4), (2.5), (2.7) and (2.9)–(2.11).

2.2. Non-dimensional parameters and equations

The governing equations for our system given by (2.1)–(2.11) are expressed in terms of dimensional variables. Here we non-dimensionalize these governing equations and all physical parameters. In the paper, all dimensionless quantities are denoted by notations without the circumflex symbol.

We introduce the Atwood number for the fluid density

(2.13)\begin{equation} A = \frac{\hat{\rho}_1-\hat{\rho}_2}{\hat{\rho}_1+\hat{\rho}_2}. \end{equation}

Since in our setting the top fluid (fluid $1$) is always heavier than the bottom fluid (fluid $2$) and the gravity points downwards, $A > 0$ and gravity always acts as a destabilizing factor. Similarly, we introduce the ‘Atwood number’ for the electric permittivity, namely

(2.14)\begin{equation} A_\epsilon = \frac{\hat{\epsilon}_1-\hat{\epsilon}_2}{\hat{\epsilon}_1+\hat{\epsilon}_2}. \end{equation}

Note that when $\hat {\epsilon }_1>\hat {\epsilon }_2$, i.e. $A_\epsilon > 0$, the top fluid (fluid $1$) has higher permittivity and, therefore, yields a higher degree of polarization in response to the external electric field than the bottom fluid (fluid $2$). When $\hat {\epsilon }_1<\hat {\epsilon }_2$, i.e. $A_\epsilon < 0$, the top fluid has lower permittivity and yields a lower degree of polarization in response to the external electric field than the bottom fluid.

The length is non-dimensionalised by the wavenumber $\hat {k}$, namely

(2.15a–c)\begin{equation} x = \hat{x} \hat{k}, \quad z = \hat{z} \hat{k}, \quad \eta = \hat{\eta}\hat{k}. \end{equation}

It is easy to check that $\boldsymbol {n} = (-\partial {\eta } / \partial {x}, 1) / \sqrt {1+(\partial {\eta } / \partial {x})^2}$ and $\boldsymbol {s} = (1, \partial {\eta } / \partial {x}) / \sqrt {1+(\partial {\eta } / \partial {x})^2}$. There are two driving forces in the system: gravity $\hat {g}$ with $\hat {g}>0$, and the electrical force caused by the external electric field characterized by $\hat {f}_E = \hat {E}_h^2 \hat {k} ({\hat {\epsilon }_1+\hat {\epsilon }_2})/({\hat {\rho }_1+\hat {\rho }_2})$. The former generates the interfacial instability and the latter suppresses it (Eldabe Reference Eldabe1989). It is the competition between these two driving forces that governs the evolution of the material interface. To quantify the relative strength of these two driving forces, we introduce the dimensionless gravity

(2.16)\begin{equation} g = \frac{\hat{g}}{ \hat{g} + \hat{f}_E }, \end{equation}

and the dimensionless horizontal electric field

(2.17)\begin{equation} E_h = \mathrm{sgn}(\hat{E}_h) \times \sqrt{\frac{\hat{f}_E}{\hat{g} + \hat{f}_E}}, \end{equation}

where $\mathrm {sgn}({\cdot })$ is the sign function. We comment that $g$ denotes the fraction of the driving forces due to gravity, and $E_h^2$ denotes the fraction of the driving forces due to the electrical force. Hence, when the system is mostly driven by gravity and the electrical force is negligible, we have $g \to 1$ and $E_h^2 \to 0$. On the other hand, when the electrical force is dominating, we have $g \to 0$ and $E_h^2 \to 1$.

The presence of the two different driving forces gives two distinctive characteristic time scales: one is the time scale $(\hat {g} \hat {k})^{-1/2}$ due to gravity alone, and the other is the time scale $(\hat {f}_E \hat {k})^{-1/2}$ due to the electrical force alone. This complicates the non-dimensionalisation of time. We introduce a combined characteristic time $\hat {t}_c = (\hat {g} \hat {k} + \hat {f}_E \hat {k})^{-1/2}$, which includes both the gravitational and electrical force time scales. This combined form avoids the characteristic time becoming zero when one of the driving forces is absent, i.e. when $\hat {g} = 0$ or $\hat {E}_h = 0$. Therefore, $\hat {t}_c$ can serve as the characteristic time of our system for all cases. Then the time, the pressure, the velocity field and the electric field are non-dimensionalised as

(2.18a–d)\begin{equation} t = \hat{t} / \hat{t}_c, \quad p_i = \hat{p}_i \hat{k}^2 \hat{t}_c^2 / (\hat{\rho}_1+\hat{\rho}_2), \quad \boldsymbol{u_i} = \boldsymbol{\hat{u}_i} \hat{k} \hat{t}_c, \quad \boldsymbol{{E}_i} = \boldsymbol{\hat{E}_i}/\hat{E}_h. \end{equation}

The dimensionless stress tensors are expressed as

(2.19)\begin{gather} \boldsymbol{\mathsf{T}}_{\boldsymbol{\mathsf{1}}} ={-}{p_1} \boldsymbol{\mathsf{I}} + \frac{1+A_\epsilon}{2} E_h^2 \boldsymbol{{E}_1} \otimes \boldsymbol{{E}_1} - \frac{1+A_\epsilon}{4} E_h^2 |\boldsymbol{{E}_1}|^2 \boldsymbol{\mathsf{I}}, \end{gather}

(2.20)\begin{gather}\boldsymbol{\mathsf{T}}_{\boldsymbol{\mathsf{2}}} ={-}{p_2} \boldsymbol{\mathsf{I}} + \frac{1-A_\epsilon}{2} E_h^2 \boldsymbol{{E}_2} \otimes \boldsymbol{{E}_2} - \frac{1-A_\epsilon}{4} E_h^2 |\boldsymbol{{E}_2}|^2 \boldsymbol{\mathsf{I}}. \end{gather}

The dimensional governing equations and boundary conditions given by (2.1)–(2.11) can be expressed in terms of the non-dimensional quantities defined by (2.15a–c)–(2.18a–d) and the results are

(2.21)\begin{gather} \frac{1+A}{2} \left(\frac{\partial \boldsymbol{u_1}}{\partial t} + \boldsymbol{u_1} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u_1} \right) ={-} \boldsymbol{\nabla} p_1 - \frac{1+A}{2} g {\boldsymbol{j}}, \end{gather}

(2.22)\begin{gather}\frac{1-A}{2} \left(\frac{\partial \boldsymbol{u_2}}{\partial t} + \boldsymbol{u_2} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u_2} \right) ={-} \boldsymbol{\nabla} p_2 - \frac{1-A}{2} g {\boldsymbol{j}}, \end{gather}

(2.23)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{{u}_i} = 0, \quad i=1,2, \end{gather}

(2.24)\begin{gather}\boldsymbol{\nabla} \times \boldsymbol{{u}_i} = 0, \quad i=1,2, \end{gather}

(2.25)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{{E}_i} = 0, \quad i=1,2, \end{gather}

(2.26)\begin{gather}\boldsymbol{\nabla} \times \boldsymbol{{E}_i} = 0, \quad i=1,2, \end{gather}

(2.27)\begin{gather}\boldsymbol{u_1}(z = \infty) = \boldsymbol{u_2}(z ={-}\infty) = 0, \end{gather}

(2.28)\begin{gather}\boldsymbol{{E}_1}(z=\infty) = \boldsymbol{{E}_2}(z={-}\infty) = \boldsymbol{i}, \end{gather}

(2.29)\begin{gather}\boldsymbol{{u}_1} \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{{u}_2} \boldsymbol{\cdot} \boldsymbol{n} = \frac{\partial \boldsymbol{X}}{\partial t} \boldsymbol{\cdot} \boldsymbol{n} \quad \mathrm{at} \ {z}={\eta}, \end{gather}

(2.30)\begin{gather}(1+A_\epsilon) \boldsymbol{{E}_1} \boldsymbol{\cdot} \boldsymbol{n} = (1-A_\epsilon) \boldsymbol{{E}_2} \boldsymbol{\cdot} \boldsymbol{n} \quad \mathrm{at}\ {z}={\eta}, \end{gather}

(2.31)\begin{gather}\boldsymbol{{E}_1} \boldsymbol{\cdot} \boldsymbol{s} = \boldsymbol{{E}_2} \boldsymbol{\cdot} \boldsymbol{s} \quad \mathrm{at} \ {z}={\eta}, \end{gather}

(2.32)\begin{gather}\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\boldsymbol{\mathsf{1}}} \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\boldsymbol{\mathsf{2}}} \boldsymbol{\cdot} \boldsymbol{n} \quad \mathrm{at} \ {z}={\eta}, \end{gather}

where $\boldsymbol {\nabla } \equiv ({\partial }/{\partial x},{\partial }/{\partial z}) = {1}/{\hat {k}} ({\partial }/{\partial \hat {x}},{\partial }/{\partial \hat {z}})$ is the non-dimensional spatial differential operator and $\boldsymbol {{X}} \equiv ({x}, {\eta }({t},{x}))$ is the location of the material interface in terms of dimensionless length.

In the following sections all analyses and results are presented in terms of non-dimensional quantities unless explicitly stated.

3.1. Governing equations in terms of potentials

Since the fluids in our system are incompressible and irrotational (see (2.23) and (2.24)), there exists a velocity potential in each fluid phase. Furthermore, the electric fields in our system are both divergence free and rotation free in the interior of the fluids (see (2.25) and (2.26)), and there also exists a voltage potential for each fluid. In this subsection we rewrite the governing equations (2.21)–(2.32) in terms of these potentials. In terms of the dimensionless velocity potential $\phi _i$, the dimensionless velocity field in the interior of each fluid phase $i$ can be expressed as $\boldsymbol {{u}_i} = -\boldsymbol {\nabla } {\phi }_i$. Therefore, (2.23) can be rewritten as the Laplace equation for ${\phi }_i$,

(3.1)\begin{equation} {\nabla}^2 {\phi}_i = 0, \quad i=1,2. \end{equation}

Similarly, in terms of the dimensionless voltage potentials ${V}_i$, the dimensionless electric fields in the interior of the fluids can be expressed as $\boldsymbol {{E}_i} = -\boldsymbol {\nabla } {V}_i$. Then (2.25) leads to a Laplace equation for ${V}_i$,

(3.2)\begin{equation} {\nabla}^2 {V}_i = 0, \quad i=1,2. \end{equation}

The boundary conditions at the far fields, i.e. (2.27) and (2.28), are also rewritten in terms of dimensionless velocity potentials ${\phi }_i$ and dimensionless voltage potentials ${V}_i$ as

(3.3)\begin{gather} \boldsymbol{\nabla} {\phi}_1 ({z} = \infty) = \boldsymbol{\nabla} {\phi}_2 ({z} ={-}\infty)= 0, \end{gather}

(3.4)\begin{gather}\boldsymbol{\nabla} {V}_1 ({z}=\infty) = \boldsymbol{\nabla} {V}_2 ({z}={-}\infty) ={-} \boldsymbol{i}. \end{gather}

Since the dimensionless voltage potentials have been normalized by the strength of the external horizontal electric field $\hat {E}_h$, in (3.4) the boundary conditions no longer explicitly depend on the strength of the electric field.

At the material interface, after expressing $\boldsymbol {{u}_i}$ and $\boldsymbol {n}$ in terms of ${\phi }_i$ and ${\eta }$, the kinematic condition (2.29) becomes

(3.5)\begin{equation} \frac{\partial{\eta}}{\partial {t}} - \frac{\partial {\phi}_i}{\partial {x}}\frac{\partial {\eta}}{\partial {x}} + \frac{\partial {\phi}_i}{\partial {z}} = 0 \quad \mathrm{at} \ {z}={\eta} \quad (i=1,2).\end{equation}

Similarly, (2.30) and (2.31) are rewritten in terms of ${V}_i$ and ${\eta }$ as

(3.6)\begin{gather} \sum_{i=1}^2 ({-}1)^{i-1} \left[ (1+ ({-}1)^{i-1} A_\epsilon) \left(\frac{\partial V_i}{\partial z} - \frac{\partial V_i}{\partial x} \frac{\partial \eta}{\partial x} \right) \right] = 0 \quad \mathrm{at}\ z=\eta, \end{gather}

(3.7)\begin{gather}\sum_{i=1}^2 ({-}1)^{i-1} \left[ \frac{\partial V_i}{\partial x} + \frac{\partial V_i}{\partial z} \frac{\partial \eta}{\partial x} \right] = 0 \quad \mathrm{at} \ z=\eta. \end{gather}

Due to the incompressible and irrotational properties of the fluids, the Euler equations, i.e. (2.21) and (2.22), can be integrated into the Bernoulli's equation

(3.8)\begin{align} &\sum_{i=1}^2 ({-}1)^{i-1} \left[ \frac{1+ ({-}1)^{i-1} A}{2} \left( -{g} {\eta} + \frac{\partial {\phi}_i}{\partial {t}} - \frac{1}{2}\left[ \left( \frac{\partial {\phi}_i}{\partial {x}} \right)^2 + \left( \frac{\partial {\phi}_i}{\partial {z}} \right)^2\right] \right) - {p}_i \right] \nonumber\\ &\quad = {h},\quad \mathrm{at} \ {z} = {\eta}, \end{align}

where ${h}$ is a constant that only depends on the initial conditions. After expressing $\boldsymbol {{E}_i}$ and $\boldsymbol {n}$ in terms of ${V}_i$ and ${\eta }$, (2.32) becomes

(3.9)\begin{align} &\sum_{i=1}^2 ({-}1)^{i-1} \left[ {$\displaystyle-{p}_i + \frac{{1+ ({-}1)^{i-1} A_\epsilon}}{2\left[1 + \left(\dfrac{\partial {\eta}}{\partial {x}}\right)^2\right]} E_h^2 \left( \frac{\partial {V}_i}{\partial {z}} - \frac{\partial {V}_i}{\partial {x}} \frac{\partial {\eta}}{\partial {x}} \right)^2$} \right. \nonumber\\ &\quad \left. \vphantom{\frac{{1+ ({-}1)^{i-1} A_\epsilon}}{2\left[1 + \left(\frac{\partial {\eta}}{\partial {x}}\right)^2\right]}} - \frac{{1+ ({-}1)^{i-1} A_\epsilon}}{4} E_h^2 \left( \left(\frac{\partial {V}_i}{\partial {x}}\right)^2 + \left(\frac{\partial {V}_i}{\partial {z}}\right)^2 \right) \right] = 0 \quad \mathrm{at}\ {z}={\eta}.\end{align}

Combining (3.8) and (3.9), we obtain

(3.10)\begin{align} & \sum_{i=1}^2 ({-}1)^{i-1} \left[\vphantom{\frac{1+ ({-}1) A_\epsilon}{1 + \left(\frac{\partial \eta}{\partial x}\right)^2}} (1+ ({-}1)^{i-1} A) \left({-}g \eta + \frac{\partial \phi_i}{\partial t} - \frac{1}{2}\left[ \left( \frac{\partial \phi_i}{\partial x} \right)^2 + \left( \frac{\partial \phi_i}{\partial z} \right)^2\right] \right) \right.\nonumber\\ &\quad \left. {$\displaystyle- \frac{1+ ({-}1)^{i-1} A_\epsilon}{1 \!+\! \left(\dfrac{\partial \eta}{\partial x}\right)^2} E_h^2 \left(\frac{\partial V_i}{\partial z} - \frac{\partial V_i}{\partial x} \frac{\partial \eta}{\partial x} \right)^2 \!+\! \frac{1+ ({-}1)^{i-1} A_\epsilon}{2} E_h^2 \left( \left(\frac{\partial V_i}{\partial x}\right)^2 \!+\! \left(\frac{\partial V_i}{\partial z}\right)^2\right)$} \right] \nonumber\\ &\qquad = 2h \quad \mathrm{at} \ z=\eta. \end{align}

Therefore, in terms of the velocity potentials $\phi _i$ and the voltage potentials $V_i$, our system is governed by (3.1)–(3.7) and (3.10), the forms of which are more convenient for theoretical analysis.

3.2. Unperturbed solutions

We first show the solutions for an unperturbed system ($a_0=0$), which is the base for our expansion. The governing equations for the unperturbed system are

(3.11)\begin{gather} \nabla^2 V_i^{(0)} = 0, \quad i=1,2; \end{gather}

(3.12)\begin{gather}(1+A_\epsilon) \frac{\partial V_1^{(0)}}{\partial z} - (1-A_\epsilon) \frac{\partial V_2^{(0)}}{\partial z} = 0 \quad \mathrm{at} \ z=0; \end{gather}

(3.13)\begin{gather}\frac{\partial V_1^{(0)}}{\partial x} - \frac{\partial V_2^{(0)}}{\partial x} = 0 \quad \mathrm{at}\ z=0; \end{gather}

(3.14)\begin{align} &\frac{1+A_\epsilon}{2} E_h^2 \left( \left(\frac{\partial V_1^{(0)}}{\partial x}\right)^2 - \left(\frac{\partial V_1^{(0)}}{\partial z}\right)^2 \right) - \frac{1-A_\epsilon}{2} E_h^2 \left( \left(\frac{\partial V_2^{(0)}}{\partial x}\right)^2 - \left(\frac{\partial V_2^{(0)}}{\partial z}\right)^2 \right) \nonumber\\ &\quad = 2h \quad \mathrm{at} \ z=0. \end{align}

These equations give the unperturbed voltage potentials $V_1^{(0)}$ and $V_2^{(0)}$,

(3.15)\begin{equation} V_1^{(0)}(t,x,z) = V_2^{(0)}(t,x,z) ={-} x. \end{equation}

Due to $a_0=0$, the governing equations give a trivial solution for the velocity potentials, namely $\phi _i^{(0)}=0$ $(i=1,2)$.

3.3. General expansion solutions

For a material interface with a small initial perturbation, i.e. $a_0 \ll 1$, we expand all physical quantities in terms of powers of $a_0$. This means that, for a function $f$, we expand it as

(3.16)\begin{equation} f=\sum_n f^{(n)}, \end{equation}

where the term $f^{(n)}$ is proportional to $a_0^n$. In particular, the material interface is expanded as

(3.17)\begin{equation} \eta(t,x) = \sum_{n=1}^\infty \eta^{(n)}(t,x). \end{equation}

The derivation for the perturbation expansion procedure is given in Appendix A. Here we show the general perturbation solutions up to any order.

The governing equations for the $n$th-order ($n \geqslant 1$) perturbation solutions are listed in Appendix A (see (A7)–(A13)). In terms of perturbation expansion, the solutions to these equations are in the form of

(3.18)\begin{gather} \eta^{(n)}(t,x) = \sum_{0 \leqslant j \leqslant n} a_{j}^{(n)}(t) \cos(jx) , \end{gather}

(3.19)\begin{gather}\phi_1^{(n)}(t,x,z) = \sum_{0 \leqslant j \leqslant n} b_{j}^{(n)}(t) \cos(jx) {\rm e}^{- jz} , \end{gather}

(3.20)\begin{gather}\phi_2^{(n)}(t,x,z) = \sum_{0 \leqslant j \leqslant n} \tilde{b}_{j}^{(n)}(t) \cos(jx) {\rm e}^{ jz} , \end{gather}

(3.21)\begin{gather}V_1^{(n)}(t,x,z) = \sum_{0 \leqslant j \leqslant n} d_{j}^{\prime (n)}(t) \sin(jx) {\rm e}^{- jz} , \end{gather}

(3.22)\begin{gather}V_2^{(n)}(t,x,z) = \sum_{0 \leqslant j \leqslant n} \tilde{d}_{j}^{\prime (n)}(t) \sin(jx) {\rm e}^{ jz} . \end{gather}

We define the following notations that will appear frequently in our expansion:

(3.23)\begin{gather} W^{-} ={-} A_\epsilon^2 E_h^2, \end{gather}

(3.24)\begin{gather}W^{+} = A_\epsilon^3 E_h^2, \end{gather}

(3.25)\begin{gather}\sigma_j = \sqrt{j \left(A g + j W^{-} \right)} = \sqrt{j \left(A g - j A_\epsilon^2 E_h^2 \right)} . \end{gather}

The coefficient $a_{j}^{(n)}(t)$ in (3.18) for the interface $\eta$ is

(3.26)\begin{equation} a_{j}^{(n)}(t) = \frac{{\rm e}^{\sigma_j t}}{2 \sigma_j } \int_0^t F_{j}^{(n)}(\tau) {\rm e}^{-\sigma_j \tau} \,{\rm d} \tau - \frac{{\rm e}^{-\sigma_j t}}{2 \sigma_j } \int_0^t F_{j}^{(n)}(\tau) {\rm e}^{\sigma_j \tau} \,{\rm d} \tau + a_0 \cosh(\sigma_1 t) \delta_{1j} \delta_{1n} , \end{equation}

which is the solution for the governing equation

(3.27)\begin{equation} \frac{{\rm d}^2 a_{j}^{(n)}}{{\rm d}t^2} - \sigma_j^2 a_{j}^{(n)} = F_{j}^{(n)}, \end{equation}

and $\sigma _j$ given by (3.25) is the eigenvalue of this equation. The source term $F_{j}^{(n)}$ is

(3.28)\begin{align} F_{j}^{(n)} &= \frac{1+A}{2} \frac{{\rm d} R_{j}^{(n)}}{{\rm d}t} + \frac{1-A}{2} \frac{{\rm d} \tilde{R}_{j}^{(n)}}{{\rm d}t} \nonumber\\ &\quad + j \left[ U_{j}^{(n)} - E_h^2 A_\epsilon S_{j}^{\prime (n)} + \frac{1-A_\epsilon^2}{2} E_h^2 T_{j}^{(n)} \right]. \end{align}

The explicit expressions for $R_{j}^{(n)}$, $\tilde {R}_{j}^{(n)}$, $S_{j}^{\prime (n)}$, $T_{j}^{(n)}$ and $U_{j}^{(n)}$ can be determined from (A16)–(A20) given in Appendix A. We comment that the right-hand side of (3.28), namely $R_{j}^{(n)}$, $\tilde {R}_{j}^{(n)}$, $S_{j}^{\prime (n)}$, $T_{j}^{(n)}$ and $U_{j}^{(n)}$, only contains solutions of order lower than $n$ (see Appendix A). Therefore, starting from the first order, $a_{j}^{(n)}$ can be solved from (3.27) order by order.

The coefficients $b_{j}^{(n)}(t)$ in (3.19) and $\tilde {b}_{j}^{(n)}(t)$ in (3.20) for fluid velocity potentials can then be determined from $a_{j}^{(n)}(t)$, i.e.

(3.29)\begin{gather} b_{j}^{(n)}(t) = \frac{1}{j} \left( \frac{{\rm d} a_{j}^{(n)}(t)}{{\rm d}t} - R_{j}^{(n)}(t)\right) , \end{gather}

(3.30)\begin{gather}\tilde{b}_{j}^{(n)}(t) ={-}\frac{1}{j} \left( \frac{{\rm d} a_{j}^{(n)}(t)}{{\rm d}t} - \tilde{R}_{j}^{(n)}(t) \right) . \end{gather}

The coefficients $d_{j}^{\prime (n)}(t)$ in (3.21) and $\tilde {d}_{j}^{\prime (n)}(t)$ in (3.22) for voltage potentials can also be determined from $a_{j}^{(n)}(t)$. The results are

(3.31)\begin{gather} d_{j}^{\prime (n)}(t) ={-} \left[ A_\epsilon a_{j}^{(n)}(t) + \frac{1}{j}\left( S_{j}^{\prime (n)}(t) - \frac{1-A_\epsilon}{2} T_{j}^{(n)}(t) \right)\right] , \end{gather}

(3.32)\begin{gather}\tilde{d}_{j}^{\prime (n)}(t) ={-} \left[ A_\epsilon a_{j}^{(n)}(t) + \frac{1}{j}\left( S_{j}^{\prime (n)}(t) + \frac{1+A_\epsilon}{2} T_{j}^{(n)}(t) \right) \right] . \end{gather}

The detailed derivation for the perturbation solutions shown in this subsection can be found in Appendix B.

3.4. Explicit expressions for perturbation solutions

Equations (3.26)–(3.32) provide us with a general procedure to obtain perturbed solutions up to any order. Here we show the expressions for the perturbed interface up to the third order in terms of $a_0$. As we will show later, to give an accurate prediction for the perturbed interface, one must include all terms up to the third order.

The shape of the interface at time $t$, up to the order of $a_0^3$, can be described by

(3.33)\begin{equation} \eta(t, x) = \eta^{(1)}(t, x) + \eta^{(2)}(t, x) + \eta^{(3)}(t, x) + O(a_0^4), \end{equation}

where the $n$th-order component $\eta ^{(n)} (t, x)$ $(n=1,2,3)$ are given below.

The first-order component is

(3.34)\begin{equation} \eta^{(1)}(t, x) = a_{1}^{(1)}(t) \cos(x), \end{equation}

with

(3.35)\begin{equation} a_{1}^{(1)}(t) = a_0 \cosh(\sigma_1 t). \end{equation}

The second-order component is

(3.36)\begin{equation} \eta^{(2)}(t, x) = a_{2}^{(2)}(t) \cos(2x), \end{equation}

with

(3.37)\begin{equation} a_{2}^{(2)}(t) = a_0^2 \left[ \left( c_1 + c_2 \right) \cosh(\sigma_2 t) - c_1 - c_2 \cosh(2 \sigma_1 t) \right]. \end{equation}

The third-order component is

(3.38)\begin{equation} \eta^{(3)}(t, x) = a_{1}^{(3)}(t) \cos(x) + a_{3}^{(3)}(t) \cos(3x), \end{equation}

with

(3.39)\begin{align} a_{1}^{(3)}(t) &= a_0^3 \left[ c_3 t \sinh(\sigma_1 t) + \left( c_4 + c_5 + c_6 \right) \cosh(\sigma_1 t) \right. \nonumber\\ &\quad \left. - c_4 \cosh(3 \sigma_1 t) - c_5 \cosh((\sigma_1-\sigma_2) t) - c_6 \cosh((\sigma_1+\sigma_2) t) \right] , \end{align}

(3.40)\begin{align} a_{3}^{(3)}(t) &= a_0^3 \left[ \left( c_7 + c_8 + c_9 + c_{10} \right) \cosh(\sigma_3 t) - c_7 \cosh(\sigma_1 t) - c_8 \cosh(3 \sigma_1 t) \right.\nonumber\\ &\quad \left. - c_9 \cosh((\sigma_1-\sigma_2) t) - c_{10} \cosh((\sigma_1+\sigma_2) t) \right]. \end{align}

Here the constants $c_m$ ($m=1,2,\ldots,10$) in (3.37), (3.39) and (3.40) are given by the following expressions:

(3.41)\begin{gather} c_1 = \frac{1}{2\sigma_2^2}(A \sigma_1^2 - W^{+}), \end{gather}

(3.42)\begin{gather}c_2 ={-}\frac{1}{2(\sigma_2^2 - 4\sigma_1^2)}(A \sigma_1^2 + W^{+}), \end{gather}

(3.43)\begin{gather}c_3 = \frac{1}{2\sigma_1}\left[ -\frac{1}{8} \sigma_1^2 + A \sigma_1^2 \left( \frac{c_2}{2} - c_1 \right) + W^{+} \left( \frac{c_2}{2} + c_1 \right) - \frac{3}{8} W^{-} \right], \end{gather}

(3.44)\begin{gather}c_4 ={-}\frac{1}{8\sigma_1^2}\left[ -\frac{1}{8} \sigma_1^2 - \frac{3 c_2}{2}A \sigma_1^2 + \frac{c_2}{2} W^{+} - \frac{1}{8} W^{-} \right], \end{gather}

(3.45)\begin{gather}c_5 = \frac{1}{2\sigma_2 (2\sigma_1 - \sigma_2)}\left[ A \sigma_1 (\sigma_1 - \sigma_2) (c_1 + c_2) - W^{+} (c_1 + c_2) \right], \end{gather}

(3.46)\begin{gather}c_6 ={-} \frac{1}{2\sigma_2 (2\sigma_1 + \sigma_2)}\left[ A \sigma_1 (\sigma_1 + \sigma_2) (c_1 + c_2) - W^{+} (c_1 + c_2) \right], \end{gather}

(3.47)\begin{gather}c_7 = \frac{3}{\sigma_3^2 - \sigma_1^2}\left[ \frac{1}{16}\sigma_1^2 - A \sigma_1^2 c_2 + W^{+} (2c_1 + c_2) + \frac{3}{16} W^{-} \right], \end{gather}

(3.48)\begin{gather}c_8 = \frac{3}{\sigma_3^2 - 9\sigma_1^2}\left[ -\frac{1}{16}\sigma_1^2 + A \sigma_1^2 c_2 + W^{+} c_2 + \frac{1}{16} W^{-} \right], \end{gather}

(3.49)\begin{gather}c_9 = \frac{3}{\sigma_3^2 - (\sigma_1-\sigma_2)^2}\left[ \frac{1}{2} A \sigma_1 \sigma_2 (c_1 + c_2) - W^{+} (c_1 + c_2) \right], \end{gather}

(3.50)\begin{gather}c_{10} = \frac{3}{\sigma_3^2 - (\sigma_1+\sigma_2)^2}\left[ -\frac{1}{2} A \sigma_1 \sigma_2 (c_1 + c_2) - W^{+} (c_1 + c_2) \right]. \end{gather}

The terms that are not displayed in (3.36) and (3.38) vanish. More specifically,

(3.51)\begin{equation} a_{1}^{(2)}(t) = a_{2}^{(3)}(t) = 0. \end{equation}

The derivation of the governing equations and the corresponding perturbation solutions for the first, second and third orders are given in Appendices C, D and E, respectively.

Taking the derivative of (3.33) with respect to $t$ gives the explicit expression for the velocity of an arbitrary point on the material interface up to the order of $a_0^3$,

(3.52)\begin{equation} \dot{\eta}(t, x) = \dot{\eta}^{(1)}(t, x) + \dot{\eta}^{(2)}(t, x) + \dot{\eta}^{(3)}(t, x) + O(a_0^4), \end{equation}

where the three leading-order components $\dot {\eta }^{(n)} (t, x)$ $(n=1,2,3)$ are as follows.

The first-order component is

(3.53)\begin{equation} \dot{\eta}^{(1)}(t, x) = \dot{a}_{1}^{(1)}(t) \cos(x), \end{equation}

with

(3.54)\begin{equation} \dot{a}_{1}^{(1)}(t) = a_0 \sigma_1 \sinh(\sigma_1 t). \end{equation}

The second-order component is

(3.55)\begin{equation} \dot{\eta}^{(2)}(t, x) = \dot{a}_{2}^{(2)}(t) \cos(2x), \end{equation}

with

(3.56)\begin{equation} \dot{a}_{2}^{(2)}(t) = a_0^2 \left[ \sigma_2 \left( c_1 + c_2 \right) \sinh(\sigma_2 t) - 2 \sigma_1 c_2 \sinh(2 \sigma_1 t) \right] . \end{equation}

The third-order component is

(3.57)\begin{equation} \dot{\eta}^{(3)}(t, x) = \dot{a}_{1}^{(3)}(t) \cos(x) + \dot{a}_{3}^{(3)}(t) \cos(3x), \end{equation}

with

(3.58)\begin{align} \dot{a}_{1}^{(3)}(t) &= a_0^3 \left[ c_3 \sinh(\sigma_1 t) + \sigma_1 c_3 t \cosh(\sigma_1 t) + \sigma_1 \left( c_4 + c_5 + c_6 \right) \sinh(\sigma_1 t) \right.\nonumber\\ &\quad - 3 \sigma_1 c_4 \sinh(3 \sigma_1 t) - (\sigma_1-\sigma_2) c_5 \sinh((\sigma_1-\sigma_2) t)\nonumber\\ &\quad \left.- (\sigma_1+\sigma_2)c_6 \sinh((\sigma_1+\sigma_2) t) \right] , \end{align}

(3.59)\begin{align} \dot{a}_{3}^{(3)}(t) &= a_0^3 \left[ \sigma_3 \left( c_7 + c_8 + c_9 + c_{10} \right) \sinh(\sigma_3 t) - \sigma_1 c_7 \sinh(\sigma_1 t) \right.\nonumber\\ &\quad - 3 \sigma_1 c_8 \sinh(3 \sigma_1 t) - (\sigma_1-\sigma_2) c_9 \sinh((\sigma_1-\sigma_2) t)\nonumber\\ &\quad \left. - (\sigma_1+\sigma_2) c_{10} \sinh((\sigma_1+\sigma_2) t)\right] . \end{align}

Here the constants $c_m$ ($m=1,2,\ldots,10$) in (3.56), (3.58) and (3.59) are given by (3.41)–(3.50).

5.1. Numerical method

To conduct validation studies, we need the numerical solution of the governing equations given by (2.21)–(2.26) with the boundary conditions given by (2.27)–(2.32). Solving these equations directly requires performing numerical computations in two dimensions. However, it is the evolution of the material interface that we are interested in, not the solution in the interior of the fluids. For our system, the evolution of the material interface can be formulated in terms of vortex dynamics that only needs to solve a one-dimensional vortex sheet that coincides with the material interface, and no computation in the interior regions is needed. This is possible because the fluids are inviscid, incompressible and perfect dielectric. Therefore, there is no vorticity or bulk charge in the interior regions of the fluids. We will take the vortex sheet method in our numerical simulations.

The vortex sheet method is one of the numerical methods designed to tackle the problems that involve multi-phase fluids (Moore Reference Moore1981; Baker, Meiron & Orszag Reference Baker, Meiron and Orszag1982; Shin, Sohn & Hwang Reference Shin, Sohn and Hwang2018). This method has been successfully adapted in the study of the classical RTI (Baker, Meiron & Orszag Reference Baker, Meiron and Orszag1980; Kerr Reference Kerr1988; Sohn Reference Sohn2004; Shin, Sohn & Hwang Reference Shin, Sohn and Hwang2022), which does not involve electric fields. For inviscid and incompressible fluids, based on the Biot–Savart law, the vortex sheet method formulates the dynamics of the entire system into the evolution of vorticity confined to the material interfaces between the fluids (Baker et al. Reference Baker, Meiron and Orszag1980; Moore Reference Moore1981; Kerr Reference Kerr1988; Sohn Reference Sohn2004). Then one can determine the evolution of the material interface by solving a boundary integral equation at the material interface.

Our system is more complicated. It involves two inviscid, incompressible and perfect dielectric fluids in the presence of both gravitational force and external electric fields. Even so, we are still able to extend the vortex method to our system. Such an extension can be achieved due to the fact that the fluid velocity field and the electric displacement field in our system share the following important properties: Both of them are rotation free and can be expressed in terms of potentials; they are both divergence free and governed by the Laplace equation in the interior of fluids; both fields are continuous in the normal direction at the material interface. Therefore, similar to the fluid velocity field, the electric displacement field in our system can also be formulated in terms of a ‘vortex sheet’ and determined from a boundary integral equation. In the case of the classical RTI, one only needs to solve the boundary integral equation for fluid vorticity. In our case, one needs to solve two coupled boundary integral equations, one for fluid vorticity and the other for ‘vorticity’ in the electric displacement field. Although two vortex sheets exist in our system (one for the fluid velocity and the other for the electric field), both vortex sheets are confined to the same material interface and, therefore, have the same location. In other words, all physical quantities are only evaluated and evolved at the one-dimensional material interface during the computation.

Here we provide the key steps for the vortex dynamics formulation of our system. The vortex strength, which is the jump in the dimensionless tangential velocity across the material interface, is defined as

(5.1)\begin{equation} \gamma = (\boldsymbol{u_2} - \boldsymbol{u_1}) \boldsymbol{\cdot} \boldsymbol{s}. \end{equation}

Since the velocities of the two fluids are not identical at the material interface ($\boldsymbol {u_1} \ne \boldsymbol {u_2}$), the motion of the material interface follows the average velocity

(5.2)\begin{equation} \bar{\boldsymbol{U}} = (\boldsymbol{u_1} + \boldsymbol{u_2}) /2. \end{equation}

Due to the incompressible and irrotational properties of the fluids, i.e. (2.23) and (2.24), the velocity of the material interface $\bar {\boldsymbol {U}}$ can be expressed in terms of the vortex strength $\gamma$ based on the Biot–Savart law,

(5.3)\begin{equation} \bar{\boldsymbol{U}}(s, t)=\frac{1}{2 {\rm \pi}} \ \text{p.v.} \int_{L} \frac{\tilde{\boldsymbol{k}} \times[{\boldsymbol{X}}(s, t)-{\boldsymbol{X}}(\tilde{s}, t)]}{|{\boldsymbol{X}}(s, t)-{\boldsymbol{X}}(\tilde{s}, t)|^{2}} \gamma(\tilde{s}, t) \,{\rm d} \tilde{s} . \end{equation}

Here $L$ denotes the whole material interface that contains an infinite number of periods, $\text {p.v.}$ denotes the Cauchy principal value of an integral, $s$ denotes the dimensionless arclength coordinate at the material interface and $\tilde {\boldsymbol {k}}$ is the unit vector perpendicular to the two-dimensional plane of the fluid system. Since the horizontal boundaries of the system are periodic, (5.3) can be simplified as

(5.4)\begin{equation} \bar{U}_x - {\rm i} \bar{U}_z = \frac{1}{4{\rm \pi} {\rm i}} \ \text{p.v.} \int_{l} \cot \left( \frac{\xi-\tilde{\xi}}{2} \right) \tilde{\gamma} \, {\rm d} \tilde{s}, \end{equation}

where $\xi = x+{\rm i}z$ is the complex notation of the interface location, $\bar {U}_x$ and $\bar {U}_z$ are the horizontal and vertical components of $\bar {\boldsymbol {U}}$, $l$ is one period of the material interface (due to the property of periodicity, it can be any period), and the notation $\widetilde {({\cdot })} \equiv ({\cdot })(\tilde {s})$ denotes the dummy variables in the integration along the interface. We comment that the velocity field given by (5.4) satisfies the far-field boundary condition given by (2.27).

The electric displacement field in our system can be formulated similarly. For the electric displacement field, we define the dimensionless electric displacement in fluid $1$ by

(5.5)\begin{equation} \boldsymbol{D_1} = \hat{\epsilon}_1\boldsymbol{\hat{E}_1}/({\hat{\epsilon}_1+\hat{\epsilon}_2})\hat{E}_h = (1+A_\epsilon) \boldsymbol{{E}_1}/2, \end{equation}

and the dimensionless electric displacement in fluid $2$ by

(5.6)\begin{equation} \boldsymbol{D_2} = \hat{\epsilon}_2\boldsymbol{\hat{E}_2}/({\hat{\epsilon}_1+\hat{\epsilon}_2})\hat{E}_h = (1-A_\epsilon) \boldsymbol{{E}_2}/2. \end{equation}

Similar to the velocity field, we define the ‘vortex strength’ for the dimensionless electric displacement field by

(5.7)\begin{equation} \beta = (\boldsymbol{D_2} - \boldsymbol{D_1}) \boldsymbol{\cdot} \boldsymbol{s}, \end{equation}

which is the jump in the tangential component of the dimensionless electric displacement field at the material interface. The electric displacement field at the material interface is defined by the average electric displacement of the two fluids

(5.8)\begin{equation} \bar{\boldsymbol{D}} = (\boldsymbol{D_1} + \boldsymbol{D_2}) /2. \end{equation}

Due to (2.25) and (2.26) and the Biot–Savart law, we have

(5.9)\begin{equation} \bar{\boldsymbol{D}}(s, t)=\frac{1}{2 {\rm \pi}} \ \text{p.v.} \int_{L} \frac{\tilde{\boldsymbol{k}} \times[{\boldsymbol{X}}(s, t)-{\boldsymbol{X}}(\tilde{s}, t)]}{|{\boldsymbol{X}}(s, t)-{\boldsymbol{X}}(\tilde{s}, t)|^{2}} \beta(\tilde{s}, t) \,{\rm d} \tilde{s} + \boldsymbol{D^{ext}}, \end{equation}

where $\boldsymbol {D^{ext}}$ is the contribution of the external electric field. Based on (5.5) and (5.6) and the boundary condition given by (2.28), the horizontal component of $\boldsymbol {D^{ext}}$ is ${D^{ext}_x} = \frac {1}{2}[{(1+A_\epsilon )}/{2} + ({1-A_\epsilon )}/{2}] = \frac {1}{2}$, namely the average electric displacement due to the external electric field, and the vertical component is $D^{ext}_z = 0$. Due to the periodicity of our system in the horizontal direction, (5.9) can then be simplified as

(5.10)\begin{equation} \bar{D}_x - {\rm i} \bar{D}_z = \frac{1}{4{\rm \pi} {\rm i}} \ \text{p.v.} \int_{l} \cot \left( \frac{\xi-\tilde{\xi}}{2} \right) \tilde{\beta} \, {\rm d} \tilde{s} + \frac{1}{2}, \end{equation}

where $\bar {D}_x$ and $\bar {D}_z$ are the horizontal and vertical components of $\bar {\boldsymbol {D}}$, respectively.

Based on the vortex dynamics formulation for the velocity and electric displacement fields, the evolution of the material interface can then be represented by the following one-dimensional integral equations along the material interface in terms of dimensionless variables:

(5.11)\begin{gather} \frac{{\rm d}\kern 0.06em {\boldsymbol{X}}}{{\rm d}t} = \bar{\boldsymbol{U}}; \end{gather}

(5.12)\begin{gather}\frac{{\rm d} \varGamma}{{\rm d}t}= 2 A \left.\int_{s^-}^{s^+}\frac{{\rm d} \bar{\boldsymbol{U}}}{{\rm d}t} \boldsymbol{\cdot} {\boldsymbol{s}} \, {\rm d}s + \left( 2 A g z - F_\gamma - F_e \right) \right|_{s^-}^{s^+} \quad \textrm{for any segment}\;(s^-, s^+); \end{gather}

(5.13)\begin{gather}\beta + 2\, A_\epsilon \bar{\boldsymbol{D}} \boldsymbol{\cdot} {\boldsymbol{s}} = 0; \end{gather}

(5.14)\begin{gather}\bar{U}_x (x,z) ={-}\frac{1}{4 {\rm \pi}}\, \text{p.v.} \int_{l} \frac{\sinh (z-\tilde{z})}{\cosh (z-\tilde{z})-\cos (x-\tilde{x})} \tilde{\gamma} \,{\rm d} \tilde{s}, \end{gather}

(5.15)\begin{gather}\bar{U}_z (x,z) = \frac{1}{4 {\rm \pi}} \, \text{p.v.} \int_{l} \frac{\sin (x-\tilde{x})}{\cosh (z-\tilde{z})-\cos (x-\tilde{x})} \tilde{\gamma} \,{\rm d} \tilde{s}; \end{gather}

(5.16)\begin{gather}\bar{D}_x (x,z) ={-}\frac{1}{4 {\rm \pi}} \, \text{p.v.} \int_{l} \frac{\sinh (z-\tilde{z})}{\cosh (z-\tilde{z})-\cos (x- \tilde{x})} \tilde{\beta} \,{\rm d} \tilde{s} + \frac{1}{2}, \end{gather}

(5.17)\begin{gather}\bar{D}_z (x,z) = \frac{1}{4 {\rm \pi}} \, \text{p.v.} \int_{l} \frac{\sin (x-\tilde{x})}{\cosh (z-\tilde{z})-\cos (x-\tilde{x})} \tilde{\beta} \,{\rm d} \tilde{s}, \end{gather}

where

(5.18)\begin{gather} \varGamma = \int_{s^-}^{s^+} \gamma(s,t) \, {\rm d}s, \end{gather}

(5.19)\begin{gather}F_\gamma ={-} \frac{A}{4}\gamma^2, \end{gather}

(5.20)\begin{gather}F_e = \frac{4 A_\epsilon}{1-A_\epsilon^2} E_h^2 \left( \bar{D}_x^2 + \bar{D}_y^2 - \frac{\beta^2}{4} \right), \end{gather}

In (5.11), $\boldsymbol {X} = (x, z)$ is the location of the material interface. In (5.12) and (5.18), $s^-$ and $s^+$ denote the dimensionless arclength coordinates of the starting and ending points of an arbitrary segment at the material interface, and $\varGamma$ is the dimensionless vorticity residing in the segment. The Atwood number for fluid densities $A$ and that for electric permittivities $A_\epsilon$ in (5.12), (5.13), (5.19) and (5.20) are defined in (2.13) and (2.14).

We comment that the formulation given by (2.21)–(2.32) and the vortex sheet formulation given by (5.11)–(5.17) are equivalent. More specifically, (5.14) and (5.15) given by the Biot–Savart law for the velocity field satisfy (2.23), (2.24), (2.27) and (2.29); (5.16) and (5.17) given by the Biot–Savart law for the electric displacement field satisfy (2.25), (2.26), (2.28) and (2.30). Equation (5.13) corresponds to the boundary condition given by (2.31). Equation (5.12) can be obtained by taking the difference between (2.21) and (2.22) and simplifying the resulting equation by the boundary condition (2.32). Equation (5.11) simply comes from the definition. Therefore, although the two formulations are different, they both solve exactly the same system and have the same solution. There is no simplification between these two formulations. However, conducting numerical simulations based on the vortex sheet method has the following three distinct advantages for our system. (1) The vortex sheet method only needs to perform one-dimensional computation along the material interface, rather than two-dimensional computation based on solving (2.21)–(2.32). This leads to a saving in computing time. (2) The key issue for the RTI is to determine how fast the fingers at the material interface grow with time. Therefore, an accurate determination of the shape and location of the material interface is very important. In the vortex sheet method, its data structure directly represents the shape and location of the material interface. On the other hand, the usual numerical methods for solving (2.21)–(2.32) in the two-dimensional physical space need to construct or estimate these quantities from the two-dimensional grid-based data. Such constructions or estimations usually contain more numerical diffusion. Therefore, the vortex sheet method can compute the velocity of fingers more accurately. (3) The domain of our system is infinite in the vertical direction. However, the numerical simulation based on a two-dimensional computation of (2.21)–(2.32) can only be performed in a finite computational domain, and one must introduce the upper and lower boundaries for the finite computational domain and impose certain boundary conditions at these numerical boundaries. Such numerical boundaries are not needed in the vortex sheet method, since the computation is solely and directly performed at the material interface.

Evaluations of (5.14)–(5.17) encounter singularities when $(\tilde {x},\tilde {z})$ approaches $(x,z)$. Krasny introduced a desingularization parameter $\delta >0$ in the equations of the Biot–Savart law (Krasny Reference Krasny1986). Linear stability analysis showed that this desingularization method diminished the short wavelength instability and yielded numerically more tractable equations (Krasny Reference Krasny1986). Applying Krasny's method, we introduce the desingularization parameter $\delta$ in the Biot–Savart law for computing the average velocity and the average electric displacement at the material interface,

(5.21)\begin{gather} \bar{U}_x (x,z) ={-}\frac{1}{4 {\rm \pi}}\, \text{p.v.} \int_{l} \frac{\sinh (z-\tilde{z})}{\cosh (z-\tilde{z})-\cos (x-\tilde{x}) + \delta^2} \tilde{\gamma} \,{\rm d} \tilde{s}, \end{gather}

(5.22)\begin{gather}\bar{U}_z (x,z) = \frac{1}{4 {\rm \pi}} \, \text{p.v.} \int_{l} \frac{\sin (x-\tilde{x})}{\cosh (z-\tilde{z})-\cos (x-\tilde{x}) + \delta^2} \tilde{\gamma} \,{\rm d} \tilde{s}; \end{gather}

(5.23)\begin{gather}\bar{D}_x (x,z) ={-}\frac{1}{4 {\rm \pi}} \, \text{p.v.} \int_{l} \frac{\sinh (z-\tilde{z})}{\cosh (z-\tilde{z})-\cos (x-\tilde{x}) + \delta^2} \tilde{\beta} \,{\rm d} \tilde{s} + \frac{1}{2}, \end{gather}

(5.24)\begin{gather}\bar{D}_z (x,z) = \frac{1}{4 {\rm \pi}} \, \text{p.v.} \int_{l} \frac{\sin (x-\tilde{x})}{\cosh (z-\tilde{z})-\cos (x-\tilde{x}) + \delta^2} \tilde{\beta} \,{\rm d} \tilde{s}. \end{gather}

The introduction of the desingularization parameter $\delta ^2$ in Krasny's method removes the singularities in the Biot–Savart law, but it also reduces the numerical solution. We have used $\delta ^2=0.02$ in all our simulations. Although the error due to the desingularization parameter $\delta ^2$ is small, by applying the technique of repeated Richardson extrapolation we further reduce the error (Richardson Reference Richardson1910). More specifically, we applied the Richardson extrapolation technique twice. All results in our validation studies are obtained with this repeated Richardson extrapolation technique. We have implemented this vortex sheet method with the repeated Richardson extrapolation technique and conducted validation studies to verify the consistency and accuracy of this numerical method. Since the derivation and the implementation details of this numerical method are lengthy, the details will be presented in a separate paper focused on the numerical method.

5.2. Validation studies of the nonlinear perturbation solutions

The central issue in the study of the RTI is to predict how the overall mixing layer between the two fluids grows with time. The overall mixing layer is characterized by the vertical distance between the tip of the spike and that of the bubble. A half of this distance is known as the overall amplitude of the fingers. This is the most important and widely considered quantity for the RTI. Researchers are interested in how the overall amplitude changes with time and at what speed it changes, namely the overall velocity. Due to the symmetry in the functional form given by (3.18), all even-order terms in the perturbation solutions do not contribute to the overall amplitude or the overall velocity. Therefore, the first-order solution and the solution up to the second order give the same prediction for the overall amplitude/velocity. Similarly, the solution up to the third order and that up to the fourth order give the same prediction for the overall amplitude/velocity. From (4.8) and (4.9), we have the overall velocity

(5.25)\begin{align} v_{ov} &= \tfrac{1}{2}(v_{bb} - v_{sp}) \nonumber\\ &= v^{(1)}(t) + v^{(3)}(t) + O(a_0^5), \end{align}

where $v^{(1)}(t)$ and $v^{(3)}(t)$ are given by (4.10) and (4.12), respectively. From these two equations, (5.25) can be rewritten as

(5.26)\begin{equation} v_{ov} = \dot{a}_{1}^{(1)}(t) + \dot{a}_{1}^{(3)}(t) + \dot{a}_{3}^{(3)}(t) + O(a_0^5), \end{equation}

where $\dot {a}_{1}^{(1)}(t)$, $\dot {a}_{1}^{(3)}(t)$ and $\dot {a}_{3}^{(3)}(t)$ are given by (3.54), (3.58) and (3.59), respectively. As we stated earlier, since the second-order solution does not enter $v_{ov}$, one needs to include at least the third-order solution to exploit the nonlinear effects of the overall velocity.

Now we compare the predictions for the overall velocity from the first-order solution (namely the linear theory) and those from the solution up to the fourth order given by (5.26) with the data from numerical simulations. Our system involves four fundamental physical parameters: the dimensionless initial amplitude $a_0$, the fraction of the driving forces due to the horizontal electric field $E_h^2$, the dielectric Atwood number $A_\epsilon$ and the Atwood number $A$. This is a four-dimensional phase space. We will set a base point in the phase space, and vary one parameter at a time in the phase space to conduct our validation studies. Since our solution is based on the perturbation expansion of the dimensionless initial amplitude $a_0$, $a_0$ should be small. We choose $a_0 = 0.1$ as the base point. In our system, two competing forces act on the material interface, the strengths of which are represented by the fraction of the driving forces due to gravity $g$ and the fraction of the driving forces due to the external horizontal electric field ${E_h}^2$. Since $g + {E_h}^2 = 1$, we choose $g = {E_h}^2 = 1/2$ as our base point. Both $|A|$ and $|A_{\epsilon }|$ vary from 0 to 1. Therefore, we choose $A = 0.5$ and $A_{\epsilon } = 0.5$ as our base point. In figures 2–5 we plot the normalized overall velocity of fingers $v_{ov}/a_0$ versus the dimensionless time $t$. In each figure we vary one of the four physical parameters aforementioned, namely $a_0$ in figure 2, ${E_h}^2$ in figure 3, $A_\epsilon$ in figure 4 and $A$ in figure 5. If we kept all three remaining parameters at the base point in each figure, there would be repeated subfigures among figures 2–5. To avoid such repetitions and maximize the information extracted from the phase space, we will choose one or two of the remaining parameters to deviate away from the base point. The rest of the remaining parameters are located at the base point.

Figure 2. Comparison between the results from the nonlinear perturbation solution given by (5.26) (blue solid curve), the linear theory (black solid curve) and the numerical simulations (blue ‘$\circ$’) for the dimensionless overall velocity of fingers $v_{ov}/{a_0}$ with various dimensionless initial amplitudes $a_0 = 0.1$, 0.2, 0.4, 0.8. The other physical parameters are listed in the main text.

Figure 3. Comparison between the results from the nonlinear perturbation solution given by (5.26) (blue solid curve), the linear theory (black solid curve) and the numerical simulations (blue ‘$\circ$’) for the dimensionless overall velocity of fingers $v_{ov}/{a_0}$ with various strengths of horizontal electric fields applied ($E_h^2 = 1/{65}, 1/{17}, 1/{5}$ and $1/{2}$). The other physical parameters are listed in the main text.

Figure 4. Comparison between the results from the nonlinear perturbation solution given by (5.26) (blue solid curve), the linear theory (black solid curve) and the numerical simulations (blue ‘$\circ$’) for the dimensionless overall velocity of fingers $v_{ov}/{a_0}$ with various dielectric Atwood numbers $A_\epsilon = -0.1$, $-0.3$, $-0.5$, $-0.7$. The other physical parameters are listed in the main text.

Figure 5. Comparison between the results from the nonlinear perturbation solution given by (5.26) (blue solid curve), the linear theory (black solid curve) and the numerical simulations (blue ‘$\circ$’) for the dimensionless overall velocity of fingers $v_{ov}/{a_0}$ with various Atwood numbers $A = 0.3$, $0.5$, $0.7$, $1.0$. The other physical parameters are listed in the main text.

In figure 2 we set $g = {E_h}^2 = 1/2$, $A = 0.6$ and $A_{\epsilon } = 0.5$ as the constant parameters, and vary $a_0$ from $0.1$ to $0.8$. For each subsequent subfigure, we double the dimensionless amplitude of the initial perturbation. In figure 2(a)–2(d), $a_0$ is $0.1, 0.2, 0.4$ and $0.8$, respectively. Figure 2 shows that the contributions from the nonlinear terms are important. By including the nonlinear correction terms, the nonlinear perturbation solution provides predictions that agree with the numerical data better than the linear theory and, thus, provides a larger range of validity for theoretical prediction. Particularly, in figure 2(d) where the dimensionless amplitude of the initial perturbation is large ($a_0=0.8$), the nonlinear effects are important even at early times. In this case, the prediction from the linear theory deviates from the numerical results early ($t \approx 0.6$), while the nonlinear perturbation solution provides a reasonable prediction up to three times longer time ($t \approx 2$). Figure 2 shows that the larger $a_0$ is, the earlier the nonlinear contribution occurs, as one expected.

In figure 3 we vary the dimensionless strength of the external horizontal electric field $E_h^2$ and plot the normalized overall velocity of fingers versus the dimensionless time. In figures 3(a)–3(d) we begin with a weak horizontal electric field and double its strength in each subsequent case until the horizontal electric field suppresses the interfacial instability. More precisely, we vary the ratio of the gravitational force to the electric-field force $g/E_h^2 = 64$, $16$, $4$ and $1$. These values correspond to the dimensionless electric-field strength $E_h^2$ being $1/{65}$, $1/{17}$, $1/{5}$ and $1/{2}$, respectively. This covers the cases from a situation in which the material interface is dominated by gravity-induced instability to that in which the interface is close to the stable state. The other parameters are $a_0 = 0.1$, $A = 0.4$ and $A_{\epsilon } =0.6$. Figure 3 shows that the regime of validity of the nonlinear perturbation solution is larger than that of the linear theory in all cases. Particularly, when $E_h^2=1/{2}$, the nonlinear perturbation solution provides much more accurate predictions than the linear theory. Figure 3 also shows that both the nonlinear perturbation solution and the linear theory are valid over a larger temporal domain as the strength of the horizontal electric field $E_h^2$ increases. This is because as $E_h^2$ increases, the material interface becomes less unstable. Therefore, the finger amplitude grows at a slower rate, which extends the regime of validity of both theories over a longer period of time.

In figure 4 we vary the dielectric Atwood number $A_\epsilon$. Since $A_\epsilon$ lies in the range $[ -1, 1 ]$ and we have presented the results for $A_\epsilon = 0.5$ and $0.6$ in figures 2 and 3, to explore more points in the phase space, we present the results for $A_\epsilon = -0.1$, $-0.3$, $-0.5$ and $-0.7$ in figures 4(a)–4(d), respectively. The other parameters are $a_0=0.2$, $A=0.5$ and $g = {E_h}^2 =1/{2}$. Figure 4 shows that the predictions of the nonlinear perturbation solution are valid in a larger regime than those of the linear theory in all cases. Particularly, when $A_\epsilon =-0.7$, the prediction of the linear theory deviates significantly from the numerical results as early as $t \approx 2$ (see figure 4d). The importance of the nonlinear behaviour is clearly shown by the numerical results and captured by the nonlinear perturbation solution in figure 4(d). As $A_\epsilon$ changes from $-0.1$ the $-0.7$, the material interface becomes less unstable, and the regime of validity of the nonlinear perturbation solution becomes larger consequently. For $A_\epsilon \leqslant -0.8$, the effects of the driving force due to the horizontal electric field become dominant such that the interface becomes stable.

In figure 5 we vary the Atwood number $A$ and present the validation study for $A = 0.3$, $0.5, 0.7$ and $1.0$ in figures 5(a)–5(d), respectively. The other parameters are $a_0=0.1$, $A_\epsilon =-0.5$ and $g = {E_h}^2 =1/{2}$. Similar to the previous comparisons, figure 5 also shows that the predictions of the nonlinear perturbation solution are valid in a larger regime than that of the linear theory. Since the increase of $A$ destabilizes the system, both the regime of validity of the nonlinear perturbation solution and that of the linear theory become smaller as $A$ increases. For $0 \leqslant A \leqslant 0.2$, the destabilizing effects from gravity are weak and the stabilizing effects from the horizontal electric field dominate, such that the material interface becomes stable.

The nonlinear theoretical prediction for the overall velocity given by (5.26) only contains the contributions up to $O(a_0^4)$. For an unstable material interface, the contributions from the terms of higher orders become important as time progresses. Then the predictions of the nonlinear perturbation solution given by (5.26) will also start to deviate from the numerical solutions. Such a deviation is a complicated function of the four parameters $a_0$, $E_h^2$, $A$ and $A_\epsilon$. Depending on where the system sits in the four-dimensional phase space, the deviation can be an overestimate or an underestimate. For example, in figure 3(c) the nonlinear perturbation solution starts to underestimate the overall velocity around $t=6$, while in figure 3(d) the nonlinear perturbation solution starts to overestimate the overall velocity around $t=18$. To verify that our solution for the overall velocity up to the fourth order, namely the terms $\dot {a}_1^{(1)}$, $\dot {a}_1^{(3)}$ and $\dot {a}_3^{(3)}$ in (5.26), is correct even in the time range in which the theoretical predictions deviate from the numerical results, we need a method for extracting the values of these terms from the data of numerical simulations. Then it allows us to compare the extracted numerical results with the theoretical predictions for $\dot {a}_1^{(1)}$, $\dot {a}_1^{(3)}$ and $\dot {a}_3^{(3)}$ given by (3.54), (3.58) and (3.59). A good agreement between these two will not only prove the correctness of our theoretical formula given by (5.26) but also confirm that the deviation between the numerical results and the theoretical predictions is indeed due to the contributions from the orders higher than those in (5.26). We achieve this by the following method.

From (3.17) and (3.18) the material interface can also be expressed in terms of Fourier modes

(5.27)\begin{align} \eta(t,x)&= \sum_{n=1}^\infty \sum_{0 \leqslant j \leqslant n} a_{j}^{(n)}(t) \cos(jx) \end{align}

(5.28)\begin{align} &= \sum_{j=0}^\infty \left(\sum_{n=j}^\infty a_{j}^{(n)}(t)\right) \cos(jx) \end{align}

(5.29)\begin{align} & = \sum_{j=0}^\infty \tilde{a}_j(t) \cos(jx), \end{align}

where

(5.30)\begin{equation} \tilde{a}_j(t) = \sum_{n=j}^\infty a_{j}^{(n)}(t) = \sum_{n=j}^\infty \tilde{c}_{j}^{(n)}(t) a_0^n . \end{equation}

We first perform the cosine Fourier transform for $\cos (jx)$ ($j=1,3$) to the numerical data for the shape of the material interface $\eta (t,x)$ at time $t$ to obtain $\tilde {a}_j (t)$. Note that due to symmetry, the even modes $\cos (jx)$ ($j=0,2,4,\ldots$) do not contribute to the overall amplitude and, consequently, not to the overall velocity. By conducting several numerical simulations for different values of $a_0$ and applying the resulting data to (5.30), we obtain a set of linear equations for the coefficients $\tilde {c}_{j}^{(n)}(t)$ (see (F5) in Appendix F). By solving this set of linear equations, we obtain the corresponding $\tilde {c}_{j}^{(n)}(t)$ from the numerical simulations. Then from the relation $a_{j}^{(n)}(t) = \tilde {c}_{j}^{(n)} (t) a_0^n$ (see (5.30)), we can extract the information about $a_{j}^{(n)}(t)$ from the data of numerical simulations. The details of the theoretical derivation for this method can be found in Appendix F. Finally, by applying the central difference approximation to the value of $a_{j}^{(n)}(t)$, we obtain the results of $\dot {a}_{j}^{(n)}(t)$ from the data of numerical simulations and present the results in figures 6 and 7.

Figure 6. Comparisons between $\dot {a}_j^{(n)}/a_0$ predicted by the nonlinear perturbation solution (solid curves) and that extracted from the numerical simulation (‘$\circ$’ symbols) for the material interface shown in figure 3(c).

Figure 7. Comparisons between $\dot {a}_j^{(n)}/a_0$ predicted by the nonlinear perturbation solution (solid curves) and that extracted from the numerical simulation (‘$\circ$’ symbols) for the material interface shown in figure 3(d).

Figure 8. Normalized dimensionless velocities of bubbles $v_{bb}/{a_0}$ and those of spikes $v_{sp}/{a_0}$ when the material interface is in a linear neutral state due to the stabilizing effects of the horizontal electric field. The horizontal magenta line is the linear prediction for both spikes and bubbles, the black curve is the prediction up to the second order for both spikes and bubbles, the blue curve is the prediction up to the third order for the bubble and the red curve is the prediction up to the third order for the spike. The blue ‘$\circ$’ and red ‘$\times$’ symbols are the numerical data for bubbles and spikes, respectively.

In figures 6 and 7 we investigate the material interfaces in figures 3(c) and 3(d). We compare the theoretical predictions for the normalized interface coefficients $\dot {a}_1^{(1)}/a_0$, $\dot {a}_1^{(3)}/a_0$ and $\dot {a}_3^{(3)}/a_0$ given by (3.54), (3.58) and (3.59) with the extracted results from numerical simulations obtained by the method we have explained above (see also (F7)–(F9)). Note that the coefficients $\dot {a}_1^{(2)}$ and $\dot {a}_2^{(3)}$ are not plotted since they are identically zero. Here $\dot {a}_2^{(2)}$ and the fourth-order coefficients are not plotted since they do not contribute to the overall velocity. Figure 6 corresponds to the system shown in figure 3(c), in which the nonlinear perturbation solution gives an underestimation at late times, and figure 7 corresponds to the system shown in figure 3(d), in which the nonlinear perturbation solution gives an overestimation at late times. In figures 6 and 7 the blue solid curves are the theoretical predictions and the blue ‘$\circ$’ symbols are the numerical results.

Figures 6 and 7 show that $\dot {a}_1^{(1)}/a_0$, $\dot {a}_1^{(3)}/a_0$ and $\dot {a}_3^{(3)}/a_0$ extracted from the data of numerical simulations are in good agreement with our theoretical results for the interface coefficients. This is true even at the late times of the simulations conducted in figures 3(c) and 3(d), at which the theoretical predictions for the overall velocity have already deviated from the numerical results. This is due to the fact that by the method that we discussed above, we are able to extract just the information about each individual coefficient $a_{j}^{(n)}(t)$ from the data of numerical simulations. One may notice that, at late times, the numerical results shown in figure 7(c) are slightly lower than the theoretical values. This is because in the numerical method one needs to introduce a numerical parameter $\delta$ to regularize the singularity on the vortex sheet (see (5.21)–(5.24) and also (12) and (13) in Sohn Reference Sohn2004). This desingularization parameter will slightly smoothen the interface. Such a reduction in the overall velocity is only noticeable for high-frequency modes and for large times. This is why it only appears at the late times of figure 7(c).

Figures 6 and 7 show good agreement between the nonlinear theoretical predictions and the numerical results for $\dot {a}_1^{(1)}$, $\dot {a}_1^{(3)}$ and $\dot {a}_3^{(3)}$ even in the time domain in which the nonlinear perturbation solution underestimates/overestimates (see figure 3c/3d). This confirms that the derived perturbation expansion coefficients given by (3.54), (3.58) and (3.59) are correct. It also confirms that the deviation between the nonlinear perturbation solution and the numerical data is due to the contributions from the higher-order terms since our solution for the overall velocity given by (5.26) is only up to the fourth order of $a_0$. Now we explain why our theoretical prediction given by (5.26) underestimates in figure 3(c) but overestimates in figure 3(d). In general, whether the nonlinear perturbation solution given by (5.26) provides an underestimation or overestimation for the overall velocity of fingers involves complicated interplay between $\dot {a}_1^{(1)}$, $\dot {a}_1^{(3)}$, $\dot {a}_3^{(3)}$ and the contributions from higher orders. Figure 6 shows that $\dot {a}_1^{(1)}$ is positive, but both $\dot {a}_1^{(3)}$ and $\dot {a}_3^{(3)}$ are negative. This is why the corresponding prediction from the nonlinear perturbation solution shown in figure 3(c) underestimates the overall velocity at late times. On the contrary, figure 7 shows that $\dot {a}_1^{(1)}$, $\dot {a}_1^{(3)}$ and $\dot {a}_3^{(3)}$ are all positive. This leads to an overestimate of the overall velocity by the nonlinear perturbation solution, as shown in the corresponding figure 3(d).

In figures 8 and 9 we present the validation study for the velocity of bubbles and that of spikes. Figure 8 is for the comparison of the normalized finger velocities, namely $v_{bb}/{a_0}$ and $v_{sp}/{a_0}$, between the numerical results (bubble:‘$\circ$’, spike: ‘$\times$’) and the theoretical predictions up to the first (magenta curve), second (black curve) and third (bubble: blue curve, spike: red curve) orders given by (4.10)–(4.12). The dimensionless physical parameters for figure 8 are $A = 1$, $A_\epsilon = 0.5$, $E_h^2 = 4/5$ and $a_0 = 0.01$. For this set of parameters, the system is linear neutral, i.e. $\sigma _1 = 0$, which is the turning point between the stable case and the unstable case in the linear analysis. In this case $F_1^{(1)}=0$ and (3.27) gives ${{\rm d}^2 a_{1}^{(1)}}/{{\rm d}t^2} = 0$. From the initial conditions $a_{1}^{(1)}(0) = a_0$ and $\dot {a}_{1}^{(1)}(0)=0$, one obtains $a_{1}^{(1)}(t) = a_0$ and $v^{(1)}(t) = \dot {a}_{1}^{(1)}(t) = 0$. Thus, the linear analysis predicts that the fingers do not grow at all. If we consider the solutions only up to the second order, the spike and the bubble exhibit oscillating motion, but grow at the same rate. Only when we expand the solutions up to the third order in terms of the initial perturbation amplitude, i.e. $a_0^3$, the theoretical prediction correctly shows that the spike grows faster than the bubble. Hence, to provide a qualitatively correct prediction for these fingers, one needs the solutions expanded up to the third order, which are displayed in (4.8) and (4.9).

Figure 9. Normalized dimensionless velocities of bubbles $v_{bb}/{a_0}$ and those of spikes $v_{sp}/{a_0}$ when a horizontal electric field of various dimensionless strengths ${E}_h^2$ is applied. Plots (a,c,e) are for bubbles; plots (b,d,f) are for spikes. Plots (a,b) are for ${E}_h^2=1/6$; plots (c,d) are for ${E}_h^2=1/3$; plots (e,f) are for ${E}_h^2=3/5$. The other physical parameters are listed in the main text. The material interface is unstable in all these cases. The blue solid curves are the predictions from the nonlinear perturbation solution given by (4.8) and (4.9). The black solid curves are the predictions from the linear theory. The ‘$\circ$’ symbols are the data from numerical simulations.

In figure 9 we plot the normalized velocities of bubbles $v_{bb}/{a_0}$ and those of spikes $v_{sp}/{a_0}$ for various strengths $E_h^2 = 1/6, 1/3, 3/5$. The other dimensionless physical parameters are $A = 0.6$, $A_\epsilon = 0.5$, $a_0=0.1$. For all cases in this figure, the system is unstable. The blue solid curves are the predictions from the nonlinear perturbation solution given by (4.8) and (4.9), and the black solid curves are the predictions from the linear theory. The ‘$\circ$’ symbols are the data from numerical simulations. The figures in the left (right) column are for bubble (spike) velocities. Figure 9 shows that the nonlinear theoretical predictions are in good agreement with the numerical results and capture the nonlinear behaviour of bubbles and that of spikes. Figure 9 also shows that as the horizontal electric field becomes stronger, both bubbles and spikes grow more slowly due to the stronger stabilization effect.