Instability of a planar fluid interface under a tangential electric field in a stagnation point flow

Abstract The interface between two immiscible fluids can become unstable under the effect of an imposed tangential electric field along with a stagnation point flow. This canonical situation, which arises in a wide range of electrohydrodynamic systems including at the equator of electrified droplets, can result in unstable interface deflections where the perturbed interface gets drawn along the extensional axis of the flow while experiencing strong charge build-up. Here, we present analytical and numerical analyses of the stability of a planar interface separating two immiscible fluid layers subject to a tangential electric field and a stagnation point flow. The interfacial charge dynamics is captured by a conservation equation accounting for Ohmic conduction, advection by the flow and finite charge relaxation. Using this model, we perform a local linear stability analysis in the vicinity of the stagnation point to study the behaviour of the system in terms of the relevant dimensionless groups of the problem. The local theory is complemented with a numerical normal-mode linear stability analysis based on the full system of equations and boundary conditions using the boundary element method. Our analysis demonstrates the subtle interplay of charge convection and conduction in the dynamics of the system, which oppose one another in the dominant unstable eigenmode. Finally, numerical simulations of the full nonlinear problem demonstrate how the coupling of flow and interfacial charge dynamics can give rise to nonlinear phenomena such as tip formation and the growth of charge density shocks.

