3.1. Subcritical tangential electric field, $E_b\lt E_b^{c}$

Following Kochurin et al. (Reference Kochurin, Zubarev and Zubareva2013) and Gao et al. (Reference Gao, Wang and Papageorgiou2024), we introduce two auxiliary functions

(3.19) \begin{align} p = \frac {\varphi +2\lambda \eta }{2},\quad q = \frac {\varphi -2\lambda \eta }{2}, \end{align}

where

(3.20) \begin{align} \lambda = \sqrt {\frac {1}{4}-\frac {E_b(1-\varepsilon )^2}{2(1+\varepsilon )}}. \end{align}

The choice of $\lambda$ follows from (2.17) with $R=1$ that gives the growth rate to be $\omega =\pm \lambda k$ in the subcritical case. Note that $p$ and $q$ correspond to the exponentially growing and decaying eigenmodes of the linearised system. To simplify the model, we assume that the initial condition only contains a linearly growing mode. Therefore, it is reasonable to assume $q = o(\delta ^2)$ and thus omit all nonlinear terms involving $q$ . Equations (3.13) then become

(3.21) \begin{align} p_t - q_t =& -\lambda \mathrm H[p_x] - \lambda \mathrm H[q_x] - \dfrac {U}{4\lambda } \big ( \big (p\mathrm H [p_x]\big )_x - \mathrm H[p p_x]_x \big ), \\[-12pt]\nonumber \end{align}

(3.22) \begin{align} p_t + q_t =& -\lambda \mathrm H[p_x] + \lambda \mathrm H[q_x] + \dfrac {E_b(1-\varepsilon )^3}{8\lambda ^2(1+\varepsilon )^2} \big ( 2\mathrm H\big [p\mathrm H[p_x]\big ]_x + \mathrm H[p_x]^2 - p_x^2 + 2(p p_x)_x \big )\nonumber \\ &- \dfrac {U}{4\lambda }\big ( \mathrm H[p p_x]_x + 2 p_x \mathrm H [p_x] - \big (p \mathrm H [p_x]\big )_x \big ). \end{align}

Adding (3.21) to (3.22) eliminates $q_t$ to give

(3.23) \begin{align} p_t + \lambda \mathrm H[p_x] = - \frac {U}{4\lambda } p_x\mathrm H[p_x] + \frac {E_b(1-\varepsilon )^3}{16\lambda ^2(1+\varepsilon )^2} \big ( 2\mathrm H\big [p\mathrm H[p_x]\big ]_x + \mathrm H[p_x]^2 + p_x^2 + 2 p p_{xx} \big ). \end{align}

An alternative equation for $\eta$ is

(3.24) \begin{align} \eta _t + \lambda \mathrm H[\eta _x] = - \frac {U}{2} \eta _x\mathrm H[\eta _x] + \frac {E_b(1-\varepsilon )^3}{8\lambda (1+\varepsilon )^2} \big ( 2\mathrm H\big [\eta \mathrm H[\eta _x]\big ]_x + \mathrm H[\eta _x]^2 + \eta _x^2 + 2 \eta \eta _{xx} \big ), \end{align}

where we have used the relation $p = 2\lambda \eta + o(\delta ^2)$ that follows from (3.19). Note that prescribing the initial condition

(3.25) \begin{align} \eta (x,0) = \delta \sin (x), \end{align}

to (3.24) is equivalent to giving initial condition

(3.26) \begin{align} \eta (x,0) = \delta \sin (x), \quad \varphi (x,0) = 2\delta \lambda \sin (x), \end{align}

to (3.13).

When $E_b = 0$ , we obtain a quadratic nonlinear model for the classical vortex-sheet problem

(3.27) \begin{align} \eta _t + \frac {1}{2}\mathrm H[\eta _x] + \frac {U}{2} \eta _x\mathrm H[\eta _x] = 0. \end{align}

Equation (3.27) differs from the one derived and studied by Caflisch & Semmes (Reference Caflisch and Semmes1990). To the best of our knowledge, (3.27) is a new vortex-sheet evolution model that has not been derived or studied previously in the literature. Note that (3.27) has similar features to the Constantin–Lax–Majda equation (Constantin et al. Reference Constantin, Lax and Majda1985), which is known to blow up in finite time. For (3.27), one can also derive the exact singularity formation time $t_s$ by noting that it is equivalent to the following complex Burgers’ equation (set $U = 1$ by rescaling):

(3.28) \begin{align} W_t + \frac {\mathrm i}{2}(1+W)W_x = 0, \end{align}

where $W=\eta _x - \mathrm i\mathrm H[\eta _x]$ . To trace the singularity, we solve (3.28) in the complex $x$ -plane by using the method of characteristics. For the initial condition $W(x,0) = \delta \mathrm e^{ ix}$ (consistent with (3.25)), the solution is written implicitly as

(3.29) \begin{align} W = \delta \mathrm e^{i \xi }, \quad \xi = x - \frac {\mathrm i}{2}\big(1+\delta \mathrm e^{ i \xi }\big)t, \end{align}

where $\xi \in \mathbb C$ labels different characteristics. The singularity first forms on $\xi = \xi ^*$ where $x_\xi = 0$ , leading to

(3.30) \begin{align} \xi ^* = 2\pi n + \mathrm i\ln \frac {\delta t}{2}, \end{align}

where $n\in \mathbb{Z}$ . The moment when the singularity intersects the real axis of the complex $x$ -plane corresponds to the physical singularity time $t_s$ , which gives

(3.31) \begin{align} 1+\frac {t_s}{2} + \ln t_s = \ln \frac {2}{\delta }. \end{align}

Replacing $\delta$ by $\delta /2$ (see the explanation below) gives the same result that Moore derived using the Birkhoff–Rott equation (Moore Reference Moore1979). It is worth mentioning that a local asymptotic form of the curvature singularity can be obtained by using this complex-variable approach. Since we aim to find the local behaviour close to the singularity $\xi ^* = \mathrm i\ln (\delta t/2)$ , we write

(3.32) \begin{align} \xi = \tilde {\xi } + \mathrm i\ln \frac {\delta t}{2}, \quad \text{for } \quad |\tilde \xi |\ll 1. \end{align}

Accordingly,

(3.33) \begin{align} \tilde \xi + \mathrm i\ln \frac {\delta t}{2} = x - \frac {\mathrm i}{2}\bigg (1 + \frac {2\mathrm e^{ i \tilde \xi }}{t}\bigg )t. \end{align}

Taylor expanding the right-hand side gives

(3.34) \begin{align} \tilde \xi + \mathrm i\ln \frac {\delta t}{2} = x - \frac {\mathrm i t}{2}- \mathrm i \bigg (1 + \mathrm i \tilde \xi - \frac {\tilde \xi ^2}{2}\bigg ) + O(\tilde \xi ^3). \end{align}

Since we are interested in solutions close to $t_s$ , we write $t = t_s - \tau$ and assume $0\lt \tau \ll 1$ . Then $\tilde \xi$ can be expressed as

(3.35) \begin{align} \tilde \xi \sim \sqrt {2\mathrm i x - \tau (1 + 2/t_s)}, \end{align}

where we have used the identity (3.31). Eventually, we have

(3.36) \begin{align} W = \delta \mathrm e^{ i\xi } = \frac {2}{t_s-\tau }\mathrm e^{ i\tilde \xi } \sim \frac {2}{t_s-\tau }\mathrm e^{i \sqrt {2i x - \tau (1 + 2/t_s)}},\quad \text{for } \quad |x|\ll 1. \end{align}

We are interested in solutions on the real axis ( $\text{Im}\{x\} = 0$ ) at the singularity time ( $\tau = 0$ )

(3.37) \begin{align} W&\sim \frac {2}{t_s} \mathrm e^{ i\sqrt {2 i x}} = \frac {2}{t_s} \mathrm e^{(-1+\text{sgn}(x) i)\sqrt {|x|}}\sim \frac {2}{t_s} \big (1-(1-\text{sgn}(x)\mathrm i)\sqrt {|x|}\big ), \\[-12pt]\nonumber \end{align}

(3.38) \begin{align} W_x&\sim \frac {1}{t_s \sqrt {|x|}}\big (-\text{sgn}(x) + \mathrm i\big ), \end{align}

Taking the real part leads to a local form of $\eta _{xx}$

(3.39) \begin{align} \eta _{xx} \sim \frac {-\text{sgn}(x)}{t_s \sqrt {|x|}}, \quad \text{for } \quad |x|\ll 1. \end{align}

This is almost the same result as that derived by Cowley et al. (Reference Cowley, Baker and Tanveer1999), except our $t_s$ is replaced by their $t_\delta$ , which is solved from

(3.40) \begin{align} \frac {t_\delta }{2} + \ln t_\delta = \ln \frac {1}{\delta }. \end{align}

Note that $t_\delta \sim t_s + o(t_s)$ for $\delta \ll 1$ . More detailed investigations on the properties of (3.27) is left for future studies, but we note here the simplicity of obtaining Moore’s result using the Hamiltonian formulation and a DNO expansion.

To investigate the effect of the electric field on the singularity time $t_s$ as $\delta \rightarrow 0$ , we applied Moore’s analysis to (3.24) and derived the following implicit expression for $t_s$ (see Appendix C for details):

(3.41) \begin{align} 1+\lambda t_s + \ln t_s + \ln \delta + \frac {1}{2}\ln \Bigg (\frac {1}{4}+\bigg (\frac {E_b(1-\varepsilon )^3}{4\lambda (1+\varepsilon )^2}\bigg )^2\Bigg ) = 0. \end{align}