Interfaces polarize in applied electric fields. Free charges brought by conduction accumulate at the boundary between phases and the electric field acting on this induced charge creates shear stresses that drag the fluids into motion (Melcher & Taylor 1969). In the case of a drop in a uniform electric field, the classic small-deformation analysis by Taylor (1966) showed that the resulting EHD flow consists of two toroidal vortices inside and a stresslet-quadrupole flow outside the drop. Depending on the electric properties of the fluids, the surface flow is directed either to the poles or to the equator. The latter case is shown in figure 1(a). In strong fields, however, this flow undergoes a plethora of instabilities that may result in drop breakup (Taylor 1964;Torza, Cox & Mason 1971;Sherwood 1988;Lac & Homsy 2007;Karyappa, Deshmukh & Thaokar 2014). Droplet disintegration can proceed in various patterns depending on fluid properties. In the case of the pole-convergent flow, the drop can develop conical tips that emit jets, which subsequently break up into droplets (Collins et al. 2008(Collins et al. , 2013Mohamed et al. 2016;Sengupta, Walker & Khair 2017). In the case of the equator-converging flow, the drop either dimples at the poles and breaks into a torus, or deforms into a pancake-like lenticular shape with a sharp edge emitting rings encircling the drop (Brosseau & Vlahovska 2017;Wagoner et al. 2020).
While EHD streaming from Taylor cones has been extensively studied (Collins et al. 2008(Collins et al. , 2013Herrada et al. 2012;Ganan-Calvo et al. 2016), the mechanisms underlying the equatorial streaming remain an open question. Noting the similarity between the EHD tip streaming and the tip streaming in flow focusing (Barrero & Loscertales 2007;Ganan-Calvo & Montanero 2009;Anna 2016), Brosseau & Vlahovska (2017) speculated, in the original paper that reported the phenomenon, that the EHD equatorial streaming arises from an interfacial instability due to a convergent flow (Pozrikidis & Blyth 2004;Blyth & Pozrikidis 2005;Tseng & Prosperetti 2015). Near a stagnation point, a perturbation of the interface may get drawn by the flow and grow into a fluid filament if viscous stresses overcome interfacial tension. Unlike flow focusing, however, EHD streaming involves both flow and electric field. In EHD tip streaming, the electric field is initially normal to the interface at the stagnation point, while in EHD equatorial streaming, the applied field is parallel to the interface at the stagnation line.
In this work, we analyse the effect of an electric field on the convergent flow instability in a configuration mimicking the EHD equatorial streaming as depicted in figure 1(b). We develop a two-dimensional model to study the dynamics of a system of two superimposed layers of fluids subject to a tangential electric field and a stagnation point flow. The convergent flow and the electric field are assumed to be independently applied, unlike the equatorial EHD instability, where the flow is generated by the electric field. Despite this simplification, the analysis provides valuable insights into mechanisms responsible for the EHD equatorial streaming such as the evolution of the convergent line instability and the emergence of charge shocks. We present the governing equations in § 2 and their non-dimensionalization in § 3. We develop a local linear stability theory in § 4. To supplement our theory, we employ the boundary element method, outlined in § 5, to perform numerical simulations as well as a numerical linear stability analysis. We compare the results from the local theory, numerical linear stability analysis and transient simulations in § 6. Finally, we conclude and discuss potential extensions of this work in § 7.

Problem definition and governing equations
We study EHD instabilities that arise at the interface S between two immiscible fluids under the combined effects of a tangential electric field E 0 and of an imposed stagnation point flow u ∞ (x), to be specified more precisely later. The two layers are labelled 1 and 2 as depicted in figure 2, with fluid 1 occupying the lower half-space. The applied electric field is uniform along the x direction, and the stagnation point is located on the interface at the origin O of the coordinate system. At equilibrium, the interface is uncharged and flat and coincides with the plane z = 0, and we consider two-dimensional dynamics in the (x, z) plane. The shape of the deformed interface is parametrized as z = ξ(x, t), and has unit normal n pointing from fluid 1 into fluid 2.
Both fluids are leaky dielectrics with electric conductivity σ , electric permittivity and dynamic viscosity μ. While we formulate the governing equations to allow for distinct viscosities, all the results presented in § 6 are for equiviscous fluids. Under the Taylor-Melcher leaky dielectric model (Melcher & Taylor 1969), any net charge in the system occurs at the location of the interface while the bulk of the fluids remains electroneutral. The Taylor-Melcher leaky dielectric model can be formally derived from more detailed electrokinetic models based on the Poisson-Nernst-Planck equations in the limit of strong electric fields and under the assumption of thin Debye layers (Schnitzer & Yariv 2015;Mori & Young 2018). As an example, these assumptions are valid for millimetre-sized drops of leaky dielectric liquids subject to electric fields of magnitude E 0 ∼ 10 3 V cm −1 according to Saville (1997). Under these assumptions, the electric potential is harmonic in both fluids: (2.1) Far from the interface, the electric field E = −∇ϕ tends to the applied uniform field: Across the interface, its tangential component is continuous while a jump can arise in the normal direction due to the mismatch in electrical properties: We have introduced the notation [[F ]] = F 2 − F 1 for the jump of any variable F defined on both sides of the interface. A surface charge density develops at the interface following Gauss's law: This surface charge density satisfies a conservation equation accounting for finite charge relaxation, Ohmic conduction from the bulk and charge convection by the flow: where ∇ s = (I − nn) · ∇ is the surface gradient operator and u is the total fluid velocity. Neglecting the effects of inertia and gravity, the velocity field and its corresponding pressure field satisfy the Stokes equations in both fluids: The velocity vector is continuous across the interface and approaches the applied flow far from the interface: and u(x) → u ∞ (x), as z → ±∞. (2.8) In the absence of Marangoni effects, the jump in hydrodynamic and electric tractions across the interface is balanced by capillary forces: with uniform surface tension γ . Hydrodynamic and electric tractions are expressed in terms of the Newtonian and Maxwell stress tensors, respectively: (2.11) Finally, the interface evolves and deforms under the local velocity field as a material surface. Defining the function g(x, t) = z − ξ(x, t), the kinematic boundary condition reads Dg (2.12) leading to the condition (u, v) are the velocity components. We also note that the surface normal is given by n = ∇g/|∇g|. Our focus is on analysing the stability of the interface near a stagnation point, and to this end we choose the background flow u ∞ = (u ∞ , v ∞ ) to be extensional along the z direction, compressional along the x direction and with a stagnation point at the origin. The strength of the background flow is characterized by the local strain rate at the stagnation point, which is denoted by A = ∂ z v ∞ , where A > 0. Three different types of background flows are used in this study and are illustrated in figure 3. The first type, depicted in figure 3(a), is simply the linear flow u ∞ = (−Ax, Az), which is used in the development of the local linear theory of § 4. The numerical scheme of § 5, however, requires periodicity in the x direction, and for this reason we also consider two periodic flow fields. The first one is the Taylor-Green vortex (Taylor & Green 1937) shown in figure 3(b) and defined as (2.16) The local strain rate at the stagnation point is then simply given by We also consider a second background flow shown in figure 3(c) and induced by a periodic array of anti-parallel point forces separated by distances L p and d along the x and z directions, respectively: where m = mê z is the strength of the point forces, x u 0 and x l 0 are the locations of upper and lower point forces, respectively, and G P is the singly periodic Green's function for the Stokes equations (Pozrikidis 1992). The local strain rate at the stagnation point is given by where k p = 2π/L p is the wavenumber associated with the unit cell. In obtaining (2.19), it is assumed that x u 0 = x l 0 = 0 and z u 0 = −z l 0 = d/2. Our numerical calculations produce identical linear stability results for a given local strain rate under the two periodic flow fields described by (2.15) and (2.18). In all transient nonlinear simulations, we use the periodic array of point forces as background flow.

Non-dimensionalization
For the system presented above, dimensional analysis yields six non-dimensional groups, three of which characterize the mismatch in physical properties between the two layers: The other three parameters can be obtained as ratios of characteristic time scales in the problem. First, we note that free charges in the bulk fluids relax on a conduction time scale defined as τ c = σ . (3.2) Note that the product RQ = τ c,1 /τ c,2 is the ratio of the charge relaxation time scales in the two liquid phases, and characterizes their responses to conduction. For instance, RQ > 1 when the lower layer is less conductive. The time scale for the deformed interface to relax to its flat configuration under capillary effects can be defined as where k −1 is the characteristic length scale associated with the deformation. In our periodic simulations and analysis, we use the length scale k −1 p = L p /2π based on the periodicity of the base flow, whereas k is the wavenumber of the plane wave in the local stability theory of § 4. The accumulation of free charges on the interface creates an electric force that can drive the fluid into motion. The time scale for deformations under this EHD flow is inversely proportional to the magnitude of electric tractions on the interface: Finally, the characteristic time scale associated with the background flow is given by the inverse of the applied strain rate at the stagnation point: Taking ratios of these time scales yields the three remaining dimensionless groups, which we define as The electric capillary number Ca E compares electric with capillary forces, while the electric Reynolds number Re E characterizes the importance of charge convection versus conduction, the two mechanisms responsible for the evolution of the interfacial charge distribution. Finally,Â is the dimensionless strain rate of the applied flow. We scale the governing equations and boundary conditions using time scale τ c,2 , length scale k −1 , pressure scale 2 E 2 0 , velocity scale (τ EHD k) −1 and characteristic electric potential E 0 k −1 . The dimensionless Stokes equations read The charge conservation equation (2.5) becomes The stress balance at the interface is written Finally, the kinematic boundary condition becomes (3.12) The remaining governing equations and boundary conditions in (2.1)-(2.3), (2.7) and (2.8) remain unchanged in their non-dimensional form, and hence are not repeated here. In the remainder of the paper, all equations and variables are presented in non-dimensional form.

Local linear stability theory
We first perform a linear stability analysis in the vicinity of the stagnation point, in the spirit of Tseng & Prosperetti (2015) who, following an approach previously proposed by Pozrikidis & Blyth (2004), considered a similar problem with inertia but no electric field. In the base state (subscript 0), the interface is flat and uncharged: ξ 0 (x) = 0, q 0 (x) = 0. The applied electric field generates a potential ϕ 0 (x, z) = −x and results in a pressure jump [[p 0 ]] = (Q − 1)/2 across the interface. The base flow is taken to be the planar extensional flow The base state variables are perturbed as We substitute these expressions into the governing equations and boundary conditions and linearize with respect to ε. Following Tseng & Prosperetti (2015), we neglect all terms in the linearization that have non-constant coefficients, which restricts our analysis to the neighbourhood of the stagnation point at (0, 0); the consequences of this approximation are discussed in § 6. The governing equations for the potential and streamfunction in the two regions are while the kinematic and dynamic boundary conditions yield Equations (4.2a,b)-(4.7) form a system of homogeneous constant-coefficient linear partial differential equations. Recall from § 3 that in the present non-dimensionalization lengths have been scaled by k −1 , where k is the wavenumber of the perturbation. We therefore seek normal-mode solutions of the formφ(x, z, t) =φ(z) exp(ix + st), with similar expressions for all the variables. Equation (4.2a,b), along with the decay properties as z → ±∞, leads toφ Applying the boundary conditions yields an algebraic system for the unknown coefficients. Setting its determinant to zero provides the dispersion relation for the growth rate s: where the dependence on wavenumber k is implicit via the electric capillary number defined in (3.6a-c). The first term on the right-hand side of (4.9) shows that the background flow is destabilizing whenÂ > 0, i.e. when the interface is aligned with the compressional axis (Tseng & Prosperetti 2015). The second and third terms describe the effects of capillary and electric stresses, respectively. A more detailed discussion of this dispersion relation is deferred to § 6.

Boundary element method and numerical stability
We complement the local linear analysis of § 4 with numerical simulations and a numerical stability analysis. We first present in § 5.1 a numerical method for the nonlinear solution of the system of governing equations (3.7)-(3.11) based on the boundary integral equations for the Laplace and Stokes equations in a periodic domain of period L p along the x direction. These simulations will provide insight into the dynamics of the system far from its base state. The methodology is similar to that of Firouznia & Saintillan (2021), and implements adaptive grid refinement to handle large local deformations and charge gradients in the nonlinear regime of growth. Subsequently in § 5.2, we utilize the same boundary element method to perform a numerical normal-mode linear stability analysis by computing the Jacobian of the dynamical system and solving for its eigenspectrum to identify fundamental modes of instability.

Boundary element method
We formulate the electric problem using the boundary integral equation for Laplace's equation (Sherwood 1988;Baygents, Rivette & Stone 1998;Lac & Homsy 2007): where the evaluation point x 0 can be anywhere in space whereas x denotes the integration point on the interface. The periodic Green's function for Laplace's equation, G P (x 0 ; x), represents the potential due to a periodic array of point sources with period L p along the x axis (Pozrikidis 2002). Taking the gradient of (5.1) with respect to x 0 and using Gauss's law (3.10), we can derive an integral equation for the jump in the normal electric field across the interface: (5.2) Given the charge distribution q, (5.2) can be used to solve for [[E n ]], from which we obtain E n 1 and E n 2 as The tangential electric field E t = −∇ s ϕ can then also be obtained by differentiating the electric potential in (5.1) in the direction tangential to the interface. Similarly, the flow problem can be formulated in boundary integral form as (Rallison & Acrivos 1978;Pozrikidis 1992) where the hydrodynamic traction jump [[ f H ]] is obtained from the dynamic boundary condition (2.9). Here, G P is the singly periodic Green's function capturing the flow due to a periodic array of point forces separated by the distance L p along the x direction, and T P is the corresponding stress tensor (Pozrikidis 1992). Note that, in the results shown below, we choose λ = 1 and therefore the double-layer potential vanishes in (5.4). The single-layer potential exhibits a logarithmic singularity when x approaches x 0 , which we isolate and treat separately with a quadrature designed for singular integrands (Stroud & Secrest 1966). Gauss-Legendre quadrature with six base points is used for non-singular elements.
The numerical algorithm for transient nonlinear simulations follows Firouznia & Saintillan (2021) and Das & Saintillan (2017b) and can be summarized as follows. At t = 0, the periodic flat interface is discretized into N elements using N grid points with locations x i (t) that move with the normal component of the interfacial velocity to satisfy the kinematic boundary condition. The interface shape is reconstructed using cubic splines based on the grid point locations, which allows for an accurate and convenient determination of geometric properties such as the normal and tangential vectors and surface curvature, and for the accurate evaluation of surface integrals. Considering the high order of accuracy of cubic spline interpolations with a reasonable grid resolution, the most challenging errors are those incurred in the numerical integration, especially in the treatment of the singular terms. The asymptotic rate of convergence of our method is found to be between 1.5 and 2. More details on error analysis are presented in § 6.1.
At every time iteration, we perform the following steps: (i) Given the current charge distribution q(x) and shape of the interface, compute [[E n (x)]] by numerically inverting (5.2) using a GMRES solver (Saad & Schultz 1986;Frayssé et al. 2005). From [[E n (x)]], obtain E n 1 and E n 2 via (5.3a,b). (ii) Determine the potential ϕ along the interface by evaluating (5.1). (iii) Differentiate the surface potential numerically along the interface in order to obtain the tangential electric field E t = −∇ s ϕ.

Numerical linear stability analysis
We also perform a numerical normal-mode linear stability analysis based on the full system of equations and boundary conditions, following a method proposed by Pozrikidis (2012) for the stability of pendant drops. We analyse the stability of the base state of a flat interface with zero charge, which allows us to parametrize the interface as z = ξ(x, t). The system of governing equations can be viewed as a dynamical system for the surface charge density q(x, t) and interface deflection ξ(x, t), which evolve according to (3.9) and (3.12), where the electric field E and velocity u are solutions of the boundary integral equations (5.2) and (5.4). The right-hand side in (5.5) is evaluated on the interface. It can be viewed as a nonlinear functional of the two variables (q, ξ) and can be calculated numerically using the algorithm of § 5.1. The flat configuration with zero charge, given by q(x, t) = ξ(x, t) = 0, is an equilibrium solution.
The numerical linear stability analysis is performed in a periodic domain of period L p . After spatial discretization of the unit period using N grid points, the dynamical system (5.5) yields a system of coupled ordinary differential equations of the form where Y and J are vectors of length 2N containing the values of the variables at the grid points: The linear stability of the equilibrium solution Y = 0 is studied by perturbing the system as Y (t) = εŶ e st . At linear order in ε 1, we obtain a linear eigenvalue problem is the Jacobian of the system. The components of the Jacobian are calculated numerically using a second-order central finite difference scheme: where each variable Y i is successively perturbed by a small amount ±δY i (corresponding to a small perturbation of charge ±δq for i = 1, . . . , N or of shape ±δξ for i = N + 1, . . . , 2N), and where J k (±δY i ) at the numerator is obtained using the boundary integral method. Once J is known, its eigenvalues s provide the growth rates of the perturbation, while its eigenvectorsŶ capture the corresponding eigenmodes of charge and shape.

Results and discussion
We present results on the stability of the system by comparing predictions from the local linear theory (LT) of § 4 and from the numerical linear stability analysis (Num-LSA) of § 5.2. Nonlinear dynamics is also explored in transient simulations (TS) using the boundary element method of § 5.1. We discuss our results in the following order. First, we analyse the behaviour of the system subject to a tangential electric field in the absence of any background flow in § 6.1. Next, in § 6.2, we study the effect of an extensional background flow when there is no electric field. Finally, we characterize in § 6.3 the interplay between the electric field and the flow when both are applied to the system simultaneously.

Effect of tangential electric field
We first consider the case where a tangential electric field is applied to the interface in the absence of any background flow. In this case, results from the LT and Num-LSA are expected to match, as the approximations of the LT only affect terms involving the applied flow. Since the base flow is stationary and the interfacial charge is zero in the equilibrium state, the effects of charge convection only arise at quadratic order and therefore have no effect on the linear stability. The growth rate predicted by LT in this case is given by Note that the dependence on wavenumber k is through the electric capillary number, with Ca −1 E ∝ k. Equation (6.1) is consistent with the results of Melcher & Schwarz (1968) in the limit of zero inertia.
In the limit of instantaneous charge relaxation (i.e. considering only Ohmic terms in (2.5)), the growth rate further simplifies to The first term on the right-hand side is always negative and captures the stabilizing effect of capillary stresses. It is proportional to k, indicating that surface tension preferentially stabilizes high wavenumbers. The last term in (6.2) captures the effect of electric stresses and can be of either sign. For the system to be electrically unstable, the following condition must be met: which means either R > 1, Q > 1 or R < 1, Q < 1. Setting s = 0 in (6.2) also provides a critical electric capillary number for instability: The system is unstable for Ca E ≥ Ca E,c (long waves), and it is stable otherwise. The maximum growth rate is reached at zero wavenumber or under vanishing surface tension (Ca E → ∞) and is given by In the case of finite charge relaxation, the dispersion relation (6.1) is a quadratic equation for the growth rate s. The roots can be shown to be imaginary only when which is incompatible with the condition of (6.3) for instability, so that the growth rate is always real in electrically unstable systems. Figure 4(a,b) shows the dominant unstable modes of deformation and charge distribution obtained via Num-LSA. It is evident that all the modes are sinusoidal as expected in the absence of flow, and the fastest-growing mode,ξ 1 , has the largest possible wavelength permitted by the computational domain. This is indeed expected based on LT. For comparison, we also calculate the growth rate of various eigenmodes numerically by performing short-time TS. The growth rate s obtained from LT, Num-LSA and TS  The nonlinear evolution of the interface shape and charge distribution is illustrated in figure 4(c, d) (also see supplementary movie available at https://doi.org/10.1017/jfm. 2021.967), showing a representative TS where the interface shape and charge distribution were initially perturbed by the dominant unstable eigenmode with a small amplitude. At short times, the sinusoidal modes amplify as the surface deflection grows and charge is brought to the interface via Ohmic conduction. As nonlinear effects become significant, the interface deflection becomes asymmetric. Electric stresses on the interface drive a flow which tends to further sweep opposite charges towards the interface tip, leading to the development of sharp charge gradients and of a pointed tip with high curvature. The tip grows unboundedly with an increasing curvature until it eventually causes our numerical method to break down. One should note that the eigenmodes have up-down mirror symmetry, meaning that both (ξ,q) and (−ξ, −q) have identical growth rates and exhibit similar dynamics at short times. This is in contrast to the nonlinear regime of evolution where the shape and charge distributions become asymmetric as evident in fate of the system is entirely determined by the balance of viscous and capillary stresses. Consequently, the only two time scales in the problem are τ γ and τ f , previously defined in (3.3) and (3.5), respectively. The LT yields the following expression for the growth rate: where s * LT has been scaled by τ −1 f instead of τ −1 c,2 , and where Ca = τ γ /τ f is the viscous capillary number and remains proportional to k −1 p . This result is consistent with the analysis of Tseng & Prosperetti (2015) in the limit of zero inertia. Figure 6(a) shows the most unstable modes of deformation obtained via Num-LSA. The modes are clearly non-sinusoidal in this case, with deflections from the flat base state occurring primarily in the neighbourhood of the stagnation point. The dominant modeξ 1 resembles a Gaussian centred around the origin, and higher-order modes involve shapes with increasing numbers of oscillations, all concentrated near x = 0. Since the modes are non-sinusoidal, the growth rate of the fastest-growing mode differs from the local prediction of (6.7). This is confirmed in figure 6(c) where the growth rates from LT and Num-LSA are compared as a function of Ca. Both methods provide similar growth rates at high capillary numbers (long wavelengths), but their predictions diverge at small values of Ca: while the local theory shows a stabilization of the system below a critical capillary number, the numerical stability analysis predicts that the system is always unstable.
The nonlinear evolution of the interface shape is illustrated in figure 6(b) (also see supplementary movie), showing a transient simulation in which the interface shape was perturbed at t = 0 by the dominant unstable eigenmode of shape. The interface deflection increases with time, and as nonlinear effects become significant the dimple in the interface narrows while the curvature at the tip increases, leading to an increase in capillary stresses which tend to resist further deformation. As shown in the inset of figure 6(b), this causes the interface deflection to saturate and reach a steady profile where capillary stresses balance viscous stresses arising from the applied flow. This is unlike the case of figure 4 for the electric field only, where the tip deformation did not saturate.

Combined effects of electric field and flow
We now turn to the general case where the system is subject to both a tangential electric field and a stagnation point flow. In this case, the LT with finite charge relaxation and charge convection by the flow yields the dispersion relation of (4.9), with sinusoidal eigenmodes for all perturbation variables. It is clear, from the form of (4.9), that the applied flow and electric field contribute additively to the growth rate: the presence of the base flow simply shifts the growth rate of (6.2) for the electric problem by an amount ofÂ. In particular, an external flow withÂ > 0 always has a destabilizing effect under the local theory approximations. If charge convection is neglected in the theory, the effect of the background flow only affects the dynamics of the system through the kinematic boundary condition and the dispersion relation reduces to As we show in figure 8 and discuss further below, this approximation results in a decrease in growth rate when compared with (4.9), suggesting that charge convection is destabilizing under the local theory. The Num-LSA and TS, however, paint a more complex picture. Recall that, in these two cases, the electric capillary number Ca E is defined based on k p = 2πL −1 p , where L p is the size of the periodic domain and sets the largest possible length scale for the unstable eigenmodes. The dominant eigenmodes of shape and charge obtained by Num-LSA for a representative case are plotted in figure 7(a,b). Similar to the case of § 6.2 with flow only, all the modes are non-sinusoidal and exhibit strong variations near the stagnation point. This is true especially of the eigenmodes of charge, which display shock-like structures at 931 A25-15 0 -π/4 -π/4 π/4 π/4 0.5 x = 0. These shocks result from the advection by the applied flow of surface charges of opposite sign on each side of the stagnation point. They are reminiscent of the nonlinear shapes of figure 4, but are even more strongly concentrated near the origin (note the different scales in figure 7a,b). Note that the charge conservation equation (2.5) does not account for surface diffusion of charge, which, if included, may regularize the profiles. These sharp gradients seen in the linear eigenmodes are yet further amplified in the nonlinear regime, as we show by performing TS with a condition given by the first unstable eigenmode with a small amplitude. The evolution of the shape and charge profiles is shown in figure 7(c,d) (also see supplementary movie): the interface deflection sharpens rapidly as charges accumulate on each side of the stagnation point, leading ultimately to failure of our numerical method. The emergence of shocks in the charge distribution has also been observed in related configurations, such as in liquid drops under applied electric fields (Lanauze, Walker & Khair 2015;Das & Saintillan 2017a,b). There, the quadrupolar Taylor flow generated by tangential electric stresses at the drop interface sweeps surface charges from the poles towards the equator, resulting in sharp gradients at that location.
To further understand the interaction between the background flow and the electric field, we study the behaviour of the system as a function of local strain rateÂ in figure 8, in a case where the interface is electrically unstable. As already discussed above, the LT predicts that the applied flow is always destabilizing, especially in the presence of charge convection. The behaviour is more complex according to the Num-LSA, showing that the background flow in fact has a stabilizing effect for 0 <Â <Â c (Â c ≈ 0.118), as the growth rate decreases from its value atÂ = 0. BeyondÂ ≥Â c , the background flow becomes destabilizing even in the presence of charge convection. However, the growth rate is always smaller when charge convection is included in the model. Interestingly, the growth rates predicted by LT and Num-LSA are in close agreement in the absence of charge convection. We discuss in § 6.4 the mechanism for these trends, which involves the subtle interplay of convection with Ohmic conduction.
The effect of the conductivity ratio R and permittivity ratio Q on the stability of the system is studied in figure 9. Although both LT and Num-LSA yield qualitatively similar trends with respect to R and Q, the maximum growth rates obtained by the two methods differ significantly, and predict opposite effects of charge convection as already observed in figure 8. Another significant difference is that according to LT the interface is only unstable above critical values of R, Q ∼ O(1), whereas it is unstable for all values of R and Q according to the numerical stability. The evident discrepancy between the results of the two methods is attributed to the local approximation made in LT, where linear terms in x were neglected in the charge conservation equation. As demonstrated by Num-LSA, these coupling terms with non-constant coefficients result in very efficient charge transport towards the stagnation point, leading to strongly localized eigenmodes unlike the Fourier modes assumed by LT.

Mechanisms of charge transport in the dominant mode of instability
In order to explain the non-monotonic role of charge convection seen in figure 9, we further analyse the various mechanisms of charge transport in the dominant mode of instability. We define the Ohmic and convective fluxes aṡ over most of the domain. According to figure 10(a), charge convection is dominant in the vicinity of the stagnation point, whereas conduction takes over further away from the origin. Close to the stagnation point,q ohm exhibits oscillations, which are stronger with increasingÂ as shown in figure 10(b) but are suppressed when charge convection is neglected.
To elucidate the underlying mechanisms for this behaviour, we analyse the respective roles of interface deflections and charge perturbations in driving Ohmic currents in the dominant eigenmode. Recall that the eigenmodes obtained by Num-LSA involve perturbations in both shapeξ and chargeq. Here we estimate the Ohmic current induced by these eigenmodes in the linear regime: (6.11) To expressÊ n 1 andÊ n 2 in terms of (q,ξ), we linearize the boundary integral equation (5.2) to find EliminatingÊ n 1 andÊ n 2 using (6.12) and (6.13) yields the following expression for the charge conduction flux:q 14) which captures the conduction response of the system to small perturbations. The first term on the right-hand side represents the Ohmic flux induced by perturbing the shape while the second term is the flux induced by the charge perturbation. Figure 11 shows the profile of the perturbation in the dominant eigenmode along with the resulting Ohmic fluxes for the same set of parameters as used in figure 10. It is evident from figure 11(b) that the Ohmic currents induced by the applied deformation and charge distribution have opposite signs over the entire domain, and thus work against each other. This explains the oscillations inq ohm observed near the stagnation point in figure 10(a), and it is this complex Ohmic response that opposes charge convection in the dominant eigenmode, leading to the stabilizing effect of convection seen in figures 8 and 9. One should note that this behaviour is independent of the sign of (1 − RQ) since there is no preferred order to the arrangement of the fluid layers in the linear regime. In other words, two systems characterized by RQ and (RQ) −1 are dynamically equivalent.

Conclusions
We have presented a theoretical and numerical model in two dimensions to study the dynamics of an interface separating two immiscible fluid layers subject to a tangential electric field and a stagnation point flow. We performed a local linear stability analysis in the vicinity of the stagnation point, which was able to recover the previous results of Melcher & Schwarz (1968) in the absence of the background flow, and of Tseng & Prosperetti (2015) in the absence of the electric field and in the limit of zero inertia. Our local theory was complemented by a numerical analysis using the boundary element method, which was also used to perform a numerical normal-mode linear stability analysis based on the complete Melcher-Taylor leaky dielectric model including charge convection. Our results show that charge convection plays a significant role in determining the dynamics of the system by altering the dominant unstable mode, in which it was shown to have a stabilizing effect. Further, we explored the dynamics of the system far from equilibrium using transient nonlinear numerical simulations and demonstrated how the coupling of flow and interfacial charge dynamics in the dominant unstable mode gives rise to strongly nonlinear effects such as the formation of high-curvature tips and of charge density shocks. In this study, the convergent flow and the electric field were assumed to be independently applied. This differs from the case of the equatorial EHD instability in drops (Brosseau & Vlahovska 2017;Wagoner et al. 2020), where the flow is also generated by the electric field. In spite of this simplification, our analysis provides valuable insights into the underlying mechanisms responsible for the EHD equatorial streaming such as the evolution of the convergent line instability and the emergence of strong charge gradients. A more detailed discussion of the relevance of the present work to describe equatorial instabilities in liquid drops is provided in Appendix A. Extensions of the present study could include considering the effect of the viscosity contrast (λ / = 1) and of equilibrium surface curvature on the behaviour of the system. Further attempts to improve the accuracy of the numerical simulations may involve implementation of shock-capturing schemes for the solution of the charge conservation equation, as well as surface reparametrization schemes for handling extreme local deformations.
Appendix A. Relevance to the equatorial streaming instability Brosseau & Vlahovska (2017) reported streaming from the equator of a drop placed in a uniform electric field. In the experimental system, RQ > 1 and the EHD flow driven by electric shear stresses on the drop interface converges at the equator (see figure 1). At the equator the applied electric field is also parallel to the drop interface. Thus, the configuration resembles the set-up considered in our theoretical study. Here, we use the LT and Num-LSA analyses to gain insight into the interface destabilization at the stagnation line of the EHD flow.
The LT of a fluid interface subject to a convergent flow predicts that the interface is always unstable (Tseng & Prosperetti 2015). However, the LT for a fluid interface subjected to a tangentially applied DC uniform field is stable for the experimental conditions according to (6.1). In the experiments the streaming was only observed at sufficiently high electric fields, Ca E ∼ O(1), which is likely due to the competition of the destabilizing and stabilizing actions of the field-driven flow and of the electric field. Figure 12 shows the theoretical prediction for the growth rate of the instability. The strain rate of the convergent flow can be estimated from the EHD flow at the equator of the drop based on Taylor's classic solution (Taylor 1966): where fluids 1 and 2 represent the drop and the suspending liquid, respectively. We obtain A ≈ 2.22 s −1 using the experimental parameters for a drop of silicone oil (ρ 1 = 960 kg m −3 , 1 / 0 = 2.8, σ 1 = 1.4 × 10 −12 S m, μ 1 = 0.048 Pa s) suspended in castor oil (ρ 2 = 961 kg m −3 , 2 / 0 = 4.6, σ 2 = 4.4 × 10 −11 S m, μ 2 = 0.69 Pa s) under an electric field of E 0 = 4 kV cm −1 and a surface tension γ = 4.5 mN m −1 . The LT does predict a threshold Ca E , which increases with viscosity ratio. According to the Num-LSA, which can be only performed for λ = 1, the interface is always unstable. However, the growth rate is vanishingly small at low Ca E and the instability may not develop on the time scale of the experiment, which is of the order of 1 s. Indeed, only above Ca E ∼ 1 does the instability grow at rate faster than 1 s −1 . Increasing the viscosity ratio further slows down the growth of the instability, and suggests that streaming would not be observed for viscosity ratios greater than one, in qualitative agreement with the experiment. Finally, we note that while the present study focuses on planar interfaces, the equatorial streaming instability in drops occurs on a curved surface and the effects of base-state interfacial curvature on the instability remain unknown.