If, however, we prescribe a different initial condition to (3.13), i.e.

(3.42) \begin{align} \eta (x,0) = \delta \sin (x), \quad \varphi (x,0) = 0, \end{align}

then $t_s$ should be governed by the same equation with $\delta$ replaced by $\delta /2$

(3.43) \begin{align} 1+\lambda t_s + \ln t_s + \ln \frac {\delta }{2} + \frac {1}{2}\ln \Bigg (\frac {1}{4}+\bigg (\frac {E_b(1-\varepsilon )^3}{4\lambda (1+\varepsilon )^2}\bigg )^2\Bigg ) = 0, \end{align}

since the initial condition now contains both growing and decaying eigenmodes with initial amplitudes of $\delta /2$ , and $t_s$ is dominated by the growing eigenmode only. In particular, (3.43) recovers Moore’s result when $E_b = 0$

(3.44) \begin{align} 1+\frac {t_s}{2} + \ln t_s = \ln \frac {4}{\delta }. \end{align}

Figure 2.

Results of model (3.24) for $U = 1, \varepsilon = 2, E_b = 1$ and $\delta =0.1$ with initial condition (3.25). Panels (a–d) plot $\eta ,\kappa , V_x$ and $V_{xx}$ at $t = 0, 0.8, 1.6, 2.4, 2.8, 2.92$ .

To examine the validity of (3.41), we solved the model (3.24) numerically by using a pseudospectral method in space and a fourth-order Runge–Kutta scheme in time. In figure 2, we plot solutions for $U = 1$ , $\varepsilon = 2$ , $E_b = 1$ (a subcritical value) and $\delta =0.1$ , with initial condition (3.25). Panels (a–d) show the interface displacement $\eta$ , the curvature $\kappa$ , the electric potential derivative $V_x$ and the second derivative $V_{xx}$ , at times $t = 0, 0.8, 1.6, 2.4, 2.8$ and $2.92$ . For this subcritical TEF, the results show that a curvature singularity appears at $x \approx -0.33$ . The electric field destroys the odd symmetry that solutions would have in its absence given odd-parity initial conditions. Meanwhile, $V_x$ nearly develops a shock at the same position as the curvature singularity, leading $V_{xx}$ to blow up. To numerically estimate $t_s$ , we used the relation

(3.45) \begin{align} \Vert \kappa ^{-1} \Vert _\infty = a(t_s-t) + b(t_s-t)^2,\quad \text{when } \quad t\rightarrow t_s, \end{align}

where $\Vert \boldsymbol{\cdot }\Vert _\infty$ is the $L^\infty$ norm evaluated for $x\in [-\pi ,\pi ]$ , and $t_s$ , $a$ and $b$ are numerically fitted from the computational data by using the least-squares method. A second approach by form fitting the numerical spectrum of $\theta$ to its asymptotic form was applied (Krasny Reference Krasny1986b ; Baker et al. Reference Baker, Caflisch and Siegel1993; Cowley et al. Reference Cowley, Baker and Tanveer1999). Here, we use the following relation:

(3.46) \begin{align} |\hat {\theta }(k,t)| = \frac {A}{k^{\nu _1}}\mathrm e^{-\nu _2 k}, \end{align}

and solve $A,\nu _1$ and $\nu _2$ by the least-squares method. Once we obtain $\nu _2$ at different times, we can extrapolate these numerical data to estimate the moment $t_s$ when $\nu _2$ vanishes. It turns out that the two approaches give nearly the same estimates for $t_s$ for most numerical solutions. We made direct comparisons of our solutions with those of Cowley et al. (Reference Cowley, Baker and Tanveer1999), running our codes with their initial conditions. Our estimates for $t_s$ are in complete agreement with deviations of 2 % or smaller. The largest errors appear for the smallest value of $\delta =0.005$ used, as also pointed out by Cowley et al. (Reference Cowley, Baker and Tanveer1999).

The dependence of $t_s$ on $\delta$ and the electric field strength $E_b$ is presented in figure 3. Figure 3(a) shows the variation of $t_s$ with $\delta$ for $U=1$ and $\varepsilon =2$ for the three subcritical electric field strengths $E_b=0, 1, 1.4$ ; the triangles represent $t_s$ obtained from the numerical solution of (3.24) with initial condition (3.25), and the solid curves are the asymptotic estimates (3.41). The agreement between the numerical results and the analytical predictions is remarkable, clearly showing that $t_s \propto \ln (1/\delta )$ as $\delta \rightarrow 0$ , in line with Moore’s classical result. On the other hand, when $E_b$ increases towards $E_b^{c}$ , which is $1.5$ for current parameters, the numerical results follow the analytic prediction qualitatively but there is a visible discrepancy. In fact, when $E_b \rightarrow E_b^{c}$ , (3.41) yields that $t_s \rightarrow 0$ for fixed $\delta$ . This can be clearly seen from figures 3(b) and 3(c) which plot $t_s$ versus $E_b$ obtained from (3.41) for $\delta = 0.1$ . It can be seen that $t_s$ monotonically increases with $E_b$ until a sudden dip appears as $E_b \rightarrow E_b^{c}$ . This indicates that the model (3.24) breaks down because the nonlinear terms overwhelm the linear terms when $\lambda \rightarrow 0$ , thus violating the basic assumption of weak nonlinearity. Full computations of the Euler equations in this regime are presented later.

Figure 3.

(a) Value of $t_s$ versus $\delta$ . Triangles and solid curves are numerical results of (3.24) and (3.41) for $U = 1$ and $\varepsilon = 2$ . The three groups of results correspond to $E_b = 0, 1, 1.4$ from the bottom to top. (b–c) Value of $t_s$ versus $E_b$ obtained from the (3.41) for $\delta = 0.1$ .

Figure 3. Long description

The image contains three line graphs labeled as Panel A, Panel B, and Panel C. Panel A shows a line graph with the vertical axis labeled as t subscript s and the horizontal axis labeled as delta. The graph includes multiple lines with triangles marking specific data points. An arrow indicates that E subscript b increases. Panel B presents a line graph with the vertical axis labeled as t subscript s and the horizontal axis labeled as E subscript b. Different curves are labeled with epsilon values of 0.2, 0.4, 0.6, and 0.8. Panel C also shows a line graph with the vertical axis labeled as t subscript s and the horizontal axis labeled as E subscript b. Curves are labeled with epsilon values of 2, 4, 6, and 8. Each panel illustrates how t subscript s varies with different parameters, showing trends and relationships in the data.

3.2. Supercritical tangential electric field, $E_b\gt E_b^{c}$

We first recast the system (3.13) into a single equation for $\eta$

(3.47) \begin{align} \eta _{tt} - c^2\eta _{xx} -\frac {U}{2}\Big ( 2\mathrm H[\eta ]\eta _{xt} + \mathrm H\big [\mathrm H[\eta ]^2\big ]_{xt} \Big )_x + \frac {E_b(1-\varepsilon )^3}{4(1+\varepsilon )^2}\Big ( 2\mathrm H[\eta ]\eta _{xx} + \mathrm H\big [\mathrm H[\eta ]^2\big ]_{xx} \Big )_x = 0, \end{align}

where

(3.48) \begin{align} c^2 = \frac {E_b(1-\varepsilon )^2}{2(1+\varepsilon )} - \frac {1}{4}\gt 0. \end{align}

In general, solutions of (3.47) contain both right- and left-travelling modes. By assuming that only one uni-directional mode dominates, the same strategy used to derive the Korteweg–de Vries equation from the Boussinesq equation, we can simplify (3.47) further. We introduce $\tau = \delta t$ and $\xi = x - ct$ , and expand $\eta$ asymptotically as

(3.49) \begin{align} \eta \sim \eta _0(\xi ,\tau ) + \delta \eta _1(\xi ,\tau ) + {\cdots} . \end{align}

The leading order of (3.47) reads

(3.50) \begin{align} -2\delta c\eta _{0,\tau \xi } + 2\delta c\gamma \,\Big ( 2\mathrm H[\eta _0]\eta _{0,\xi \xi } + \mathrm H\big (\mathrm H[\eta _0]^2\big )_{\xi \xi } \Big )_\xi = O(\delta ^3), \end{align}

where

(3.51) \begin{align} \gamma = \frac {E_b(1-\varepsilon )^3}{8c(1+\varepsilon )^2} + \frac {U}{4}. \end{align}

Keeping leading-order terms and integrating (3.50) once with respect to $\xi$ , we have a simplified equation for $\eta _0$ , namely,

(3.52) \begin{align} \eta _{0,\tau } - \gamma \Big ( 2\mathrm H[\eta _0]\eta _{0,\xi \xi } + \mathrm H\big [\mathrm H[\eta _0]^2\big ]_{\xi \xi } \Big ) = 0. \end{align}

In arriving at (3.52) we have used the property $\mathrm H[fg] = \mathrm H[f]g+ f\mathrm H[g] + \mathrm H [\mathrm H[f]\mathrm H[g] ]$ . The constant of integration is zero because we require conservation of mass

(3.53) \begin{align} \int _0^{2\pi } \eta (\xi ,\tau )\,\mathrm d\xi = 0. \end{align}

We numerically solve (3.52) by using a pseudospectral method in space and a fourth-order Runge–Kutta scheme in time, as before. In the results of figure 4 we set $\gamma = 1$ and prescribe the initial condition $\eta _0(\xi ,0) = 0.1\sin (\xi )$ . A parity analysis of (3.52) shows that $\eta$ remains an odd function of $\xi$ . Successive solutions from $t = 0$ to $t = 2.276$ are plotted in figure 4(a), showing that the interface gradually steepens on the left side of the crest. The curvature at $t=2.276$ is plotted in figure 4(b), indicating the appearance of a curvature singularity. Close to the singularity time, we can seek a self-similar solution

(3.54) \begin{align} \eta _0(x,t) \sim (t_s - t)^a g (\zeta ), \quad \zeta = \frac {\xi }{(t_s-t)^b} = O(1), \end{align}

where $a$ and $b$ are positive constants. Substituting (3.54) into (3.52), we have the balance

(3.55) \begin{align} a - 2b = -1, \end{align}

and

(3.56) \begin{align} \Vert \kappa \Vert _\infty \sim \Vert \eta _{0,xx}\Vert _\infty \sim (t_s - t)^{a-2b} = (t_s - t)^{-1}. \end{align}

Considering this relation is true only when $t_s-t\ll 1$ , we can still use the assumed relation (3.45) to estimate $t_s$ . In figure 4(c), we plot the inverse curvature magnitude $\Vert \kappa \Vert _\infty ^{-1}$ and the least-squares fitted quadratic curve. The estimated singularity time is $t_s = 2.35$ .

Figure 4.

Numerical solution of (3.52) with $\gamma = 1$ and initial condition $\eta _0(\xi ,0) = 0.1\sin (\xi )$ . (a) Interfaces at $t = 0, 0.5, 1, 1.5, 2$ and $2.276$ from bottom to top. (b) Curvature profile at $t = 2.276$ . (c) Inverse curvature amplitude versus time. The blue curve is the numerical result, and the red dotted line is the fitted curve using the quadratic relation (3.45). The red cross sign labels the estimated singularity time at $t = 2.35$ .

Note that (3.52) is invariant under the transformation

(3.57) \begin{align} \xi \rightarrow -\xi ,\quad \gamma \rightarrow -\gamma . \end{align}

Physically, changing the sign of $\xi$ by reversing the direction of the $\xi$ -axis, implies that $c$ and $U$ also change signs, leading to $\gamma \to -\gamma$ . In particular, for initial condition $\eta (\xi ,0) = 0.1\sin (\xi )$ and $\gamma =-1$ , the interface steepens on the right side of the crest. In general, the transformation

(3.58) \begin{align} \eta _0\rightarrow \delta \eta _0, \quad t\rightarrow \frac {t}{\delta |\gamma |}, \end{align}

rescales the initial amplitude to $1$ , and $\gamma$ to $1$ if it is positive, or $-1$ otherwise. Therefore, the singularity time is inversely proportional to $\delta |\gamma |$ . For initial condition $\eta _0(\xi ,0) = \delta \sin (\xi )$ , the result in figure 4 gives an estimate of the singularity time to be

(3.59) \begin{align} t_s = \frac {0.235}{\delta }\bigg |\frac {E_b(1-\varepsilon )^3}{8c(1+\varepsilon )^2} + \frac {U}{4}\bigg |^{-1}. \end{align}

Note that (3.52) is derived under the assumption that only one uni-directional component dominates. For model (3.13), this is guaranteed by giving the following initial conditions:

(3.60) \begin{align} \eta (x,0) = \delta \sin (x),\quad \varphi (x,0) = -2\delta c \cos (x). \end{align}

If, however, the initial conditions are

(3.61) \begin{align} \eta (x,0) = \delta \sin (x),\quad \varphi (x,0) = 0, \end{align}

there coexist two counter-travelling modes. Depending on the values of $U$ and $\varepsilon$ , the singularity first develops in the component travelling at a speed of $\text{sgn} ((1-\varepsilon )U )|c|$ . With $U\gt 0$ , which is the case we focus on, a curvature singularity forms and travels to the left (right) for $\varepsilon \gt 1$ ( $\varepsilon \lt 1$ ). Thus, an estimate of the singularity time is

(3.62) \begin{align} t_s = \frac {0.47}{\delta }\bigg ( \bigg |\frac {E_b(1-\varepsilon )^3}{8c(1+\varepsilon )^2}\bigg | + \frac {|U|}{4} \bigg )^{-1}. \end{align}

Note that Moore’s analysis is not adequate for model (3.52) because this method requires $t_s = o(1/\delta )$ (see Appendix C). For model (3.52), however, it is clear that $t_s\propto 1/\delta$ , which makes the asymptotic expansion non-uniformly valid when $t\rightarrow t_s$ .

Figure 5.

Typical plots of $|\gamma |^{-1}$ versus $E_b$ for $cU\gt 0$ : (a) $\varepsilon \lt 1$ ; (b) $\varepsilon \gt 1$ .

Figure 5. Long description

Two line graphs depict the relationship between inverse gamma and Eb for different values of epsilon. Panel A: The line graph shows the inverse of gamma on the vertical axis and Eb on the horizontal axis. The graph includes multiple lines representing different values of epsilon, specifically 0.1, 0.3, 0.5, and 0.7. Each line shows a trend where the inverse of gamma initially increases with Eb, reaches a peak, and then decreases. Panel B: The line graph shows the inverse of gamma on a logarithmic scale on the vertical axis and Eb on the horizontal axis. The graph includes multiple lines representing different values of epsilon, specifically 2, 2.5, and 3. Each line shows a trend where the inverse of gamma initially increases with Eb, reaches a peak, and then decreases, with the peaks becoming more pronounced as epsilon increases.

Depending on the values of $c, U$ , and $\varepsilon$ , $|\gamma |^{-1}$ either remains bounded or blows up at some specfici values of $E_b$ . In general, we have the following results that arise directly from (3.51):

(3.63) \begin{align} \max \{|\gamma |^{-1}\} = \begin{cases} \infty ,& \text{at}\,\, E_b = \dfrac { (1+\varepsilon )^2(1\pm \sqrt {\varepsilon })^2}{(1-\varepsilon )^4}, \,\,\text{when}\,\, c U\gt 0, \varepsilon \gt 1,\\[12pt] \infty ,& \text{at}\,\, E_b = \dfrac { (1+\varepsilon )^2(1\pm \sqrt {\varepsilon })^2}{(1-\varepsilon )^4}, \,\,\text{when}\,\, c U\lt 0, \varepsilon \lt 1,\\[12pt] 2(1+\varepsilon ),& \text{at}\,\,E_b = \dfrac { (1+\varepsilon )}{(1-\varepsilon )^2}, \,\,\text{when}\,\, c U\gt 0, \varepsilon \lt 1,\\[12pt] \dfrac {2(1+\varepsilon )}{\varepsilon },& \text{at}\,\,E_b = \dfrac { (1+\varepsilon )}{(1-\varepsilon )^2}, \,\,\text{when}\,\, c U\lt 0, \varepsilon \gt 1. \end{cases} \end{align}

It is worth mentioning that the two cases $c U\gt 0$ and $c U\lt 0$ can be transformed to each other by using the transformations

(3.64) \begin{align} U \rightarrow - U,\quad \varepsilon \rightarrow \frac {1}{\varepsilon },\quad E_b \rightarrow \varepsilon E_b, \end{align}

which can be interpreted as turning the two-layer system upside down. Since the definition of $E_b$ relies on the permittivity constant of the lower fluid, a factor $\varepsilon$ emerges when we flip the system. Equation (3.63) indicates that the curvature singularity may be significantly delayed when $E_b$ is close to the values that lead $ |\gamma |^{-1}$ to blow up. Strictly speaking, however, when $|\gamma | = O(\delta )$ , the cubic nonlinearity becomes equally important as the quadratic one, meaning that (3.52) is no longer an asymptotically valid model. Although the precise behaviour of solutions in this regime is unclear, i.e. whether the curvature singularity can be suppressed by a possible saturation arising from the cubic nonlinearity, it is certain that $t_s$ becomes significantly large, typically $O(1/\delta ^2)$ , and we shall confirm this assertion later in solutions of the Euler equations. Figure 5 plots some typical $|\gamma |^{-1}$ curves for $c U\gt 0$ . It is interesting to note that these curves indicate that $t_s$ decreases with $E_b$ when it becomes sufficiently large. When $E_b \rightarrow \infty$

(3.65) \begin{align} t_s \propto \frac {1}{\sqrt {E_b}}. \end{align}

At this stage the legitimacy of this estimate is unclear, as is also the prediction that $t_s\rightarrow 0$ when $E_b\rightarrow E_b^{c}$ . For similar reasons to those mentioned in § 5.2, i.e. when $|\gamma | = O(1/\delta )$ , the leading quadratic nonlinearity becomes as important as the linear terms, which may violate the basic weakly nonlinear assumption used to derive the model (3.13). Nevertheless, we shall check the validity of these asymptotic estimates in comparisons with full Euler simulations.

4.1. The $\theta -s_\alpha$ Lagrangian formulation

We parameterise the interface by a Lagrangian parameter $\alpha \in [0,1]$ . It is well known that for a two-dimensional material curve, its tangent angle $\theta (\alpha ,t)$ and arclength derivative $s_\alpha (\alpha ,t)$ are governed by the following $\theta -s_\alpha$ formulation (see Yang Reference Yang1992; Hou et al. Reference Hou, Lowengrub and Shelley1994, Reference Hou, Lowengrub and Shelley1997; Guan & Vanden-Broeck Reference Guan and Vanden-Broeck2025):

(4.1) \begin{align} \theta _{t} = \frac {N_{\alpha }}{s_{\alpha }} + \frac {\theta _{\alpha } T}{s_{\alpha }}, \\[-12pt]\nonumber \end{align}

(4.2) \begin{align} s_{\alpha t} = T_{\alpha } - \theta _{\alpha } N, \end{align}

where $T$ and $N$ represent the tangential and normal velocities that are defined below. Using the relation

(4.3) \begin{align} \frac {z_\alpha }{s_\alpha } = \mathrm e^{ i\theta },\quad z(\alpha ,t)=x(\alpha ,t)+\mathrm{i}\eta (\alpha ,t), \end{align}

one can construct the interface by the following integration:

(4.4) \begin{align} z(\alpha ,t) = \int _0^\alpha \mathrm e^{ i\theta (\alpha ',t)} s_\alpha (\alpha ',t)\,\mathrm d\alpha ' + z_0. \end{align}

The function $z_0= z(0,t) = x_0 + \mathrm i\eta _0$ is the complex coordinate of the leftmost Lagrangian marker, hence

(4.5) \begin{align} z_{0,t} = \big (T(0,t)+\mathrm i N(0,t)\big )\mathrm e^{i\theta (0,t)}. \end{align}

The tangential and normal velocities of the lower and upper fluid particles on the interface are

(4.6) \begin{align} T_{\!j} = -(-1)^j\frac {U}{2}\cos \theta + \frac {\partial \phi _{\!j}}{\partial \boldsymbol{t}},\quad N = (-1)^j\frac {U}{2}\sin \theta + \frac {\partial \phi _{\!j}}{\partial \boldsymbol{n}},\quad \text{for} \,\, j = 1 \,\,\text{and}\,\,2. \end{align}

Note that the normal velocity is continuous across the interface, so we can drop the subscript. On the other hand, $T_1 \neq T_2$ , so there is freedom to define $T$ in (4.1) and (4.2). Following Baker et al. (Reference Baker, Meiron and Orszag1982) and Grue et al. (Reference Grue, Friis, Palm and Rusaas1997), we write $T$ as a linear combination of $T_1$ and $T_2$

(4.7) \begin{align} T = \beta T_1 + (1-\beta ) T_2, \end{align}

where $\beta \in [0,1]$ . Note that $\beta = 0$ or $1$ corresponds to choosing the upper or the lower fluid particles as the Lagrangian markers, respectively. In our simulations, we find that $\beta = 0.5$ is the best choice for subcritical electric fields (this is also the usual choice in the Birkhoff–Rott equation). For supercritical electric fields, however, there is no general way to choose the optimal value of $\beta$ . Here by optimal, we mean the specific value of $\beta$ such that $s_\alpha$ takes a global minimum value at the singular point among all choices. In this case, we test different values of $\beta$ and take the one giving the smallest energy error (4.52).

Defining $\varphi _{\!j}(\alpha ,t) = \phi _{\!j} (x(\alpha ,t),\eta (\alpha ,t),t )$ , we have

(4.8) \begin{align} \varphi _{j,t} = \frac {D\phi _{\!j}}{Dt} = \phi _{j,t} + T\frac {\partial \phi _{\!j}}{\partial \boldsymbol{t}} + N\frac {\partial \phi _{\!j}}{\partial \boldsymbol{n}},\quad \text{for} \quad j = 1 \,\,\text{and}\,\,2, \end{align}

where $D/Dt$ denotes the material derivative associated with the Lagrangian markers. The dimensionless Bernoulli equation is expressed as

(4.9) \begin{align} \varphi _{t} - T\bigg [\frac {\partial \phi _1}{\partial \boldsymbol{t}}- R\frac {\partial \phi _2}{\partial \boldsymbol{t}}\bigg ] - N\bigg [\frac {\partial \phi _1}{\partial \boldsymbol{n}} - R\frac {\partial \phi _2}{\partial \boldsymbol{n}}\bigg ] + \frac {T_1^2+N^2}{2} - R\frac {T_2^2+N^2}{2} \nonumber \\ + \frac {1-R}{F_r^2}\eta + \frac {E_b}{2}\big [\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\varSigma }\boldsymbol{\cdot }\boldsymbol{n}\big ]_2^1 = 0, \end{align}

where $\varphi (\alpha ,t) = \varphi _1(\alpha ,t)-R\varphi _2(\alpha ,t)$ . Using relations (4.6), we have

(4.10) \begin{align} \varphi _{t} &+ \frac {1-2\beta }{2}\big (T_1^2+\textit{RT}_2^2\big ) - \big (1-\beta -R\beta \big )T_1 T_2 - \frac {1-R}{2}N^2\nonumber \\ &+ \frac {U}{2}(1+R)\big (T\cos \theta -N\sin \theta \big ) + \frac {1-R}{F_r^2}\eta + \frac {E_b}{2}\big [\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\varSigma }\boldsymbol{\cdot }\boldsymbol{n}\big ]_2^1 = 0, \end{align}

where

(4.11) \begin{align} \big [\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\varSigma }\boldsymbol{\cdot }\boldsymbol{n}\big ]_2^1 = \bigg (\frac {\partial V_1}{\partial \boldsymbol{n}}-\sin \theta \bigg )^2 - \bigg (\frac {\partial V_1}{\partial \boldsymbol{t}}+\cos \theta \bigg )^2 \nonumber \\ - \varepsilon \bigg (\frac {\partial V_1}{\partial \boldsymbol{n}}-\sin \theta \bigg )^2 + \varepsilon \bigg (\frac {\partial V_1}{\partial \boldsymbol{t}}+\cos \theta \bigg )^2. \end{align}

Therefore, we have four time-evolution equations (4.1), (4.2), (4.5) and (4.10) for the unknowns $\theta$ , $s_\alpha$ , $z_0$ and $\varphi$ .

4.2. The Cauchy integral formula and derivation of integral equations of the second kind

In each fluid layer we define the complex velocity

(4.12) \begin{align} w_{\!j} = \phi _{j,x} - \mathrm i\phi _{j,y},\quad \text{for} \,\, j = 1 \,\,\text{and}\,\,2. \end{align}

They are analytic functions and satisfy Cauchy’s integral formula

(4.13) \begin{align} w_{\!j}(z_0) = \frac {1}{\pi \mathrm i} \mathrm{PV}\oint _{C_{\!j}} \frac {w_{\!j}(z)}{z-z_0}\,\mathrm dz,\quad \text{for} \,\, j = 1 \,\,\text{and}\,\,2, \end{align}

where PV denotes the Cauchy principal value, $z$ and $z_0$ are points on the interface and $C_1$ and $C_2$ are contours along the interface such that the fluid $1$ and $2$ are on their left side. Alternatively, we use the transformation $\zeta = \mathrm e^{-i z}$ to cast (4.13) into

(4.14) \begin{align} w_{\!j}(\zeta _0) = \frac {1}{\pi \mathrm i} \mathrm{PV}\oint _{\varOmega _{\!j}} \frac {w_{\!j}(\zeta )}{\zeta -\zeta _0}\,\mathrm d\zeta ,\quad \text{for} \,\, j = 1 \,\,\text{and}\,\,2. \end{align}

The new contours $\varOmega _1$ and $\varOmega _2$ are the images of $C_1$ and $C_2$ , and become closed under the transformation. Using (4.3) and the identities

(4.15) \begin{align} w_{\!j} \,\mathrm e^{ i\theta } = \frac {\partial \phi _{\!j}}{\partial \boldsymbol{t}} - \mathrm i \frac {\partial \phi _{\!j}}{\partial \boldsymbol{n}},\quad \text{for} \,\, j = 1 \,\,\text{and}\,\,2, \end{align}

we can rewrite (4.14) as

(4.16) \begin{align} \frac {\partial \phi _{1}}{\partial \boldsymbol{t}}(\alpha ') &- \mathrm i \frac {\partial \phi _{1}}{\partial \boldsymbol{n}}(\alpha ') = \frac {1}{\pi } \int _0^1 \frac {\mathrm e^{ i\theta (\alpha ')}s_{\alpha }(\alpha )}{1 - \mathrm e^{ i(z(\alpha )-z(\alpha '))}}\bigg ( \frac {\partial \phi _{1}}{\partial \boldsymbol{t}}(\alpha ) - \mathrm i \frac {\partial \phi _{1}}{\partial \boldsymbol{n}}(\alpha ) \bigg )\,\mathrm d\alpha , \\[-12pt]\nonumber \end{align}

(4.17) \begin{align} \frac {\partial \phi _{2}}{\partial \boldsymbol{t}}(\alpha ') &- \mathrm i \frac {\partial \phi _{2}}{\partial \boldsymbol{n}}(\alpha ') = -\frac {1}{\pi } \int _0^1 \frac {\mathrm e^{ i\theta (\alpha ')}s_{\alpha }(\alpha )}{1 - \mathrm e^{ i(z(\alpha )-z(\alpha '))}}\bigg ( \frac {\partial \phi _{2}}{\partial \boldsymbol{t}}(\alpha ) - \mathrm i \frac {\partial \phi _{2}}{\partial \boldsymbol{n}}(\alpha ) \bigg )\,\mathrm d\alpha , \end{align}

where we have dropped the PV for brevity. Similarly, the voltage potentials satisfy the following analogous equations:

(4.18) \begin{align} \frac {\partial V_{1}}{\partial \boldsymbol{t}}(\alpha ') &- \mathrm i \frac {\partial V_{1}}{\partial \boldsymbol{n}}(\alpha ') = \frac {1}{\pi } \int _0^1 \frac {\mathrm e^{ i\theta (\alpha ')}s_{\alpha }(\alpha )}{1 - \mathrm e^{ i(z(\alpha )-z(\alpha '))}}\bigg ( \frac {\partial V_{1}}{\partial \boldsymbol{t}}(\alpha ) - \mathrm i \frac {\partial V_{1}}{\partial \boldsymbol{n}}(\alpha ) \bigg )\,\mathrm d\alpha , \\[-12pt]\nonumber \end{align}

(4.19) \begin{align} \frac {\partial V_{2}}{\partial \boldsymbol{t}}(\alpha ') &- \mathrm i \frac {\partial V_{2}}{\partial \boldsymbol{n}}(\alpha ') = -\frac {1}{\pi } \int _0^1 \frac {\mathrm e^{ i\theta (\alpha ')}s_{\alpha }(\alpha )}{1 - \mathrm e^{ i(z(\alpha )-z(\alpha '))}}\bigg ( \frac {\partial V_{2}}{\partial \boldsymbol{t}}(\alpha ) - \mathrm i \frac {\partial V_{2}}{\partial \boldsymbol{n}}(\alpha ) \bigg )\,\mathrm d\alpha . \end{align}

These integral equations must be solved with appropriate boundary conditions. Since $\varphi$ is given as initial condition, we have

(4.20) \begin{align} \frac {\partial \phi _1}{\partial \boldsymbol{t}} - R \frac {\partial \phi _2}{\partial \boldsymbol{t}} = \frac {\varphi _\alpha }{s_\alpha }. \end{align}

The normal velocity $N$ is continuous across the interface, thus giving

(4.21) \begin{align} \frac {\partial \phi _1}{\partial \boldsymbol{n}} - \frac {\partial \phi _2}{\partial \boldsymbol{n}} = U\sin \theta . \end{align}

For the voltage potentials, we define $V(\alpha ,t) = V_1 (x(\alpha ,t),\eta (\alpha ,t),t ) = V_2 (x(\alpha ,t), \eta (\alpha ,t),t )$ . The tangential derivatives must be identical, meaning

(4.22) \begin{align} \frac {\partial V_{1}}{\partial \boldsymbol{t}} = \frac {\partial V_{2}}{\partial \boldsymbol{t}} = \frac {\partial V}{\partial \boldsymbol{t}}. \end{align}

For their normal derivatives, (2.10) gives

(4.23) \begin{align} \frac {\partial V_1}{\partial \boldsymbol{n}} - \varepsilon \frac {\partial V_2}{\partial \boldsymbol{n}} = (1-\varepsilon )\sin \theta . \end{align}

These are appropriate boundary conditions for (4.16)–(4.19). We illustrate how to embed them into the boundary integral equations for (4.16) and (4.17). Separating the real and imaginary parts, we have

(4.24) \begin{align} \frac {\partial \phi _1}{\partial \boldsymbol{n}}(\alpha ') &= \int _0^1 \bigg (\mathcal A(\alpha ',\alpha ) \frac {\partial \phi _1}{\partial \boldsymbol{n}}(\alpha ) - \mathcal B(\alpha ',\alpha )\frac {\partial \phi _1}{\partial \boldsymbol{t}}(\alpha )\bigg )\,\mathrm d\alpha , \\[-12pt]\nonumber \end{align}

(4.25) \begin{align} \frac {\partial \phi _2}{\partial \boldsymbol{n}}(\alpha ') &= \int _0^1 \bigg (\mathcal B(\alpha ',\alpha ) \frac {\partial \phi _2}{\partial \boldsymbol{t}}(\alpha ) - \mathcal A(\alpha ',\alpha ) \frac {\partial \phi _2}{\partial \boldsymbol{n}}(\alpha )\bigg )\,\mathrm d\alpha , \\[-12pt]\nonumber \end{align}

(4.26) \begin{align} \frac {\partial \phi _1}{\partial \boldsymbol{t}}(\alpha ') &= \int _0^1 \bigg (\mathcal A(\alpha ',\alpha ) \frac {\partial \phi _1}{\partial \boldsymbol{t}}(\alpha ) + \mathcal B(\alpha ',\alpha ) \frac {\partial \phi _1}{\partial \boldsymbol{n}}(\alpha )\bigg )\,\mathrm d\alpha , \\[-12pt]\nonumber \end{align}

(4.27) \begin{align} \frac {\partial \phi _2}{\partial \boldsymbol{t}}(\alpha ') &= -\int _0^1 \bigg (\mathcal A(\alpha ',\alpha ) \frac {\partial \phi _2}{\partial \boldsymbol{t}}(\alpha ) + \mathcal B(\alpha ',\alpha ) \frac {\partial \phi _2}{\partial \boldsymbol{n}}(\alpha )\bigg )\,\mathrm d\alpha , \end{align}

where

(4.28) \begin{align} \mathcal A(\alpha ',\alpha ) = \frac 1\pi \text{Re}\bigg (\frac {\mathrm e^{ i\theta (\alpha ')}s_{\alpha }(\alpha )}{1 - \mathrm e^{ i(z(\alpha )-z(\alpha '))}} \bigg ),\quad \mathcal B(\alpha ',\alpha ) = \frac 1\pi \text{Im}\bigg (\frac {\mathrm e^{ i\theta (\alpha ')}s_{\alpha }(\alpha )}{1 - \mathrm e^{ i(z(\alpha )-z(\alpha '))}} \bigg ). \end{align}

We introduce

(4.29) \begin{align} \mathcal T = \frac {\partial \phi _1}{\partial \boldsymbol{t}} - R\frac {\partial \phi _2}{\partial \boldsymbol{t}}, \quad \mathcal U = \frac {\partial \phi _1}{\partial \boldsymbol{t}} + \frac {\partial \phi _2}{\partial \boldsymbol{t}},\quad \mathcal J = \frac {\partial \phi _1}{\partial \boldsymbol{n}} - \frac {\partial \phi _2}{\partial \boldsymbol{n}}, \quad \mathcal Q = \frac {\partial \phi _1}{\partial \boldsymbol{n}} + R\frac {\partial \phi _2}{\partial \boldsymbol{n}}, \end{align}

to express the original unknowns as

(4.30) \begin{align} \frac {\partial \phi _1}{\partial \boldsymbol{t}} = \frac {\mathcal T+R \mathcal U}{1+R}, \quad \frac {\partial \phi _2}{\partial \boldsymbol{t}} = \frac {\mathcal U - \mathcal T}{1+R}, \quad \frac {\partial \phi _1}{\partial \boldsymbol{n}} = \frac {\mathcal Q + R \mathcal J}{1+R}, \quad \frac {\partial \phi _2}{\partial \boldsymbol{n}} = \frac {\mathcal Q - \mathcal J}{1+R}. \end{align}

Multiplying (4.25) by $R$ and adding it to (4.24), we have

(4.31)

Adding (4.26) to (4.27), we obtain

(4.32)

Similarly, after introducing

(4.33) \begin{align} \mathcal S = \frac {\partial V_1}{\partial \boldsymbol{n}} - \varepsilon \frac {\partial V_2}{\partial \boldsymbol{n}}, \quad \mathcal W = \frac {\partial V_1}{\partial \boldsymbol{n}} + \frac {\partial V_2}{\partial \boldsymbol{n}}, \end{align}

(4.18) and (4.19) can be rewritten as

(4.34)

(4.35)

The normal derivatives of the voltage potentials can be calculated from

(4.36) \begin{align} \frac {\partial V_1}{\partial \boldsymbol{n}} = \frac {\mathcal S + \varepsilon \mathcal W}{1+\varepsilon }, \quad \frac {\partial V_2}{\partial \boldsymbol{n}} = \frac {\mathcal W - \mathcal S}{1+\varepsilon }. \end{align}

Therefore, we obtain four decoupled Fredholm integral equations of the second kind. In particular, (4.31) and (4.25) are automatically solved when $R = 1$ . It is worth mentioning that the current boundary integral formulation can be readily applied to flows with flat or uneven solid boundaries. Interested readers are referred to Guan & Vanden-Broeck (Reference Guan and Vanden-Broeck2025) for more details.

4.3. Numerical implementation

We introduce $M$ Lagrangian markers

(4.37) \begin{align} \alpha _m = \frac {m}{M},\quad m = 0,1,{\cdots} ,M-1. \end{align}

All spatial derivatives and integrations are calculated using fast Fourier transforms (FFTs). When constructing the interface from (4.4), $x$ and $\eta$ are calculated via

(4.38) \begin{align} x(\alpha ,t) &= \int _0^\alpha \cos \big (\theta (\alpha ',t)\big ) s_\alpha (\alpha ',t)\,\mathrm d\alpha ' + x_0, \\[-12pt]\nonumber \end{align}

(4.39) \begin{align} \eta (\alpha ,t) &= \int _0^\alpha \sin \big (\theta (\alpha ',t)\big ) s_\alpha (\alpha ',t)\,\mathrm d\alpha ' + \eta _0. \end{align}

Note that directly evaluating (4.38) may violate the periodicity condition

(4.40) \begin{align} x(1,t) - x(0,t) = 2\pi , \end{align}

because this requires the mean of $\cos (\theta (\alpha ,t) ) s_\alpha (\alpha ,t)$ to be exactly $2\pi$ . Following Ambrose et al. (Reference Ambrose, Camassa, Marzuola, McLaughlin, Robinson and Wilkening2022), we rewrite (4.38) as

(4.41) \begin{align} x(\alpha ,t) = \int _0^\alpha \mathcal P\big (\cos \big (\theta (\alpha ',t)\big ) s_\alpha (\alpha ',t)\big )\,\mathrm d\alpha ' + 2\pi \alpha + x_0, \end{align}

where $\mathcal P$ is the zero-mean operator

(4.42) \begin{align} \mathcal P\big (f(\alpha )\big ) = f(\alpha ) - \int _0^1 f(\alpha )\,\mathrm d\alpha . \end{align}

Once the interface is computed, the boundary integral (4.31)–(4.35) can be solved. Note that when $\alpha \rightarrow \alpha _m$

(4.43) \begin{align} \frac {\mathrm e^{ i\theta (\alpha _m)}s_{\alpha }(\alpha )}{1 - \mathrm e^{ i(z(\alpha )-z(\alpha _m))}} = \frac {\mathrm i}{\alpha -\alpha _m} + \frac {\mathrm i s_{\alpha \alpha }(\alpha _m)}{2 s_\alpha (\alpha _m)} + \frac {s_\alpha (\alpha _m)\mathrm e^{ i\theta (\alpha _m)}}{2} + \frac {\theta _\alpha (\alpha _m)}{2} + O(\alpha -\alpha _m). \end{align}

Therefore, $\mathcal B(\alpha _m,\alpha )$ is singular at $\alpha = \alpha _m$ . For a smooth periodic function $f(\alpha )$

(4.44) \begin{align} \int _0^1 \mathcal B(\alpha _m,\alpha ) f(\alpha )\,\mathrm d\alpha = \int _{\alpha _m-\frac 12}^{\alpha _m+\frac 12} \bigg [\mathcal B(\alpha _m,\alpha )f(\alpha ) - \frac {f(\alpha _m)}{\pi (\alpha -\alpha _m)}\bigg ]\,\mathrm d\alpha \nonumber \\ + \underbrace {\frac {f(\alpha _m)}{\pi } \mathrm{PV}\int _{\alpha _m-\frac 12}^{\alpha _m+\frac 12} \frac {1}{\alpha -\alpha _m} \,\mathrm d\alpha }_{0}, \end{align}

where the integrals are understood to be principal valued, and we have used periodicity to change the lower and upper limits. The first integrand on the right-hand side of (4.44) becomes regular and tends to $f_\alpha /\pi$ when $\alpha \rightarrow \alpha _m$ . The trapezoid rule below provides spectral accuracy

(4.45) \begin{align} \int _0^1 \mathcal B(\alpha _m,\alpha ) f(\alpha )\,\mathrm d\alpha \approx \frac {1}{M}\sum _{n = 0,n \neq m}^{M-1} \mathcal B(\alpha _m,\alpha _n) f(\alpha _n) + \frac {f_\alpha (\alpha _m)}{\pi M}. \end{align}

Therefore, the discretised boundary integral (4.31)–(4.35) can be written in compact form as follows:

(4.46) \begin{align} \begin{bmatrix} \mathbb I - \dfrac {1-R}{1+R}\mathbb A & 0 & 0 & 0\\[9pt] 0 & \mathbb I + \dfrac {1-R}{1+R}\mathbb A & 0 & 0\\[9pt] 0 & 0 & \mathbb I + \dfrac {1-\varepsilon }{1+\varepsilon }\mathbb A & 0\\[9pt] 0 & 0 & 0 & \mathbb I - \dfrac {1-\varepsilon }{1+\varepsilon }\mathbb A \end{bmatrix} \begin{bmatrix} \mathcal Q^m\\[9pt] \mathcal U^m\\[9pt] \mathcal W^m\\[9pt] \Big (\dfrac {\partial V}{\partial \boldsymbol{t}}\Big )^m \end{bmatrix}\nonumber \\[5pt] = \begin{bmatrix} \dfrac {2R}{1+R}\mathbb A \mathcal J^m - \mathbb B \mathcal T^m\\[11pt] \dfrac {2}{1+R}\mathbb A \mathcal T^m + \mathbb B \mathcal J^m\\[11pt] \dfrac {2}{1+\varepsilon } \mathbb A \mathcal S^m\\[11pt] \dfrac {1}{1+\varepsilon } \mathbb B \mathcal S^m \end{bmatrix} + \begin{bmatrix} -\dfrac {1}{\pi M} \mathcal (T_{\alpha })^m\\[11pt] \dfrac {1}{\pi M} \mathcal (J_{\alpha })^m\\[11pt] 0\\[11pt] \dfrac {1}{\pi M(1+\varepsilon )} \mathcal (S_{\alpha })^m \end{bmatrix}\!, \end{align}

where the superscript $m$ denotes evaluation at $\alpha =\alpha _m, (m = 1,2,{\cdots} ,M-1)$ , $\mathbb I$ is the $M \times M$ identity matrix and $\mathbb A$ and $\mathbb B$ are defined as

(4.47) \begin{align} \mathbb A_{mn} &= \begin{cases} \dfrac {\mathcal A(\alpha _m,\alpha _n)}{M},\,\,\,\,\quad & m \neq n,\\[9pt] \dfrac {s_\alpha (\alpha _m)\cos \big (\theta (\alpha _m)\big )+\theta _\alpha (\alpha _m)}{2M},\,\,\,\,\quad & m = n, \end{cases} \\[-12pt]\nonumber \end{align}

(4.48) \begin{align} \mathbb B_{mn} &= \begin{cases} \dfrac {\mathcal B(\alpha _m,\alpha _n)}{M}, & m \neq n,\\[9pt] \dfrac {s_{\alpha \alpha }(\alpha _m)}{2M s_\alpha (\alpha _m)} + \dfrac {s_\alpha (\alpha _m)\sin \big (\theta (\alpha _m)\big )}{2M}, & m = n. \end{cases} \end{align}

Note that the first two diagonal blocks degenerate to the identity matrix when $R=1$ , reflecting the fact that $\mathcal Q^m$ and $\mathcal U^m$ are already solved in this case. For $\mathcal W^m$ and $(\partial V/\partial \boldsymbol{t})^m$ , we use the GMRES solver with a residual tolerance $10^{-11}$ .

Time integration of (4.1), (4.2), (4.5) and (4.10) is performed with a fourth-order Runge–Kutta method. Higher-order time-marching algorithms have been tested but we found that the fourth-order method balances accuracy and efficiency quite well.

To control aliasing errors we implemented and used both a Krasny filter $\mathcal K_1$ (Krasny Reference Krasny1986b )

(4.49) \begin{align} \mathcal K_1\big (\hat f(k)\big ) = \begin{cases} \hat f(k), & |\hat f(k)| \gt 10^{-11},\\[5pt] 0, & |\hat f(k)| \leqslant 10^{-11} ,\end{cases} \end{align}

and a $36$ th-order exponentially decaying filter $\mathcal K_2$ (Hou & Li Reference Hou and Li2007)

(4.50) \begin{align} \mathcal K_2\big (\hat f(k)\big ) = \hat f(k) \exp \big (-36 (k/k_{\textit{max}})^{36}\big ), \end{align}

where $\hat f(k)$ denotes the Fourier transform of $f(\alpha )$ , and $k_{\textit{max}}$ is the Nyquist wavenumber. Here, $\mathcal K_1$ is applied at every time after performing the FFT, and $\mathcal K_2$ is used at the end of each Runge–Kutta step.

Calculations were typically started with $M = 64$ and $\Delta t = 0.005$ . When the amplitude of the last $20\,\%$ of $\hat \theta (k)$ exceeds $10^{-8}$ , the solution is spectrally interpolated by doubling or quadrupling $M$ . Simultaneously, $\Delta t$ is halved or quartered to prevent aliasing. This process continues until $M$ reaches a predefined value, typically $8192$ .

The Hamiltonian $\mathcal H$ given in § 3 is a conserved quantity, and it was used to monitor our numerical accuracy. An equivalent form of $\mathcal H$ is

(4.51) \begin{align} \mathcal H(t) &= \frac {1}{2}\int _0^1 \varphi _1 \frac {\partial \phi _1}{\partial \boldsymbol{n}}s_\alpha \,\mathrm d\alpha - \frac {U}{2} \int _0^1 \varphi _1 s_\alpha \sin \theta \,\mathrm d\alpha - \frac {R}{2}\int _0^1 \varphi _2 \frac {\partial \phi _2}{\partial \boldsymbol{n}}s_\alpha \,\mathrm d\alpha \nonumber \\ &- \frac {\textit{RU}}{2} \int _0^1 \varphi _2 s_\alpha \sin \theta \,\mathrm d\alpha + \frac {1-R}{2 F_r^2}\int _0^1 \eta ^2s_\alpha \cos \theta \,\mathrm d\alpha \nonumber \\ &+ \frac {E_b}{2} \int _0^1 V \frac {\partial V_1}{\partial \boldsymbol{n}}s_\alpha \,\mathrm d\alpha - \frac {E_b\varepsilon }{2} \int _0^1 V \frac {\partial V_2}{\partial \boldsymbol{n}}s_\alpha \,\mathrm d\alpha . \end{align}

We define its relative error

(4.52) \begin{align} \mathcal H_r(t) = \bigg |\frac {\mathcal H(t) - \mathcal H(0)}{\mathcal H(0)}\bigg |, \end{align}

and monitor its value during the simulation to anticipate loss of accuracy.

5.1. Curvature singularity

We begin with subcritical values of $E_b$ so that the system is linearly unstable and ill posed due to K–H instability. When $E_b = 0$ , we obtain the classical curvature singularity of vortex sheets (Moore Reference Moore1979; Krasny Reference Krasny1986b ; Cowley et al. Reference Cowley, Baker and Tanveer1999), as shown in figure 6 where the initial condition (5.2) was used. The interface and curvature remain odd functions of $x$ during the simulation. They grow and develop the curvature singularity rapidly. Interestingly, the singularity emerges at two symmetrical positions $x = \pm 0.254$ , instead of at $x = 0$ . This is a typical phenomenon when the initial amplitude is moderate, as reported by Krasny (Reference Krasny1986b ) and Cowley et al. (Reference Cowley, Baker and Tanveer1999).

In figure 7, we show the solution with $E_b = 1$ , with initial condition (5.2) as above. The electric field mitigates the K–H instability. The curvature singularity is postponed but eventually develops at $x \approx -0.35$ . Note that the odd symmetry of the interface and curvature is broken. In figure 8 we plot the first- and second-order tangential derivatives of the voltage potential $V$ with respect to the arclength. The behaviour is similar to that shown in figure 2, with $V_s$ and $V_{ss}$ developing a local shock and a singularity, respectively, as $t_s$ is approached.

Figure 6.

Numerical solutions for $E_b = 0$ and $\mu =0.1$ with initial condition (5.2). (a) Interfaces at $t = 0, 1,2,3$ and $3.4925$ . (b) Curvature profiles at $t = 0,1,2,3$ . (c) Curvature at $t = 3.4925$ .

In this section we study the effects of subcritical and supercritical electric fields on curvature singularities. We fix $\varepsilon = 2$ , which implies a critical electric field strength $E_b^{c} = 1.5$ . Initial conditions (5.2)–(5.4) are considered with $\mu = 0.1$ , which give an almost sinusoidal interface with maximum elevation $0.10014$ at $t = 0$ .

Figure 7.

Numerical solutions for $E_b = 1, \varepsilon = 2$ and $\mu =0.1$ with initial condition (5.2). (a) Interfaces at $t = 0,1,2,3,4$ and $4.5$ . (b) Curvature profiles at $t = 0,1,2,3,4$ . (c) Curvature at $t = 4.4975$ .

Figure 7. Long description

The image contains three line graphs labeled (a), (b), and (c). Panel A: A line graph shows the evolution of interfaces at different times. The x-axis is labeled x and ranges from -3 to 3. The y-axis is labeled η and ranges from -0.2 to 0.15. The graph displays multiple curves representing the interface at various times, with an arrow indicating the direction of increasing time. Panel B: A line graph illustrates the curvature profiles at a specific time. The x-axis is labeled x and ranges from -3 to 3. The y-axis is labeled κ and ranges from -0.2 to 0.2. The graph shows multiple curves representing the curvature profiles, with an arrow indicating the direction of increasing time. Panel C: A line graph depicts the curvature at a specific time. The x-axis is labeled x and ranges from -3 to 3. The y-axis is labeled κ and ranges from -10 to 10. The graph shows a single curve representing the curvature.

Figure 8.

Numerical solutions for $E_b = 1, \varepsilon = 2$ and $\mu =0.1$ with initial condition (5.2): (a) $V_s$ at $t = 0,1,2,3,4$ and $4.5$ ; (b) $V_{ss}$ at $t = 0,1,2,3,4$ ; (c) $V_{ss}$ at $t = 4.4975$ .

Figure 9.

Typical curvature amplitude $\Vert \kappa \Vert _\infty$ versus time $t$ for different values of $E_b$ .

For other values $E_b\lt E_b^{c}$ , similar solutions are observed but the time to the singularity increases significantly with $E_b$ . In figure 9, we display typical results for $E_b=0, 0.5, 1, 1.3, 1.5$ that show the evolution of the $L^\infty$ norm of the curvature. Interestingly, when $E_b = E_b^{c}=1.5$ , a similar growing behaviour of the curvature amplitude is observed, even though linear theory indicates that there is no unstable eigenmode.

Figure 10.

Contour plots of (a) $\eta$ and (b) $\theta$ for $E_b = 5$ , $\varepsilon = 2$ and $\mu =0.1$ with initial condition (5.2). The dashed line in (b) is the characteristic line (5.6).

For supercritical fields $E_b\gt E_b^c$ the system becomes neutrally stable and well posed in the linear sense. In figure 10 we set $E_b = 5$ and use the initial condition (5.2). Panels (a) and (b) exhibit contour plots of $\eta$ and $\theta$ , respectively, in the $x{-}t$ plane. The solution resembles standing waves with an approximate temporal period $2\pi /c = 8.22$ , where $c = 0.764$ . After one temporal period, the interface develops some visible distortions, but the amplitude is well bounded by its initial value of $0.10014$ . This demonstrates the stabilising effect of the electric fields. On the other hand, a cusp of $\theta$ is gradually formed and propagates leftward at approximately a speed $-c$ , consistent with the analysis carried out in § 3 that addressed which component develops the curvature singularity first. In panel (b), we plot a typical characteristic line

(5.6) \begin{align} x = -c t + 2\pi , \end{align}

which captures well the trajectory of the cusp in the $x-t$ plane.

In figure 11, we compare the results of the full BIM computation (solid curves) with the model (3.13) (dots). Panel (a) depicts the interface amplitude of $\eta$ by evaluating its $L^\infty$ norm $\Vert \eta \Vert _\infty$ (this is clearly bounded), and panel (b) shows the spatial variation of $\eta$ at the final computed time $t = 13.88$ . Panel (c) tracks the evolution of the curvature amplitude $\Vert \kappa \Vert _\infty$ , which is seen to exhibit oscillatory and growing behaviour, and panel (d) shows $\kappa$ at the final computed time $t = 13.88$ , showing the emergence of a curvature singularity. All results demonstrate excellent agreement between the weakly nonlinear model and the fully nonlinear solutions of the Euler-electrostatics equations. Panel (e) shows the relative error $\mathcal H_r$ of the Hamiltonian (see (4.52) for its definition), which oscillates near round-off level for most of the simulation. Results of the BIM and the model still match well, although a small discrepancy emerges as a singularity is encountered. The structure of the singularity is investigated in panel (e) where the spectrum of $\theta$ is plotted at successive times as the singularity is approached. The red line indicates an algebraically decaying rate proportional to $k^{-1.5}$ , a typical result of the vortex-sheet curvature singularity (Moore Reference Moore1979; Krasny Reference Krasny1986b ; Cowley et al. Reference Cowley, Baker and Tanveer1999).

Figure 11.

Numerical results for $E_b = 5$ , $\varepsilon = 2$ and $\mu =0.1$ with initial condition (5.2). (a) Interface amplitude $\Vert \eta \Vert _\infty$ versus $t$ . (b) Value of $\eta$ at $t = 13.88$ . (c) Curvature amplitude $\Vert \kappa \Vert _\infty$ versus time $t$ . (d) Local blow-up of $\kappa$ at $t = 13.88$ . (d) Relative error of the Hamiltonian $\mathcal H_r$ . (e) Spectrum of $\theta$ at $t = 4,8,12,13.6$ and $13.88$ from left to right. The red line is proportional to $k^{-1.5}$ . In (a–c), solid lines and dots are results of the BIM and model (3.13).

Figure 12 shows the results for $E_b = 5$ and the initial condition (5.4). Panel (a) displays four interfacial solutions at increasing times. The profiles are shifted upwards for comparison. The solid curves, circles and asterisks are results of the BIM, model (3.13) and model (3.52), respectively. The agreement is excellent, with no visible discrepancy. In panel (b), we plot the profiles of $\theta$ at $t = 14.45$ , just before a cusp develops. The discrepancy between the three solutions can now be clearly seen. As explained in § 3.2, given an odd initial condition for $\eta$ , model (3.52) preserves the odd symmetry with respect to the vertical line $x-ct = 0$ . Therefore, the corresponding $\theta$ displays even symmetry. On the other hand, neither (3.13) nor the full Euler-electrostatics equations have this property. This discrepancy is due to the absence of high-order corrections in (3.52). Nevertheless, (3.52) still provides an excellent estimate of the singularity time $t_s$ , as will be shown in the next subsection. Panel (c) plots the curvature amplitude $\Vert \kappa \Vert _\infty$ versus time $t$ for the three solutions, exhibiting finite-time blow-up behaviour in all cases. In panel (d) we plot the relative error $\mathcal H_r$ of the Hamiltonian calculated from the BIM, and we see that it is smaller than $2\times 10^{-11}$ throughout the simulation.

5.2. Singularity time

In what follows we display some typical numerical results on the singularity time $t_s$ and compare them with the asymptotic estimates obtained in § 3. Note that for the linearised system, the initial conditions (5.3) and (5.4) correspond to a temporally growing eigenmode for subcritical electric fields and a right-travelling eigenmode for supercritical electric fields. They align with the basic assumptions that are used to derive models (3.24) and (3.52). On the other hand, the initial condition (5.2) represents a mixture of growing and decaying eigenmodes, or left-travelling and right-travelling eigenmodes. Therefore, the singularity time is dominated by one of these two modes whose initial amplitudes are halved, and the asymptotic estimations remain valid if we replace $\delta$ by $\delta /2$ , as observed earlier. For convenience, we summarise the theoretical model predictions from § 3 in table 1.

Table 1.

Theoretical predictions of the singular time $t_s$ and the asymptotic behaviours for the three initial conditions and different ranges of $E_b$ .

Table 1. Long description

A table with three columns and four rows, including a header row. The columns are labeled Range of Eb, Initial conditions, and Asymptotic behaviours. The rows compare different ranges of Eb and their corresponding initial conditions and asymptotic behaviors. Row 1: Range of Eb: (5.2), Initial conditions: (5.3), Asymptotic behaviours: delta to 0. Row 2: Range of Eb: Eb less than Eb sub c, Initial conditions: (3.43), Asymptotic behaviours: t s proportional to ln(1/delta). Row 3: Range of Eb: Eb greater than Eb sub c, Initial conditions: (3.62), Asymptotic behaviours: t s proportional to 1/delta. Row 4: Range of Eb: (5.4), Initial conditions: (3.41), Asymptotic behaviours: t s proportional to Eb sub b to the power of negative 1/2.

Figure 12.

Results for $E_b = 5$ and the initial condition (5.4) in the frame of reference moving at a speed $c = 0.7638$ . (a) Interfaces at $t = 0, 5, 10$ and $14.45$ . (b) Profiles of $\theta$ at $t = 14.45$ . (c) Curvature amplitude $\Vert \kappa \Vert _\infty$ versus time $t$ . (d) Relative error $\mathcal H_r$ of the Hamiltonian. Solid curves, circles and asterisks represent solutions of BIM, (3.13) and (3.52).

Figure 12. Long description

Panel A: A line graph shows the evolution of the interface over time. The horizontal axis is labeled x minus c t and ranges from negative 3 to positive 3. The vertical axis is labeled eta and ranges from negative 0.10 to positive 0.25. The graph includes multiple lines representing different time points, with an arrow indicating the direction of increasing time. Panel B: A line graph displays the profiles of theta at different times. The horizontal axis is labeled x minus c t and ranges from negative 1.5 to positive 4.5. The vertical axis is labeled theta and ranges from negative 0.10 to positive 0.30. An inset zooms in on the region from negative 0.3 to positive 0.3 on the horizontal axis. Panel C: A line graph illustrates the curvature amplitude versus time. The horizontal axis is labeled t and ranges from 0 to 15. The vertical axis is labeled with a mathematical symbol and ranges from 0.1 to 100 on a logarithmic scale. Panel D: A line graph shows the relative error of the Hamiltonian over time. The horizontal axis is labeled t and ranges from 0 to 15. The vertical axis is labeled H subscript r and ranges from 10 to the power of negative 15 to 10 to the power of negative 11 on a logarithmic scale.

Figure 13.

Value of $t_s$ versus $E_b$ obtained from the BIM (circles) and asymptotic estimates (curves) for $\varepsilon =0.5$ . The initial conditions are those matching the theoretical predictions in table 1 with $\mu = 0.1$ . The black dashed line identifies the critical electric field strength $E_b^{c} = 3$ .

Figure 13 displays some typical results of the variation of $t_s$ with $E_b$ for $\varepsilon = 0.5$ . For the parameters used the critical electric field strength is $E_b^{c} = 3$ . Solid curves show results for the four asymptotic predictions represented by their labels, and circles are results obtained from the full BIM simulations with the matched initial conditions shown in table 1. The numerical results show that for subcritical values of the electric field, the singularity time $t_s$ increases monotonically with $E_b$ . This behaviour is captured well by the asymptotic estimates (3.41) and (3.43), except in a narrow region adjacent to $E_b = E_b^{c}=3$ where a drastic decrease of $t_s$ is predicted. As mentioned previously, this is due to the underlying model (3.24) becoming non-uniformly valid near $E_b^c$ . The numerical results of the BIM simulations (open circles) show that for initial condition (5.2), $t_s$ smoothly increases and passes through the critical line $E_b = E_b^{c}$ . For supercritical electric fields, (3.63) yields that $t_s$ reaches a finite maximum value at $E_b = (1+\varepsilon )/(1-\varepsilon )^2 = 6$ , which perfectly matches the numerical results. When $E_b\gt 6$ , $t_s$ slowly decreases with $E_b$ , following the theoretical prediction closely.

Figure 14.

Comparisons of $t_s$ obtained from the BIM (circles) and asymptotic estimates (solid curves) for $\varepsilon =2$ . The initial condition for the BIM is (5.4) with $\mu = 0.1$ . The black dashed lines label the two roots, $E_b = 1.544$ and $52.46$ , of $\gamma = 0$ . The corresponding critical electric field strength $E_b^c = 1.5$ .

When $\varepsilon \gt 1$ , $t_s$ exhibits a totally different behaviour. Figure 14 shows a typical result for $\varepsilon = 2$ . The circles represent the results of the BIM simulations for initial condition (5.4) with $\mu = 0.1$ , and the solid curves correspond to the asymptotic estimate (3.59). Following (3.63), $|\gamma |^{-1}$ blows up at $E_b = 1.544$ and $E_b = 52.46$ , suggesting that the curvature singularity is significantly postponed in parameter regions adjacent to these two values. This explains well the observed local peak of $t_s$ around $E_b = 1.544$ , even though the full BIM simulation predicts a curvature singularity for $E_b=1.544$ at $t_s\approx 19.4$ . As discussed earlier, the predicted unbounded growth $t_s\to \infty$ by the model (3.24) reflects its breakdown in the vicinity of $E_b = 1.544$ . However, it is interesting to note that for solutions centred around $E_b = 52.46$ , the corresponding $t_s$ for the full BIM simulations does indeed increase rapidly and follows the theoretical prediction closely. We have not included data for $t_s$ for the full BIM simulations in the interval $E_b\in (35,80)$ , because the curvature singularity develops very slowly and makes the precise estimate of $t_s$ challenging. It is worth mentioning that given small-amplitude initial conditions, solutions for $E_b = 52.46$ are rather stable. Using only $64$ Fourier modes, we ran the fully nonlinear simulations up to $t = 1000$ for $\mu = 0.1$ and $0.2$ with initial condition (5.4). The spectrum of $\theta$ for $\mu = 0.2$ broadens, with the largest amplitude among the last $20\,\%$ Fourier modes at $t=1000$ falling below $10^{-10}$ , and the amplitude for $\mu = 0.1$ being of the order of round-off level. All solutions exhibit perfect travelling-wave forms. Results for $\mu = 0.1$ are shown in figure 15 where the profiles of the interface and the curvature at $t = 0$ and $1000$ (the latter profile is shifted along the $x$ -direction for comparison) are compared.

Figure 15.

Numerical solutions for $E_b = 52.46$ and $\varepsilon = 2$ at $t = 0$ (curve) and $t = 1000$ (circles). (a) Interface $\eta$ . (b) Curvature $\kappa$ . The solution at $t = 1000$ is shifted along the $x$ -direction for comparison.

In the limit $E_b\rightarrow \infty$ , the asymptotic estimates (3.59) and (3.62) predict that $t_s\propto E_b^{-1/2}$ . Figure 16 presents some typical results for $\varepsilon =0.5$ employing a $\log$ – $\log$ scale, and using initial condition (5.4) with $\mu = 0.1$ . As $E_b$ increases up to $10^4$ , the computed estimates of $t_s$ show excellent agreement with the theoretical prediction. The numerical results converge to the predicted power law scaling $E_b^{-1/2}$ scaling, which is indicated by the red line.

Figure 16.

Value of $t_s$ versus $E_b$ in a $\log$ – $\log$ scale. Circles and curves are results of the BIM and the asymptotic estimation (3.59) for $\varepsilon = 0.5$ and $\mu = 0.1$ with initial condition (5.4). The red line is proportional to $E_b^{-1/2}$ .

Finally, we examine the dependence of $t_s$ on the initial amplitude. When $\mu \ll 1$ , the initial conditions (5.2)–(5.4) give rise to almost sinusoidal interface profiles with amplitudes nearly equal to $\mu$ . For instance, $\mu = 0.1$ corresponds to an almost sinusoidal interface with an amplitude of $0.10014$ . Therefore, we directly compare our numerical results and the asymptotic estimations without distinguishing $\mu$ and $\delta$ . In figure 17 we present the numerical results of the BIM simulations and the three model (3.13), (3.24) and (3.52). The results for each are represented by the circles, plus signs, triangles and crosses, respectively. The solid curves are theoretical predictions corresponding to their labels. The initial conditions are those matched in table 1, and they are numerically transformed into the initial conditions for $\eta$ and $\varphi$ when solving the model equations, as explained in § 5. In figure 17(a) we set $E_b = 1$ and $\varepsilon = 2$ (subcritical electric field). The $t_s$ values obtained from the BIM simulations and the various model equations agree well and follow the theoretical predictions consistently. Figure 17(b) displays results for $E_b = 5$ and $\varepsilon = 0.5$ (supercritical TEF) using a $\log$ – $\log$ scale. The numerical results exhibit a clear inverse proportionality with $\delta$ , and match the theoretical predictions precisely when $\delta \lt 0.1$ . For moderate values of $\delta$ , nonlinear effects cause $t_s$ to slightly deviate from the theoretical predictions, with an underestimation for (3.59) and an overestimation of (3.62).

Figure 17.

Singular time $t_{s}$ versus $\delta$ for $E_b = 1, \varepsilon = 2$ (a), and $E_b = 5, \varepsilon = 0.5$ (b). The circles, plus signs, triangles and crosses are numerical results obtained from the BIM simulations, (3.13), (3.24) and (3.52), respectively. The solid curves are theoretical estimates represented by their labels. The initial conditions are those matching the theoretical predictions in table 1.