1. Introduction
An emulsion is a dispersion of drops of one fluid in another. Emulsions are ubiquitous in nature and are extensively used in chemical, coating, food, cosmetics, agricultural and medical industries to achieve products with desired properties. The shelf life of products such as paints, lotions, insecticides, medicines, etc., directly depends on the stability of the emulsions. Usually, emulsions are stabilised using compatible surfactants, also called emulsifiers. While the stability of emulsions is desired in most applications, in some applications, naturally occurring emulsions need to be phase-separated, for example in oil and petroleum industries. Interesting applications utilising the convection currents generated in droplets due to an electric field are emerging. For example, Penkova et al. (Reference Penkova, Pan, Hodjaoglu and Vekilov2006) show that the electric-field-driven flow in a droplet affects the nucleation of protein clusters. Chung & Oliver (Reference Chung and Oliver1990) showed that electric-field-induced circulation results in enhanced heat transfer in a spherical droplet translating in a dielectric medium. Electric field also affects the wetting characteristics which can be harnessed for various applications in microfluids (see Romero Herreros Reference Romero Herreros2014). Electric field can also be employed to control the size of droplets and bubbles from a nozzle as shown in Notz & Basaran (Reference Notz and Basaran1999) and Sunder & Tomar (Reference Sunder and Tomar2013). Application of an electric field across an emulsion has been shown to affect the stability of emulsions (Barnes Reference Barnes1994; Kilpatrick Reference Kilpatrick2012; Goodarzi & Zendehboudi Reference Goodarzi and Zendehboudi2019). For instance, in the oil refining industry, where the conventional methods for oil dehydration, namely centrifugation, chemical treatment and gravity separation, have shown only limited success, electric field has been recently used effectively for separation of water from emulsified oil (Mhatre & Thaokar Reference Mhatre and Thaokar2015). The effect of electric field on a liquid drop has been studied extensively since the time of Rayleigh (Reference Rayleigh1882). Garton & Krasucki (Reference Garton and Krasucki1964) observed in their experiments that bubbles in insulating liquids deform in prolate shapes and showed that with an increase in electric field beyond a critical value, the bubbles can become unstable. Curious flow patterns around a droplet were observed in the experimental investigations of Allan & Mason (Reference Allan and Mason1962), Bungenberg de Jong & Hoskam (Reference Bungenberg de Jong and Hoskam1941) and O'Konski & Harris (Reference O'Konski and Harris1957), and were explained in the seminal work of Taylor (Reference Taylor1966) by incorporating a small but finite conductivity for both the drop and the ambient fluid. Taylor (Reference Taylor1966) not only could explain the origin of the circulatory flow, but also showed that the flow generated by the tangential electric stresses at the drop surface can lead to oblate deformation of the droplets, where only prolate drop deformations are expected for perfectly conducting or perfectly dielectric fluids. The leaky dielectric theory, based on these results, is presented in detail in Melcher & Taylor (Reference Melcher and Taylor1969) and has been extensively used since then to explain the behaviour of an isolated droplet in an emulsion. Saville (Reference Saville1997) provides a detailed review of electrohydrodynamics showing the connection between the electrokinetic and the leaky dielectric theories (also see Zholkovskij, Masliyah & Czarnecki Reference Zholkovskij, Masliyah and Czarnecki2002).
Subsequent to the deformation predictions in the low-Reynolds-number limit by Taylor (Reference Taylor1966), several improvements on the solution have been proposed. Second-order corrections to drop deformation were given by Ajayi (Reference Ajayi1978), followed by the effect of charge convection on drop deformation as discussed in Feng & Scott (Reference Feng and Scott1996), Shutov (Reference Shutov2002), Xu & Homsy (Reference Xu and Homsy2006), Vlahovska (Reference Vlahovska2011), Fernandez (Reference Fernandez2013) and Lanauze, Walker & Khair (Reference Lanauze, Walker and Khair2015). The primary objective of those investigations was to resolve the discrepancy between some of the experimental results of Torza, Cox & Mason (Reference Torza, Cox and Mason1971) and the theoretical predictions obtained using Taylor's linear theory (Taylor Reference Taylor1966). Taylor's theory is valid only in the limit of small electric capillary number, defined as the ratio of the electric stresses and the capillary pressure, 
$Ca_E=\varepsilon E_\infty ^2 R/\gamma$, where 
$\varepsilon$ is the permittivity of the medium, 
$E_\infty$ is the strength of the externally applied electric field, 
$R$ is the radius of the droplet and 
$\gamma$ is the surface tension coefficient. Drop deformation at higher 
$Ca_E$ has also been investigated, mostly numerically (Sherwood Reference Sherwood1988; Feng & Scott Reference Feng and Scott1996; Feng Reference Feng1999; Basaran Reference Basaran2002; Lac & Homsy Reference Lac and Homsy2007; Supeene, Koch & Bhattacharjee Reference Supeene, Koch and Bhattacharjee2008; Fernandez Reference Fernandez2013; Wang, Wang & Qiu Reference Wang, Wang and Qiu2014; Das & Saintillan Reference Das and Saintillan2017a), but also theoretically by Ajayi (Reference Ajayi1978). Stability and breakup of perfectly dielectric droplets under large deformation have been investigated by Cheng & Chaddock (Reference Cheng and Chaddock1984) using energy arguments, and that of charged droplets in air by Brazier-Smith, Jennings & Latham (Reference Brazier-Smith, Jennings and Latham1971), and for leaky dielectric droplets in a leaky dielectric ambient fluid by Feng & Scott (Reference Feng and Scott1996) and Lac & Homsy (Reference Lac and Homsy2007). Giglio et al. (Reference Giglio, Gervais, Rangama, Manil, Huber, Duft, Müller, Leisner and Guet2008) studied the onset of Rayleigh instability in charged conducting microdroplets. The transient dynamics of droplets in electric field has been computationally studied by Sherwood (Reference Sherwood1988), Esmaeeli & Sharifi (Reference Esmaeeli and Sharifi2011), Lanauze, Walker & Khair (Reference Lanauze, Walker and Khair2013), Zhang, Zahn & Lin (Reference Zhang, Zahn and Lin2013) and Das & Saintillan (Reference Das and Saintillan2017b). The importance of the ratio of charge relaxation time scale to the time scale of other processes such as charge convection due to fluid flow or drop deformation due to local shear has been discussed in Saville (Reference Saville1971), Collins et al. (Reference Collins, Sambath, Harris and Basaran2013) and Lanauze et al. (Reference Lanauze, Walker and Khair2015). The effect of electric field on compound droplets in double emulsions has also been studied theoretically and experimentally (Gouz & Sadhal Reference Gouz and Sadhal1989; Tsukada et al. Reference Tsukada, Mayama, Sato and Hozawa1997; Ha & Yang Reference Ha and Yang1999; Behjatian & Esmaeeli Reference Behjatian and Esmaeeli2013; Soni, Juvekar & Naik Reference Soni, Juvekar and Naik2013; Abbasi et al. Reference Abbasi, Song, Kim and Lee2017). Most of these studies of compound droplets have focused on the relative deformation of the inner and the outer droplets, and also the possibilities of breakup and release of the inner droplet into the ambient fluid. Gouz & Sadhal (Reference Gouz and Sadhal1989) used a bipolar coordinate system to examine various configurations of compound droplets under which a compound droplet is stable in a translational flow for a given set of suitable electrical conductivities and permittivities of the constituent fluids. Ha & Yang (Reference Ha and Yang1999) used a domain perturbation method, similar to Taylor (Reference Taylor1966), to predict the equilibrium shapes of the inner and the outer droplets and also computed the rheological response of a dilute double emulsion of compound droplets. Recently, Abbasi et al. (Reference Abbasi, Song, Kim and Lee2017), using a level set method, studied the various modes under which a double emulsion may break up. Abbasi et al. (Reference Abbasi, Song, Kim and Lee2017) showed that the inner droplet may undergo large oblate deformation whereas the outer droplet undergoes prolate deformation, and thus can lead to bursting of the compound droplet ejecting the inner droplet fluid into the ambient. They showed that, if the inner droplet is not concentrically placed, it may migrate to the outer surface of the compound droplet and be ejected out for certain conductivity and permittivity ratios.
The rheological response of an emulsion has been shown to be altered in the presence of an electric field. Where a considerable effort has been invested in studying the behaviour of isolated droplets in electric fields, which has led to a better understanding of the rheology of dilute emulsions in electric fields (Ha & Yang Reference Ha and Yang1999; Vlahovska Reference Vlahovska2011; Sengupta, Walker & Khair Reference Sengupta, Walker and Khair2017; Mandal et al. Reference Mandal, Sinha, Bandopadhyay and Chakraborty2018), the effect of interactions between the droplets on the rheology and stability of dense emulsions needs to be understood better. There have been some studies of the interaction between droplets and particles suspended in a perfect insulator in the presence of an externally applied electric field (Latham & Roxburgh Reference Latham and Roxburgh1966; Brazier-Smith et al. Reference Brazier-Smith, Jennings and Latham1971; Arp, Foister & Mason Reference Arp, Foister and Mason1980). Due to the relevance to the oil industry and atmospheric science, interactions between conducting droplets suspended in an insulating fluid (such as in water-in-oil systems and water droplets in clouds) in the presence of an electric field have been studied extensively (Pearce Reference Pearce1954; Brazier-Smith Reference Brazier-Smith1971; Atten Reference Atten1993; Zhang, Basaran & Wham Reference Zhang, Basaran and Wham1995; Mhatre & Thaokar Reference Mhatre and Thaokar2015) and more recently by Sorgentone et al. (Reference Sorgentone, Kach, Khair, Walker and Vlahovska2020). A curious phenomenon of electrocoalescence, or rather the phenomenon of non-coalescence, of charged droplets in an electric field has also revived some interest lately in the study of interactions between a pair of droplets in the presence of an electric field (Aryafar & Kavehpour Reference Aryafar and Kavehpour2007; Ristenpart et al. Reference Ristenpart, Bird, Belmonte, Dollar and Stone2009; Anand, Juvekar & Thaokar Reference Anand, Juvekar and Thaokar2019; Sunder & Tomar Reference Sunder and Tomar2020). For uncharged droplets in perfectly dielectric fluids, the droplets always attract each other due to dipole–dipole electrostatic interactions. Interactions between a pair of leaky dielectric droplets suspended in a leaky dielectric fluid were first studied theoretically by Sozou (Reference Sozou1975). Sozou (Reference Sozou1975) used bispherical harmonics, considering the two drops to be of the same size, to solve for the electric field and the velocity field due to the tangential electric stresses at the drop interface. The key conclusion was that the non-uniformity in the electric field in the region between the droplets increases as the droplets approach distances shorter than two diameters of the droplets and affects the nature of deformation and coalescence dynamics of the droplets. Zhang et al. (Reference Zhang, Basaran and Wham1995) used a population dynamics approach to study the interaction between droplet pairs, but their study was conducted using perfectly conducting drops dispersed in a perfect dielectric liquid. Using a boundary integral method, Baygents, Rivette & Stone (Reference Baygents, Rivette and Stone1998) performed axisymmetric simulations of interactions between two leaky dielectric droplets in leaky dielectric ambient fluid with the axis of symmetry aligned with the electric field. They showed that the hydrodynamic forces dominate the electrostatic forces in the limit of small Reynolds number and dictate the stability of emulsions. An emulsion is rendered unstable when the droplets attract, whereas the emulsion is stable when the droplets repel. Droplets attract each other if the flow due to the electric stresses at the droplet surface is from the poles towards the equator, thus resulting in a velocity field which pulls the droplets together along the axis, with the externally applied electric field direction aligned with the axis. On the other hand, droplets repel each other when the flow is from the equator to the poles. The flow direction can be simply predicted by comparing the drop to ambient electric conductivity and permittivity ratios, 
$\sigma _r$ and 
$\varepsilon _r$, respectively. If 
$\sigma _r<\varepsilon _r$, flow is from poles to equator and thus droplets are expected to attract, whereas when 
$\sigma _r>\varepsilon _r$ flow is from the equator to the poles. Baygents et al. (Reference Baygents, Rivette and Stone1998) showed that the relative velocity between the droplets of perfectly dielectric systems is governed by 
$h^{-4}$, where 
$h$ is the separation distance between the droplets. They showed that, in contrast, for leaky dielectric fluids, droplets attract or repel at a rate proportional to 
$h^{-2}$ due to the hydrodynamic forces. Tomar et al. (Reference Tomar, Gerlach, Biswas, Alleborn, Sharma, Durst, Welch and Delgado2007) proposed a CLSVOF-based method to simulate two-phase leaky dielectric electrohydrodynamics and showed good agreement with the results for drop–drop interaction of Baygents et al. (Reference Baygents, Rivette and Stone1998). Lin, Skjetne & Carlson (Reference Lin, Skjetne and Carlson2012) used a phase-field method to simulate interaction between a pair of droplets in the presence of an electric field. The effect of viscosity ratio on the time scale of coalescence was discussed, but the scaling for the relative velocity between the droplets was not examined.
In the present study, we investigate drop–drop interactions in single and double emulsions in the presence of an electric field using a two-phase flow solver based on a volume of fluid method. The scaling for the relative velocity between the interacting droplets is studied first and it is shown that it can vary as 
$h^{-2}$ or 
$h^{-4}$ (where 
$h$ is the centre-to-centre distance between the droplets) depending upon the Reynolds number of the flow generated by the tangential electric stresses. Then, we present simulations for low-viscosity fluids (with the ambient to drop viscosity ratio as one), for which the droplets attract initially but eventually repel each other. For certain electric properties of the drops and the ambient fluid, we show that the droplets attract, but the coalescence phenomenon is delayed due to the slow drainage of the film between the droplets and a drop-pair doublet is formed similar to the particle pairs formed in Arp & Mason (Reference Arp and Mason1977), chain of droplets in Holto, Berg & Lundgaard (Reference Holto, Berg and Lundgaard2009) and cell–cell (vesicle) doublets in Zimmermann & Vienken (Reference Zimmermann and Vienken1982). These nonlinear effects due to the finite Reynolds number have not been explored before. Finally, we study the behaviour of compound droplets in an electric field and the effect of electric field on the drop–drop interactions between compound droplets.
The paper is organised as follows. The problem definition and numerical scheme used in this study are discussed in § 2. Results and discussions of the simulations of single and compound droplets in single and double emulsions, respectively, are presented in § 3. Important conclusions from the study are discussed in § 4.
2. Problem formulation
We study the interaction dynamics of a pair of leaky dielectric single and compound droplets in a leaky dielectric ambient fluid under the influence of an externally applied electric field using an axisymmetric formulation. The initial configuration of the emulsion droplets for the computations is shown in the schematic in figure 1. The dashed line marks the axis of symmetry (along the horizontal direction). The extent of the domain in the radial and axial directions is 
$5 R_o$ and 
$10 R_o$, respectively, where 
$R_o$ is the radius of the outer shell of the compound droplet (see figure 1). The chosen computational domain is large enough to neglect the effect of the boundary on the dynamics of the droplets at the time scales considered in this study. The initial separation distance between the droplets is 
$4 R_o$. An external electric field is applied along the axial direction and is imposed using 
$\phi = \phi _0$ as the boundary condition on the right-hand boundary, where 
$\phi$ denotes the electric potential (marked in figure 1). The left-hand boundary is electrically grounded using the boundary condition 
$\phi = 0$. Neumann boundary condition for the electric potential is imposed on the top boundary (
$r = 5R_o$), 
$\partial \phi /\partial r = 0$, where 
$\partial /\partial r$ denotes the derivative in the radial direction. Symmetry boundary conditions are imposed on the axis of symmetry for the velocity field and the electric potential, and no-slip velocity boundary conditions are imposed on the other boundaries.
Schematic (not to scale) showing a pair of compound droplets interacting in the presence of an externally applied electric field. The radii of the outer and inner droplets are 
$R_o$ and 
$R_i$, respectively. Various electric, transport and thermodynamic properties are marked in the schematic.
We consider the droplets to be neutrally buoyant with densities 
$\rho _1 = \rho _2$. Here, the properties of the outer fluid are denoted by subscript 
$1$, whereas the fluid constituting the shell of the compound droplet is denoted using the subscript 
$2$ (as marked in figure 1). The properties of the fluid constituting the dispersed phase of the compound droplet are considered to be the same as those of the outer fluid and are denoted by subscript 
$1$. Viscosities of the fluids 
$1$ and 
$2$ are given by 
$\mu _1$ and 
$\mu _2$, respectively. Electric permittivities and conductivities are 
$\varepsilon _1$ and 
$\sigma _1$, respectively, for the ambient fluid and 
$\varepsilon _2$ and 
$\sigma _2$ for the outer-drop fluid. The two-phase system can be categorised by the following ratios of different properties: 
$\rho _r = \rho _2/\rho _1$, 
$\mu _r = \mu _2/\mu _1$, 
$\varepsilon _r = \varepsilon _2/\varepsilon _1$ and 
$\sigma _r = \sigma _2/\sigma _1$. For the simulations presented in this study, we assume 
$\rho _r = 1$ and 
$\mu _r = 1$. The fluids involved are assumed to be incompressible and Newtonian, and the surface tension coefficient of the interface between the two fluids is given by 
$\gamma _{12} = \gamma$. In the following, we discuss the governing equations and the numerical formulation employed for the simulations.
The two-phase system discussed above is governed by the modified incompressible Navier–Stokes equations using the one-fluid formulation as
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}
\begin{gather}\rho (C) \left( \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{uu} \right) ={-}\boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot} [\mu (C) ( \boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{u}^T ) ] + \gamma \kappa \boldsymbol{n} \delta_s + \boldsymbol{F}_{el}, \end{gather}
where 
$\boldsymbol {u}$ is the velocity field, 
$p$ is the pressure and the indicator function 
$C$ is given by the volume fraction that takes a value of zero in fluid 1 and one in fluid 2. The density and viscosity are functions of the volume fraction and are given by 
$\rho ( C ) = \rho _2 C + \rho _1 (1-C)$ and 
$\mu ( C ) = \mu _2 C + \mu _1 (1-C)$. Surface tension force is modelled as a volumetric force given by 
$\gamma \kappa \boldsymbol {n} \delta _s$, where 
$\gamma$ is the surface tension coefficient, 
$\kappa$ is the curvature at the interface, 
$\boldsymbol {n}$ is the normal vector and 
$\delta _s$ is the smoothed Dirac delta function. Similarly, electric field force 
$\boldsymbol {F}_{el}$ can be modelled as a volumetric force given by (Melcher & Taylor Reference Melcher and Taylor1969)
\begin{equation} \boldsymbol{F}_{el} ={-} \tfrac{1}{2} ( \boldsymbol{E} \boldsymbol{\cdot} \boldsymbol{E} ) \boldsymbol{\nabla} \varepsilon + q_v \boldsymbol{E}, \end{equation}
where 
$\boldsymbol {E}$ is the electric field and 
$q_v$ is the volumetric free charge. The first term in the above equation is the dielectrophoretic force and the second term is the electrophoretic force arising due to the presence of free charges. The electric force can also be written in terms of the divergence of the Maxwell stress tensor as follows:
\begin{equation} \boldsymbol{F}_{el} = \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_e, \end{equation}
where the Maxwell stress tensor 
${\boldsymbol{\mathsf{T}}}_e$ is given by
\begin{equation} {\boldsymbol{\mathsf{T}}}_e = \varepsilon \left( \boldsymbol{EE} - \frac{E^2}{2}\boldsymbol{I} \right). \end{equation}
Since the time scales at which the magnetic induction effects become important are significantly smaller than the flow time scales in this study, we ignore the electromagnetic coupling and assume that the electric field 
$\boldsymbol {E}$ is curl-free 
$( \boldsymbol {\nabla } \times \boldsymbol {E} = 0 )$. Thus, we can write 
$\boldsymbol {E} = -\boldsymbol {\nabla } \phi$. We solve the following Poisson equation with the boundary conditions discussed earlier:
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \varepsilon \boldsymbol{E} = \boldsymbol{\nabla} \boldsymbol{\cdot} (\varepsilon (- \boldsymbol{\nabla} \phi) ) = q_v. \end{equation}The charge conservation equation is given by
\begin{equation} \frac{\partial q_v}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} q_v + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{J} = 0. \end{equation}
Here, 
$\boldsymbol {J}= ( \sigma \boldsymbol {E} )$ is the current density. The electric permittivity and conductivity, in the one-fluid formulation, are given by 
$\varepsilon = \varepsilon _2 C + \varepsilon _1 (1-C)$ and 
$\sigma = \sigma _2 C + \sigma _1 (1-C)$, respectively. The evolution of the interface is governed by an advection equation in terms of the volume fraction (
$C$):
\begin{equation} \frac{\partial \textit{C}}{\partial t}+ \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} C = 0. \end{equation}
We note that when the electric conductivities of the two fluids are such that the charge relaxation time constants, 
$\varepsilon _1/\sigma _1$ and 
$\varepsilon _2/\sigma _2$, are small in comparison to the flow time scales, the charge conservation equation reduces to 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {J} = 0$. Thus, the bulk is essentially charge-free and the free charge accumulation only occurs at the interface where there is change in the relaxation time constant. Therefore, both the dielectrophoretic forces and the electrophoretic forces are surface forces acting only at the interface. The surface tension force and the surface electric force are modelled as volumetric forces in the one-fluid formulation using a modified version of the continuum surface force framework proposed by Brackbill, Kothe & Zemach (Reference Brackbill, Kothe and Zemach1992).
The above one-fluid formulation can be shown to impose the following boundary conditions implicitly. Using Gauss’s law at the interface yields the jump condition in the normal component of the electric field:
\begin{equation} \|\varepsilon \boldsymbol{E} \boldsymbol{\cdot} \boldsymbol{n}\|=q_s, \end{equation}
where 
$\|\cdot \|$ indicates the jump in the quantity across the interface in the direction of the normal vector 
$\boldsymbol {n}$. Here, 
$q_s$ is the surface charge density of free charge and 
$\boldsymbol {n}$ is the normal at the interface. Since no phase change is involved, we have velocity continuity at the interface:
\begin{equation} \|\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}\|=0. \end{equation}Normal stress balance condition at the interface yields
\begin{equation} \|p\|=\gamma\kappa + \|\boldsymbol{n} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_v \boldsymbol{n}\| + \|\boldsymbol{n} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_e \boldsymbol{n}\|, \end{equation}
where 
$\|p\|$ denotes the jump in pressure, 
$\gamma \kappa$ is the component of the pressure jump due to surface tension, 
$\|\boldsymbol {n} \boldsymbol {\cdot } {\boldsymbol{\mathsf{T}}}_v \boldsymbol {n}\|$ is the jump in the normal component of the viscous stress 
$({\boldsymbol{\mathsf{T}}}_v=\mu (C)(\boldsymbol {\nabla } \boldsymbol {u}+\boldsymbol {\nabla } \boldsymbol {u}^T))$ across the drop surface and 
$\|\boldsymbol {n} \boldsymbol {\cdot } {\boldsymbol{\mathsf{T}}}_e \boldsymbol {n}\|$ is the jump in the normal component of the Maxwell stress tensor defined in (2.5). In the absence of any variation in the surface tension coefficient, tangential stress balance is given by
\begin{equation} \|\boldsymbol{t} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_v \boldsymbol{n}\| + \|\boldsymbol{t} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_e \boldsymbol{n}\|. \end{equation}Conservation of the surface charge is ensured by accounting for advection of the interface, along with the change in area of the interface due to motion. Thus, surface charge conservation equation is given by
\begin{equation} \frac{\partial q_s}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}_s q_s = q_s \boldsymbol{n} \boldsymbol{\cdot} (\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\nabla} )\boldsymbol{u}-\|\sigma\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{E}\|+\|q_v \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{n}\|, \end{equation}
where the term 
$q_s \boldsymbol {n} \boldsymbol {\cdot } (\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {\nabla })\boldsymbol {u}$ accounts for the change in charge density due to the change in the interfacial area. Variation in surface charge density due to the difference in the current density across the interface is given by the term 
$\|\sigma \boldsymbol {n}\boldsymbol {\cdot } \boldsymbol {E}\|$. The term 
$\|q_v \boldsymbol {v} \boldsymbol {\cdot } \boldsymbol {n}\|$ corresponds to the transport of charge from the bulk to the interface by the fluid motion near the interface as suggested by (22
$^{\prime }$) in Saville (Reference Saville1997).
In the one-fluid formulation, we solve a volumetric charge conservation equation which also accounts for the surface charge conservation in a thin diffused layer (over one to three grid cells) at the interface. The efficacy of the one-fluid formulation has been demonstrated in Tomar et al. (Reference Tomar, Gerlach, Biswas, Alleborn, Sharma, Durst, Welch and Delgado2007) and López-Herrera, Popinet & Herrada (Reference López-Herrera, Popinet and Herrada2011). The algorithm used in the present case is from the latter study which employs a charge conservation approach.
The governing equations for mass, momentum and charge conservation along with the Poisson equation for the electric field potential are solved numerically in the one-fluid formulation. Geometric volume of fluid, coupled with level set method, has been used previously for several electrohydrodynamic problems (see e.g. Tomar et al. Reference Tomar, Gerlach, Biswas, Alleborn, Sharma, Durst, Welch and Delgado2007; Sunder & Tomar Reference Sunder and Tomar2020). Here, in order to investigate the coalescence outcomes of compound droplets in electric field, we employ an open source, adaptive mesh refinement (AMR)-based, parallel-flow solver: Basilisk (Popinet Reference Popinet2015). Basilisk, similar to Gerris (Popinet Reference Popinet2003, Reference Popinet2009; Tomar et al. Reference Tomar, Fuster, Zaleski and Popinet2010), is based on an Oct-Tree AMR. López-Herrera et al. (Reference López-Herrera, Popinet and Herrada2011) simulated electrohydrodynamics by modifying Gerris and compared the accuracy of the numerical scheme with that of Tomar et al. (Reference Tomar, Gerlach, Biswas, Alleborn, Sharma, Durst, Welch and Delgado2007). A similar numerical scheme has been used in Basilisk as well. A geometric volume of fluid method is implemented in Basilisk with the interface approximated by a piecewise linear reconstruction. In Basilisk, normal and curvature at the interface are computed using the height function approach, similar to that in Popinet (Reference Popinet2009). Geometric advection of volume of fluid (equation (2.8)) is carried out using a second-order Godunov scheme (Bell, Colella & Glaz Reference Bell, Colella and Glaz1989) and a balanced-force discretisation, for both electric stresses and surface tension force, is used to minimise spurious currents (Francois et al. Reference Francois, Cummins, Dendy, Kothe, Sicilian and Williams2006). A collocated grid discretisation is used and an approximate projection scheme is used to solve the Navier–Stokes equations. Details of the discretisation of the governing equations ((2.2), (2.6), (2.7) and (2.8)) along with detailed validations can be found in Popinet (Reference Popinet2015) and López-Herrera et al. (Reference López-Herrera, Popinet and Herrada2011).
We use AMR to enhance the computational efficiency while maintaining the required accuracy of the simulations. In the present study, we employ AMR using a cost function based on the gradient of the volume fraction, thus using a very fine mesh near the interface while maintaining a relatively coarser mesh elsewhere (see Popinet (Reference Popinet2015) for details of the AMR method). To understand the grid requirements for the simulations involving compound droplets presented in this study, we perform a grid independence study. We use the Taylor deformation parameter, defined as 
$D = (L-B)/(L+B)$, where 
$L$ is the major axis and 
$B$ is the minor axis of the deformed compound droplet, to compare numerical convergence with the increase in refinement. We perform simulations using different grid sizes and check the temporal evolution of 
$D$ as shown in figure 2. The inset shows a zoomed view for minimum value of deformation. The curves corresponding to 
$R_o/\varDelta x_{min} = 102$ and 
$205$ are very close with a maximum difference of 
$0.6\,\%$. Here, 
$\varDelta x_{min}$ is the minimum grid size used in the present study. For the simulations presented in this study, the entire interface is resolved using 
$\varDelta x_{min}$ and the grid is progressively coarsened away from the droplet. Based on the grid independence study, we use 
$R_o/\varDelta x_{min} = 102$ for the simulations presented in this study.
Temporal variation of deformation of the shell under a steady electric field, using five different grid refinements. The parameters considered for the simulations are electric capillary number 
$Ca_E = 0.2$, 
$Oh = 1.41$, 
$\sigma _r = 0.1$, 
$\varepsilon _r = 2$ and 
$Re_E=0.772$.
The following non-dimensional parameters govern the dynamics of droplets due to electrostatic forces and the flow generated by them. Electric capillary number 
$Ca_E = \mu _1 U_e/\gamma$, where 
$\mu _1$ is the viscosity of the ambient continuous phase, 
$\gamma =\gamma _{12}$ is the surface tension coefficient at the interface of fluids 1 and 2 and 
$U_e$ is the velocity scale defined based on the flow generated by the electric forces. The scaling for velocity generated due to the electric forces can be obtained by balancing viscous stresses with Maxwell stresses at the interface. The scale for the velocity is thus given by 
$U_e \sim E_\infty ^2 \varepsilon _1 R_o/\mu _1$. We note here that the above scaling with an additional factor of 
$(1 - \tau _e^2/\tau _e^1)\sigma _r/(\sigma _r + 2)^2$ would yield a more appropriate scaling following the solution for velocity given by (27) in Taylor (Reference Taylor1966), where 
$\tau _e^1 = \varepsilon _1/\sigma _1$ and 
$\tau _e^2 = \varepsilon _2/\sigma _2$ are the charge relaxation time constants for fluid 1 (ambient) and fluid 2 (drop), respectively. Indeed, in our simulations for 
$\tau _e^1 = \tau _e^2$, no electrohydrodynamic flow is observed (not shown here). Thus, the above scaling for velocity, in some cases, may be an order-of-magnitude different from the observed values. However, for non-dimensionalising velocity, we use 
$U_e$ defined above as also proposed in Baygents et al. (Reference Baygents, Rivette and Stone1998), since it allows the delineation of the effects of electrical parameters while comparing results for different values of 
$\sigma _r$ and 
$\varepsilon _r$. Nevertheless, later when we compute scaling for electrostatic and hydrodynamic forces, as well as flow, interface and charge relaxation time scales, we account for the effect of other electrical properties of the fluids on velocity scaling. Using 
$U_e$, the electric capillary number can now be written as 
$Ca_E = E_\infty ^2 \varepsilon _1 R_o/\gamma$.
The ratio of the viscous capillary time scale, 
$\tau _\mu ^\gamma = \mu _1 R_o/\gamma$, and the Rayleigh inertial capillary time scale, 
$\tau _\rho ^\gamma = \sqrt {\rho _1 R_o^3/\gamma }$, yields a measure of the role of viscosity in the interfacial dynamics, called the Ohnesorge number: 
$Oh_1 = \mu _1/\sqrt {\rho _1 R_o \gamma }$ for the ambient fluid; this can be defined similarly for the drop fluid: 
$Oh_2 = \mu _2/\sqrt {\rho _2 R_o \gamma }$. The Reynolds number, based on 
$U_e$, can be defined as 
$Re = \rho _1 R_o^2 E_\infty ^2 \varepsilon _1/\mu _1^2 = Ca_E/Oh_1^2$. An electrical equivalent of the hydrodynamic Reynolds number can also be defined by taking the ratio of charge convection to conduction: 
$Re_E = U_e \tau _e^1 (1 + \varepsilon _r)/(R_o (1 + \sigma _r))$ as suggested in Feng (Reference Feng1999). Feng (Reference Feng1999) showed that the electric Reynolds number, 
$Re_E$, defined above has an effect on the extent of deformation of the droplet.
Using the definition of the electric capillary number 
$Ca_E$, we can non-dimensionalise the electric field intensity with 
$\sim \sqrt {\gamma /\varepsilon _1 R_o}$ (as suggested by Feng (Reference Feng1999)), and thus the non-dimensionalised externally applied electric field intensity can now be written as 
$E_\infty ^*\sim \sqrt {Ca_E}$. Another important parameter requiring attention is the amount of free charge per unit area that develops at the interface: 
$Q_s \sim 3\varepsilon _1 (1 - \tau _2^e/\tau _1^e){E}_\infty /(\sigma _r + 2)$ (computed from an exact solution for a spherical droplet). The ratio of the charge relaxation time constants for the two fluids 
$(\tau _2^e/\tau _1^e)$ in the above expression essentially indicates the nature of the charge distribution on the surface of the drop. If the ratio is greater than one, the ambient fluid carries away the charge rapidly and a negative charge develops at the interface (in the northern hemisphere of the droplet with the axis of the droplet aligned with the applied electric field); whereas if the ratio is less than one, the drop fluid supplies positive charge in response to the applied electric field, at a more rapid rate than can be carried away by the ambient fluid, and thus a positive charge accumulates at the surface. The positive or negative charge determines the direction of the tangential stress and thus also the sense of circulation produced by the electrohydrodynamic forces, which in turn determines whether the nature of the interaction between the droplets is attractive or repelling (see Baygents et al. Reference Baygents, Rivette and Stone1998). At low Reynolds number, hydrodynamic viscous forces are stronger and drive the interaction between the droplets, whereas at higher Reynolds numbers, as shown in the next section, the relative velocity during attraction is governed by the electrostatic interaction between the droplets.
In the next section, we discuss the effects of an externally applied electric field on a pair of single and compound droplets.
3. Results and discussion
In this section, we study the interaction dynamics of a pair of droplets, single and compound, suspended in another immiscible fluid under the action of an externally applied electric field. We perform two-dimensional axisymmetric simulations with outer droplet radius 
$R_o$ and inner droplet radius 
$R_i$, as shown in the schematic in figure 1. First, we present results for droplets in a single emulsion (
$R_i = 0$). Based on the interaction dynamics of a pair of droplets, we comment on the stability of an emulsion in the presence of an electric field. Subsequently, we discuss interaction dynamics between a pair of compound droplets.
3.1. Single emulsion in an externally applied electric field
We validate the numerical model employed in this study by comparing the deformation of a droplet with the theoretical predictions of Taylor (Reference Taylor1966) and also previous computational studies (Tomar et al. Reference Tomar, Gerlach, Biswas, Alleborn, Sharma, Durst, Welch and Delgado2007; López-Herrera et al. Reference López-Herrera, Popinet and Herrada2011). The size of the domain is the same as that used previously (Tomar et al. Reference Tomar, Gerlach, Biswas, Alleborn, Sharma, Durst, Welch and Delgado2007; López-Herrera et al. Reference López-Herrera, Popinet and Herrada2011). All the computations have been performed up to the steady state. Figure 3 shows variation in 
$D$ with 
$\sigma _r$ for 
$\varepsilon _r = 10$, 
$Ca_E = 0.18$ and 
$Oh = 3.16$. We observe that the droplet deforms into an oblate shape (
$D < 0$) for low values of 
$\sigma _r$ and into a prolate shape (
$D > 0$) for higher values of 
$\sigma _r$. In the presence of an electric field, depending on the electrical properties of the two phases, the drop deforms either into an oblate or a prolate shape. Taylor (Reference Taylor1966) provided an analytical expression for the deformation of a drop as a function of the fluid properties and the electric field intensity, in the limit of small deformation, as
\begin{equation} D=\frac{9}{16}\frac{Ca_E}{( 2 + \sigma_r )^2} \left[ 1 + \sigma_r^2 - 2\varepsilon_r + \frac{3}{5} ( \sigma_r - \varepsilon_r ) \frac{2+3\mu_r}{1+\mu_r} \right], \end{equation}
where 
$D = (L-B)/(L+B)$ is the Taylor deformation parameter, with 
$L$ being the length of the droplet along the direction of the electric field and 
$B$ the width of the droplet in the plane perpendicular to the applied electric field.
Comparison of deformation 
$(D)$ of a single droplet with theoretical and computational studies. Other non-dimensional parameters for the computations are 
$Ca_E = 0.18$, 
$Oh = 3.16$ and 
$\varepsilon _r = 10$.
The accuracy of Taylor's linear theory has been discussed in several previous studies (Ajayi Reference Ajayi1978; Feng Reference Feng1999), especially in the context of not-so-good agreement with some of the experiments of Torza et al. (Reference Torza, Cox and Mason1971) for low-conductivity fluids. In particular, the deformation of the drop is expected to be proportional to 
$E_\infty ^2$ from the theory, but the experiments of Torza et al. (Reference Torza, Cox and Mason1971) showed only qualitative agreement and the extent of deformation in some cases was over twofold different from the predictions of the theory. Feng (Reference Feng1999) suggested that the possible origin of this discrepancy could be due to the neglect of charge convection in Taylor's theory (see (2.3) in Feng Reference Feng1999). At low, but finite, electric Reynolds numbers 
$(Re_E)$, the deviation was shown to be significant especially for very low-conductivity fluids such as silicone and castor oils. In our simulations, we account for charge convection effects. We observe that for small drop deformation there is good agreement between the computational results and Taylor's theory (as shown in figure 3), and only at larger 
$D$ is there deviation from the theory essentially due to the breakdown of the assumption of linearity. Comparison with Ajayi's second-order theory shows good agreement with the simulation results. We note that we have focused here on cases where electrical conductivities of the fluids investigated are not as low as those considered by Feng (Reference Feng1999). We will address the physics of the interaction between droplets in very poorly conducting emulsions (e.g. silicone oil droplets in castor oil) in a separate study.
Figure 4 shows the drop shapes and the streamlines of the flow generated by the tangential electric forces. Due to the accumulation of the free charge at the drop interface, a tangential electric force acts at the interface that leads to fluid flow inside the droplet as well as in the ambient fluid, as shown in figure 4 for three different values of 
$\sigma _r$. Other parameters for all three cases shown in figure 4 are 
$Ca_E = 0.18$, 
$Oh = Oh_1 = 3.16$ and 
$\varepsilon _r = 10$. For 
$\sigma _r = 1.81$ and 
$\sigma _r = 5.1$, the flow is from poles to equator, resulting in the oblate shape for 
$\sigma _r = 1.81$ but the droplet remains spherical for 
$\sigma _r = 5.1$. As indicated by the deformation curve shown in figure 3, a droplet can deform prolately even for 
$\sigma _r < \varepsilon _r$ (for 
$\sigma _r > 5.1$ when 
$\varepsilon _r = 10$ and 
$\mu _r = 1$) when the electrostatic forces normal to the interface are stronger compared with the hydrodynamics forces. For 
$\sigma _r = 5.1$, 
$\varepsilon _r = 10$ and 
$\mu _r = 1$, the droplet remains spherical. This no-deformation case corresponds to the electrical and fluid properties for which the Maxwell stresses are exactly balanced by the normal viscous stress and the surface tension forces for a spherically shaped droplet. Nevertheless, non-zero tangential electric stresses generate a flow from poles to equator as shown in figure 4(b). The switchover of the flow from poles-to-equator to equator-to-poles occurs when 
$\sigma _r > \varepsilon _r$ as shown in figure 4(c) for 
$\sigma _r = 14$ and 
$\varepsilon _r = 10$. All the computations shown here are in good agreement with the analytical results, as well as with the previous computational studies (Tomar et al. Reference Tomar, Gerlach, Biswas, Alleborn, Sharma, Durst, Welch and Delgado2007; López-Herrera et al. Reference López-Herrera, Popinet and Herrada2011).
Streamlines and variation in the drop shape for different conductivity ratios. The electric Reynolds number 
$Re_E$ is (a) 2.228, (b) 1.027 and (c) 0.417. Other parameters are 
$\varepsilon _r = 10$, 
$Ca_E = 0.18$ and 
$Oh = 3.16$.
We now investigate interaction dynamics of two droplets in a single emulsion in the presence of an electric field. The computational domain size for these simulations is 
$17 R_o \times 34 R_o$. Charge accumulation at the interface of initially neutrally charged droplets results in the formation of electric dipoles and the droplets interact electrostatically through dipole–dipole electrostatic forces. In the large-separation limit, the electrostatic force between the droplets can be shown to be proportional to 
$E_\infty ^2/h^4$, where 
$h$ is the centre-to-centre spacing between the droplets. On the other hand, flow generated by the electric stress would generate a drag proportional to the relative velocity in the low-
$Re$ limit. Thus, the hydrodynamic drag force on one droplet due to the flow generated at the other is proportional to 
$E_\infty ^2/h^2$. However, at higher 
$Re$, the hydrodynamic drag is proportional to velocity squared, and thus a drag force proportional to 
$E_\infty ^4/h^4$ should be expected. Simulations presented in this study are for low 
$Ca_E$ values that correspond to low 
$Re = Ca_E/Oh^2 \sim O(1)$ and therefore we expect the linear drag law to be valid.
As discussed earlier, the direction of the flow depends upon the sign of 
$(1 - \tau _e^2/\tau _e^1)$, and therefore is expected to also govern the hydrodynamic interactions between the droplets in an emulsion. Based on the convection currents, the droplets are expected to move towards each other for 
$\tau _e^2/ \tau _e^1 > 1$, or equivalently 
$\sigma _r < \varepsilon _r$, whereas they will move away from each other for 
$\tau _e^2/\tau _e^1 < 1$, that is, 
$\sigma _r > \varepsilon _r$. Here, we present two cases with parameters chosen from Baygents et al. (Reference Baygents, Rivette and Stone1998) such that the droplets attract in the first case with 
$\sigma _r = 6$ and 
$\varepsilon _r = 8$ (
$\tau _e^2/\tau _e^1 = \varepsilon _r/\sigma _r = 1.33 > 1$) as shown in figure 5 and therefore are expected to render the emulsion unstable, whereas the droplets repel each other for 
$\sigma _r = 1.04$ and 
$\varepsilon _r = 0.2$ (
$\tau _e^2/\tau _e^1 \sim 0.2 < 1$) as shown in figure 6, and the electric field is expected to stabilise the emulsion.
Drop interaction along with the streamline patterns observed at different non-dimensional time 
$(t/\tau _\rho ^\gamma )$ for 
$\sigma _r = 6$ and 
$\varepsilon _r = 8$ at 
$Ca_E = 1.5$ and 
$Oh = 0.63$. The electric Reynolds number is 
$Re_E = 4.86$.
Drop interaction along with the streamline patterns observed at different non-dimensional time 
$(t/\tau _\rho ^\gamma )$ for 
$\sigma _r = 1.04$ and 
$\varepsilon _r = 0.2$ at 
$Ca_E = 1.5$ and 
$Oh = 0.63$. The electric Reynolds number is 
$Re_E = 2.223$.
In the first case with 
$\tau _e^2/\tau _e^1 > 1$, charge develops much faster 
$(t/\tau _\rho ^\gamma \sim 0.02\text {--}0.1 \ll 1)$ than the morphological evolution of the droplet 
$t/\tau _\rho ^\gamma \sim 2$. Also, the morphological evolution is much faster than the drop migration time scale 
$t/\tau _\rho ^\gamma \sim 40$. Thus, the drops first deform and subsequently without deforming much, either migrate towards each other or move away depending upon the ratio 
$\tau _e^2/ \tau _e^1$. The streamlines for 
$t/\tau _\rho ^\gamma \sim 1$ are similar (see figures 5a and 6a) for both cases due to the initial prolate deformation. Subsequently, the electric stresses result in a flow which is opposite for the two cases as shown in figures 5(b) and 6(b). Droplets in figure 5 coalesce, whereas those shown in figure 6 move apart.
Charge distributions for the two cases are presented in figure 7. As discussed earlier, the signs of the charges developing on the drop surfaces are opposite for the two cases shown in figure 7(a,b). Free charge distribution on the interface for a spherically shaped droplet is given by 
$Q_s \cos (\theta )$, where the angle 
$\theta$ is measured from the horizontal direction aligned with the applied electric field (as marked in the insets of figure 8). Thus, the charge is maximum at the poles and zero at the equator. For the first case, negative charge develops on the northern hemisphere of the droplet with the axis of the droplet aligned with the electric field, whereas positive charge appears for the second case. Figure 8 shows the variation in the charge density on the drop surface for different time instances. Initially when the droplets are far apart, charge develops symmetrically on the northern and the southern hemispheres, similar to the isolated droplet case. However, when the interfaces of the droplets come closer than one diameter, due to non-uniformity in the electric field in the region between the droplets, charge increases in the regions where the droplets face each other.
Variation in volumetric charge density 
$(q_v)$ along the interface for (a) 
$\sigma _r = 6$ and 
$\varepsilon _r = 8$ (
$Re_E=4.86$) and (b) 
$\sigma _r = 1.04$ and 
$\varepsilon _r = 0.2$ (
$Re_E=2.223$) at 
$Ca_E = 1.5$ and 
$Oh = 0.63$.
Temporal variation in volumetric charge density 
$(q_v)$ along the interface and the charge distribution around the interface for (a) 
$\sigma _r = 6$ and 
$\varepsilon _r = 8$ (
$Re_E=4.86$) and (b) 
$\sigma _r = 1.04$ and 
$\varepsilon _r = 0.2$ (
$Re_E=2.223$) at 
$Ca_E = 1.5$ and 
$Oh = 0.63$.
The dynamics of the droplets shown in the sequence of images in figures 5 and 6 can be understood by comparing the time scales for the different phenomena active during the interaction between the droplets. The time scale of charge accumulation at the drop surface is given by 
$\tau _e^1/\tau _\rho ^\gamma$, which is much shorter due to relatively high conductivities of the fluids considered in the simulations presented in this study. The time scale for the development of the flow due to the tangential electric stresses is given by the time scale for the diffusion of the momentum 
$\tau _{\mu } \sim h^2 \rho _1/\mu _1$. Thus, 
$\tau _{\mu }/\tau _\rho ^\gamma \sim h^2/(Oh\,R_o^2)$ yields the time scale for the influence of the electric-stress-driven flow near one drop to reach the other drop. If the droplets are far apart, it will take longer for the flow to develop and thus electrostatic forces may initially influence the droplets before being overtaken by the hydrodynamic forces. Similarly, for low 
$Oh$ values, even for droplets as close as 
$h/R_o \sim 5$, 
$\tau _{\mu }/\tau _\rho ^\gamma$ will be large and the droplets would first respond to the attractive electrostatic forces. Figure 9 demonstrates the above phenomenon for low 
$Oh = 0.02$ and high 
$Re_E=44.131$. The electric Reynolds number is much higher relative to previous cases, thus indicating that the convection is dominant mode of charge transport compared to conduction for the present case, essentially due to the higher ambient velocities for low viscosities (low 
$Oh$). For 
$Oh = 0.63$, droplets repel each other for 
$\sigma _r = 1.04$ and 
$\varepsilon _r = 0.2$. However, for 
$Oh = 0.02$, as shown in figure 9, we observe that the droplets initially attract each other and subsequently when the flow has developed and the hydrodynamic forces are stronger, droplets repel each other.
Motion of a pair of droplets for low Ohnesorge number 
$(Oh = 0.02)$ at different non-dimensional time (
$t/\tau _\rho ^\gamma$) for 
$\sigma _r = 1.04$ and 
$\varepsilon _r = 0.2$ using 
$Ca_E = 1.5$, 
$Re_E = 44.131$ and 
$Oh = 0.02$.
Thus, to determine the nature of the interaction between the droplets, it is not sufficient to look at the ratio of the relaxation time constants (or the direction of the electrically driven convective flow) alone. Baygents et al. (Reference Baygents, Rivette and Stone1998) also suggested that, in the low-Reynolds-number limit, the droplets interact hydrodynamically when the separation distance is large and they are essentially driven by electrostatic forces when the separation distance is of the order of two droplet diameters. In what follows, we estimate the hydrodynamic and electrohydrodynamic forces, and use the ratio of the forces to determine the nature of the dynamics of a pair of droplets.
Assuming linear drag law for low-
$Re$ cases discussed here, we can estimate the net hydrodynamic force on the droplet as (cf. § 4.9 in Batchelor Reference Batchelor2000)
\begin{equation} F_h \sim 4 {\rm \pi}\mu_1 U_e R_o \frac{1 + 3/2 \mu_r}{1 + \mu_r} (\boldsymbol{u} - \boldsymbol{v}), \end{equation}
where 
$\boldsymbol {v}$ is the velocity of the centre of mass of the droplet, non-dimensionalised by 
$U_e$, and the magnitude of the ambient flow velocity 
$\boldsymbol {u} \sim (R_o/h)^2 (\sigma _r - \varepsilon _r)/(\sigma _r + 2)^2$ defines the non-dimensionalised velocity of the ambient fluid. The factor 
$(\sigma _r - \varepsilon _r)/(\sigma _r + 2)^2$ is the correction for the velocity scale as discussed in the previous section. Electrostatic force between the droplets can be estimated by considering the interactions between the dipoles defined by the sum of the net free and bound charges on the drop surface. For spherical droplets in an electric field, the total charge in the northern hemisphere is given by
\begin{equation} Q_{net} = 3 {\rm \pi}\varepsilon_0 E_\infty R_o^2 \frac{(\sigma_r - 1)}{(\sigma_r + 2)}, \end{equation}
where 
$\varepsilon _0$ is the permittivity of free space. An estimate of the corresponding electrostatic force between the droplets (dipole–dipole interaction), placed at a centre-to-centre separation distance of 
$h$, can be obtained as
\begin{equation} F_e \sim \frac{2 Q_{net}^2}{4 {\rm \pi}\varepsilon_1} \frac{R_o^2 (R_o^2 - 3 h^2)}{h^2(h^2 - R_o^2)^2}. \end{equation}
The above electrostatic force in the limit of large 
$h$ results in 
$E_\infty ^2/h^4$ dependence discussed in Baygents et al. (Reference Baygents, Rivette and Stone1998). The relative strengths of the electrostatic to hydrodynamic forces can be estimated from the following expression:
\begin{equation} \frac{F_e}{F_h} \sim \frac{9}{8} \left(\frac{\varepsilon_0}{\varepsilon_1}\right)^2 \frac{(\sigma_r-1)^2}{(\sigma_r - \varepsilon_r)} \left(\frac{(1+\mu_r)^2}{2+3\mu_r}\right)\frac{1-3(h/R_o)^2}{((h/R_o)^2-1)^2}. \end{equation}
Thus, from the above expression, the key parameters that govern the relative velocity between the droplets are the ratio of electric conductivities 
$\sigma _r$, initial separation distance between the droplets 
$h/R_o$, the ratio of charge relaxation time 
$\tau _e^2/\tau _e^1 = \varepsilon _r/\sigma _r$ and also the permittivity of the ambient fluid 
$\varepsilon _1$. Permittivity of the ambient fluid has a detrimental effect on the electrostatic interactions due to dielectric screening effects.
Further, response of the droplets to the hydrodynamic forces is also governed by the relaxation time that accounts for the inertia of the droplet. The relaxation time for a droplet in a flow (cf. Batchelor Reference Batchelor2000), 
$\tau _R$, can be written as
\begin{equation} \frac{\tau_R}{\tau_\rho^\gamma} = \frac{2\rho_r}{3 Oh}\frac{(1+\mu_r)}{(2+3\mu_r)}. \end{equation}
In the simulations presented in the current study, we have 
$\rho _r = \mu _r = 1$. Thus, 
$\tau _R/\tau _\rho ^\gamma = 4/(15 Oh)$. This indicates that at low 
$Oh$, a droplet will not relax to the ambient flow immediately, and thus the droplet may first respond to electrostatic attraction before being driven by the hydrodynamic forces as shown in figure 9. For 
$Oh \sim 0.02$, we observe that the droplets initially tend to come together but as they approach each other, the flow becomes stronger and they eventually repel each other. The relaxation time constants 
$t_R/t_\rho ^\gamma$ for 
$Oh = 0.63$ and 
$Oh = 0.02$ are 0.42 and 13.33, respectively. Therefore, for 
$Oh = 0.63$, the droplets quickly relax to the ambient fluid velocity, whereas for 
$Oh = 0.02$, electrostatic forces initially drive the droplet before the hydrodynamic forces set in. Moreover, as discussed earlier, time for flow development 
$\tau _{\mu }/t_\rho ^\gamma$ (
${\sim }1250$ for 
$Oh = 0.02$) is slower for low-viscosity fluids and therefore adds to the effect of slow relaxation time of the droplets. Further, for even lower-viscosity fluids, the droplets can be expected to actually coalesce due to the electrostatic attraction instead of repelling each other due to the hydrodynamic forces (not shown here). Thus, the mutual interaction between the drops of an emulsion is expected to become complex due to the presence of inertial effects in addition to the viscous and electrostatic forces.
For the cases shown in figures 5 and 6, 
$t_R/t_\rho ^\gamma \sim 0.42$ corresponding to a fast deformation of the droplets in response to the electrostatic forces and 
$\tau _{\mu }/t_\rho ^\gamma \sim 39.68$ which corresponds to the time scale of migration of the droplets. The order of this time scale is in good agreement with the observations from the simulations. For these cases, 
$F_e/F_h$ is 
$1.44$ and 
$0.00022$, thus clearly suggesting that electrostatic forces are dominant in the first case for 
$\sigma _r = 6$ and 
$\varepsilon _r = 8$, whereas for 
$\sigma _r = 1.04$ and 
$\varepsilon _r = 0.2$ hydrodynamic forces are dominant. Thus, it is expected that the velocity of attraction will be proportional to 
$1/h^4$ for the former case and to 
$1/h^2$ for the latter. We present in figure 10 the variation in the relative centre-of-mass velocity with the separation distance between the droplets on a log–log scale. As expected, relative velocity between the droplets varies as 
$1/h^4$ for approaching droplets as shown in figure 10(a) and as 
$1/h^2$ for repelling droplets as shown in figure 10(b).
Velocity variation with respect to centre-to-centre distance between the drops for (a) 
$\sigma _r = 6$ and 
$\varepsilon _r = 8$ 
$(Re_E = 4.86)$ and (b) 
$\sigma _r = 1.04$ and 
$\varepsilon _r = 0.2$ 
$(Re_E = 2.223)$ at 
$Ca_E = 1.5$ and 
$Oh = 0.63$.
Other possible interesting behaviour is also observed for a certain set of parameters. For instance, we show that the droplets may approach each other, come close, but not coalesce and instead form a doublet as shown in figure 11. For 
$\sigma _r = 25$ and 
$\varepsilon _r = 2$, when the droplets are initially at a separation of 
$h/R_o = 4$, we obtain 
$F_e/F_h \sim 4.7$ where the electrostatic forces are attractive but the hydrodynamic forces are repulsive. Since electrostatic forces are stronger, droplets attract but the drainage of the film between the droplets is delayed leading to the formation of a drop–drop doublet the drops of which being separated by a thin film of the ambient fluid as shown in the zoomed-in view in figure 11. The electric Reynolds number 
$Re_E = 0.092$ is much smaller compared to the other cases presented in this study. A low 
$Re_E$ indicates much stronger influence of the electrostatic forces compared to the hydrodynamic drag which should have led to coalescence of the droplets. However, due to the nature of the hydrodynamic flow around the droplets, a thin film of the ambient fluid is captured. To resolve the thin film, we use a more refined grid with 
$\varDelta x_{min} = R_0 /410$ for the current case. The gap at 
$t/\tau _\rho ^\gamma \sim 44$ is less whereas a dimple forms at the centre of the droplets around 
$t/\tau _\rho ^\gamma \sim 50$ and remains stable at later times, as shown in figure 11 for 
$t/\tau _\rho ^\gamma \sim 70$. We note that once the flat film is fully formed at 
$t/\tau _\rho ^\gamma \sim 50$, further evolution leads to the formation of a dimple essentially due to the drop shape relaxation during which the thin film is pulled together, thus resulting in the formation of the dimple. The minimum thickness of the film is towards the periphery (
${\sim }0.0073R_o$) and in the dimple region it is 
$0.061R_o$. The streamlines shown in figure 11 indicate that the ambient fluid is driven towards the thin film from the surface of the droplets essentially slowing the drainage of the film. In the simulations, we find that these doublets do not change significantly with time indicating a steady-state behaviour in the time scales of interest. Since we are not incorporating the effect of van der Waals forces, the drainage time is dictated by the slow viscous time scale which in this case is further slowed down due to the ambient fluid flow towards the film. The formation of a dimple in the centre is similar to that observed during coalescence of bubbles or the impact of a droplet on a flat surface. The above phenomenon of non-coalescence and formation of a drop–drop doublet is different from the non-coalescence phenomenon observed for charged droplets, where charged droplets initially approach each other, form a conical neck, exchange charge and then repel each other (Aryafar & Kavehpour Reference Aryafar and Kavehpour2007; Ristenpart et al. Reference Ristenpart, Bird, Belmonte, Dollar and Stone2009; Anand et al. Reference Anand, Juvekar and Thaokar2019; Sunder & Tomar Reference Sunder and Tomar2020). Although there is no direct mention of the above phenomenon in experimental observations, there is evidence in certain experiments. For example, figures 3 and 4(a) of Zimmermann & Vienken (Reference Zimmermann and Vienken1982) show two vesicles forming doublets during cell–cell fusion in the presence of electric field. Also, figure 5 of Holto et al. (Reference Holto, Berg and Lundgaard2009) suggests the formation of a chain of water droplets in an oil and indicating that the droplets align with the electric field but do not coalesce in spite of being in the immediate neighbourhood of each other. We note that the electric stability of the ambient medium will also play a role in determining the complete physics of interaction, since electric field intensity between the droplets is high which can lead to an electric breakdown leading to a more complex phenomenon involving local heating as discussed in the case of particles by Arp & Mason (Reference Arp and Mason1977). We hope that this observation would motivate some careful experiments for drop–drop coalescence in this regime.
Enlarged view of a pair of non-coalescing droplets at different non-dimensional time (
$t/\tau _\rho ^\gamma$) for 
$\sigma _r = 25$ and 
$\varepsilon _r = 2$ 
$(Re_E=0.092)$ at 
$Ca_E = 0.2$ and 
$Oh = 0.5$ and the flow patterns around the drop. Also shown are the zoomed-in views (representing the same area) at 
$t/\tau _\rho ^\gamma = 44$ and 70 and streamlines at 
$t/\tau _\rho ^\gamma = 70$.
3.2. Double emulsion in an externally applied electric field
In this section, we first simulate the deformation of a single compound droplet, constituting a double emulsion in an externally applied electric field and compare the results with those of Abbasi et al. (Reference Abbasi, Song, Kim and Lee2017). As shown in figure 12, the results are in good agreement for the deformation of the constituent droplets. A compound droplet with shell radius 
$R_o$ and inner droplet radius 
$R_i$ is suspended in a continuous medium and an externally applied electric field acts along the horizontal direction. The core radius 
$R_i$ is equal to 
$\beta R_o$. An axisymmetric computational domain of size (
$3R_o\times 6R_o$) is chosen with the compound drop at the centre of the domain. An external electric field is applied such that the right-hand boundary of the domain is considered at a higher potential and the left-hand boundary is grounded as discussed earlier in § 2. The deformation of the outer droplet is similar to that expected for single droplets for low aspect ratios 
$\beta$. For the inner droplet, the deformation will depend on the reduced electric field intensity in the outer droplet and the ratio of the permittivities and conductivities. However, for higher aspect ratios some deviation is expected (see Abbasi et al. Reference Abbasi, Song, Kim and Lee2017). The stability of a compound droplet was investigated numerically by Abbasi et al. (Reference Abbasi, Song, Kim and Lee2017), where it was suggested that if 
$\varepsilon _r > \sigma _r$, the compound drop is unstable and the inner droplet would exit the outer droplet, whereas if 
$\varepsilon _r < \sigma _r$ the inner droplet migrates to the centre of the outer droplet. Abbasi et al. (Reference Abbasi, Song, Kim and Lee2017) also investigated droplet bursting due to excessive deformation of the inner droplet, especially when the aspect ratio 
$\beta$ is large. In figure 12(a), 
$\varepsilon _r$ is kept constant at 2 and simulations are performed over a range of 
$\sigma _r$ varying from 0.1 to 20. When 
$\sigma _r < \varepsilon _r$, we can observe that both the interfaces experience large deformation. When 
$\sigma _r > \varepsilon _r$, the magnitude of deformation at both the interfaces decreases which is depicted by the flattening of the curves. In figure 12(b), 
$\sigma _r$ is kept constant at 2.5 and simulations are performed over a range of 
$\varepsilon _r$ varying from 0.01 to 50. When 
$\varepsilon _r < \sigma _r$ (for a sufficiently large 
$\sigma _r$ for a given 
$\varepsilon$ as discussed previously in the case of an isolated single droplet), the core shows oblate deformation whereas the shell shows prolate deformation. For 
$\varepsilon _r > \sigma _r$, the shell now shows oblate deformation and the core shows prolate deformation. We note that the deformations of the inner and outer droplets obtained from the simulations presented in this study are in good agreement with those of the simulations of Abbasi et al. (Reference Abbasi, Song, Kim and Lee2017).
Comparison of deformation of the core and the shell of a compound droplet at 
$Ca_E = 0.2$ and 
$Oh = 1.41$ for (a) 
$\varepsilon _r = 2$ and (b) 
$\sigma _r = 2.5$.
Interactions between compound droplets in a double emulsion can lead to various interesting configurations for different values of 
$\sigma _r$ and 
$\varepsilon _r$. To investigate this, we consider a pair of compound droplets, each with outer radius 
$R_o$ and inner droplet radius 
$R_i$, separated by a centre-to-centre distance of 
$4R_o$. The shell is three times the size of the core (
$\beta = 1/3$). An axisymmetric computational domain of size (
$5R_o \times 10R_o$) subjected to an externally applied electric field along the axial direction is considered (see figure 1). The flow parameters used for the simulations are 
$Ca_E = 0.2$ and 
$Oh = 0.5$. In figure 13(a), for 
$\sigma _r = 11$ and 
$\varepsilon _r = 16.5$ (
$\sigma _r < \varepsilon _r$), droplets attract and come together and so do the inner droplets, leading to the formation of a single larger compound droplet. The outer droplet (shell) deforms prolately, whereas the deformation of the inner droplet is small. Electric field intensity inside the outer droplet is 
${\sim }3 \boldsymbol {E}_\infty /(\sigma _r + 2)$ and thus, for higher 
$\sigma _r$, it is expected to reduce significantly. The effective 
$Ca_E$ for the inner droplets is thus much smaller due to the reduced electric field intensity and increased capillary pressure because of smaller radius. Therefore, the inner droplets deform significantly less. Interestingly, during coalescence the inner droplets approach faster than the centroid of the outer droplet. Since the ratio of the electrostatic to hydrodynamic forces is 
$F_e/F_h \sim 3.42$, electrostatic forces dominate and drive the droplets together, and thus also drag the inner droplets along. Subsequently, hydrodynamic forces on the inner droplets drag them further towards the periphery. This suggests the possibility of the inner droplet escaping the outer droplet as is shown later for 
$\sigma _r = 10$ and 
$\varepsilon _r = 27.5$. For 
$\sigma _r = 3$ and 
$\varepsilon _r = 0.15$ (
$\sigma _r > \varepsilon _r$), shown in figure 13(b), the outer droplets move apart at a very slow rate, and therefore do not show much variation in terms of droplet shapes with time for the duration investigated in this study. Further, since 
$\sigma _r > \varepsilon _r$, the inner droplets are also expected to remain at the centre of the compound droplets. This is in agreement with the observations from the simulations of Abbasi et al. (Reference Abbasi, Song, Kim and Lee2017), where for 
$\sigma _r < \varepsilon _r$, an off-centred inner droplet is driven by a flow in the outer droplet from the centre towards the pole, and the inner droplet is ejected out of the outer droplet. On the other hand, if the flow generated in the outer droplet is such that the flow is from the poles towards the centre along the centreline of the droplet (as observed for 
$\sigma _r > \varepsilon _r$), the inner droplet would be pushed back to the centre of the outer droplet. Figure 14 shows the variation in the relative velocity of the compound droplets with the separation distance. Oscillations in the relative velocity of the outer droplet as well as for the inner droplets are observed in the repulsion case, as shown in figure 14(b). Since the magnitudes of the velocities are relatively small, the oscillations are more prominently visible in the repulsion case. Since 
$F_e/F_h$ is 3.42 for 
$\sigma _r = 10$ and 
$\varepsilon _r = 27.5$, a scaling of 
$1/h^4$ is found for the approach velocity of the droplets (see figure 14a); whereas for 
$\sigma _r = 3$ and 
$\varepsilon _r = 0.15$, 
$F_e/F_h\sim 0.26$, and droplets repel each other, governed essentially by the hydrodynamic forces resulting in a relative velocity that decreases as 
$1/h^2$ as shown in figure 14(b).
Drop deformation and translation at different non-dimensional time (
$t/\tau _\rho ^\gamma$) for (a) 
$\sigma _r = 11$ and 
$\varepsilon _r = 16.5$ 
$(Re_E=1.167)$ and (b) 
$\sigma _r = 3$ and 
$\varepsilon _r = 0.15$ 
$(Re_E=0.23)$ using the flow parameters 
$Ca_E = 0.2$ and 
$Oh = 0.5$. We also plot the flow patterns in and around the core and the shell.
Velocity variation with respect to centre-to-centre distance between the drops for (a) 
$\sigma _r = 11$ and 
$\varepsilon _r = 16.5$ 
$(Re_E=1.167)$ and (b) 
$\sigma _r = 3$ and 
$\varepsilon _r = 0.15$ 
$(Re_E=0.23)$ at 
$Ca_E = 0.2$ and 
$Oh = 0.5$. The motion of the inner core droplet is also depicted using a dash-dotted line.
For 
$\sigma _r = 6$ and 
$\varepsilon _r = 45$, as shown in figure 15(a), the inner droplets escape the compound droplets from the rear ends, as the two compound droplets approach each other. Whereas, for 
$\sigma _r = 10$ and 
$\varepsilon _r = 27.5$, the inner droplets are released from the approaching ends of the two compound droplets before they coalesce (see figure 15b), thus resulting in the breakup of the double emulsion into a single emulsion, which as expected for 
$\sigma _r < \varepsilon _r$ would result in an unstable emulsion in the presence of an electric field. For 
$\sigma _r = 6$ and 
$\varepsilon _r = 45$, 
$F_e/F_h \sim 0.12$; whereas for 
$\sigma _r = 10$ and 
$\varepsilon _r = 27.5$, 
$F_e/F_h \sim 0.87$. Therefore, in the former case, electrostatic forces are much weaker compared to hydrodynamic forces and the small inertial lag of the inner droplet leads to a biased ejection of the inner droplet from the portion of the periphery of the droplet that faces away from the other compound droplet. On the other hand, for 
$\sigma _r = 10$ and 
$\varepsilon _r = 27.5$, electrostatic forces drive the droplets together initially, but are not sufficiently strong as was the case for 
$\sigma _r = 11$ and 
$\varepsilon _r = 16.5$, and therefore the inner droplets eject from the portion of the outer droplets facing each other. We note that the hydrodynamic response of the droplets can be weakened by reducing 
$Oh$ leading to a possible entrapment of the inner droplets even for the case with 
$\sigma _r = 10$ and 
$\varepsilon _r = 27.5$.
Drop deformation and translation at different non-dimensional time (
$t/\tau _\rho ^\gamma$) for (a) 
$\sigma _r = 6$ and 
$\varepsilon _r = 45$ 
$(Re_E=5.257)$ and (b) 
$\sigma _r = 10$ and 
$\varepsilon _r = 27.5$ 
$(Re_E=2.073)$ using the flow parameters 
$Ca_E = 0.2$ and 
$Oh = 0.5$. We also plot the flow patterns in and around the core and the shell.
Finally, we present three cases showing interesting dynamics of the droplets. In the first one, the outer droplets approach each other but do not coalesce but form a drop–drop doublet, as shown in figure 16(a) for 
$\sigma _r = 25$ and 
$\varepsilon _r = 2$. As discussed earlier for the single droplets with the same 
$\sigma _r$ and 
$\varepsilon _r$, this interesting physics is due to the strong initial electrostatic attraction and subsequent delay in the drainage of the thin film between the two droplets due to the nature of the flow around the droplets as shown by the streamlines in figure 16(a) for 
$t/\tau _\rho ^\gamma = 100$. Since electrostatic forces are much stronger than hydrodynamic forces as suggested by 
$F_e/F_h \sim 4.71$, droplets are initially driven together and the flow around the droplets is developed by the motion of the droplets as indicated by the streamlines at 
$t/\tau _\rho ^\gamma = 43$, but the flow lines reverse after the development of the flow generated by the tangential electrical stresses as shown by the streamlines at 
$t/\tau _\rho ^\gamma = 100$. However, for 
$\sigma _r = 30$ and 
$\varepsilon _r = 10$, 
$F_e/F_h \sim 7.9$, indicating an even stronger electric field that drives the droplets to approach each other and eventually coalesce (see figure 16b). These cases clearly indicate that the simple rule depending on the sign of 
$(\sigma _r - \varepsilon _r)$ for determining the interactions between the droplets is insufficient and the relative strengths of electrostatic and hydrodynamic forces, along with inertial effects governed by 
$Oh$, need to be accounted for. Interestingly, the inner droplets for 
$\sigma _r = 30$ and 
$\varepsilon _r = 10$ do not coalesce even when simulations are run up to 
$t/\tau _\rho ^\gamma \sim 400$. For 
$\sigma _r = 7$ and 
$\varepsilon _r = 10$ shown in figure 16(c), the inner droplets form a configuration similar to that in figure 16(b) but eventually the inner droplets repel (since 
$\sigma _r < \varepsilon _r$) and are ejected.
Drop deformation and translation at different non-dimensional time (
$t/\tau _\rho ^\gamma$) for (a) 
$\sigma _r = 25$ and 
$\varepsilon _r = 2$ 
$(Re_E=0.092)$, (b) 
$\sigma _r = 30$ and 
$\varepsilon _r = 10$ 
$(Re_E=0.284)$ and (c) 
$\sigma _r = 7$ and 
$\varepsilon _r = 10$ 
$(Re_E=1.1)$ using the flow parameters 
$Ca_E = 0.2$ and 
$Oh = 0.5$. We also plot the flow patterns in and around the core and the shell.
4. Summary
In the present study, we have numerically investigated the interaction dynamics of a pair of droplets in a uniform electric field. A leaky dielectric model is employed for this study with the inner and the continuous phase being the same, and the outer shell of the compound droplet being different. We note that the stability of a compound droplet as well the relative motion between a pair of droplets depend essentially on the charge relaxation time constant ratio of the ambient and the outer shell of the compound droplet (
$\tau _e^2/\tau _e^1$) or equivalently the ratio 
$\sigma _r/\varepsilon _r$. This parameter describes the sign of the charge that develops at the droplet interface and thus also the sense of the circulation of the flow generated due to the tangential electric stresses. If the flow is from poles to equator in the ambient fluid, it implies that the flow in the outer shell is from the centre towards the pole. The flow generated in the ambient indicates a hydrodynamic drag on the neighbouring droplet thus driving the two droplets together and making the emulsion unstable. Also, due to the flow inside the outer shell of the compound droplet being from the centre towards the pole, any deviation of the inner droplet from the centre of the compound droplet would result in a hydrodynamic drag that will drive the inner droplet towards the ambient fluid. We show that this simple rule of thumb is not useful when the inertial forces need to be accounted for. Despite the key role played by the hydrodynamic forces governed by the strength of the electric field 
$Ca_E$ and the ratio of charge relaxation time constants 
$\tau _e^2/\tau _e^1$, we observed that for lower values of Ohnesorge number, 
$Oh$, electrostatic dipole–dipole forces can modify certain features of the interaction dynamics. We showed that for low 
$Oh$, droplets that are expected to repel each other due to the hydrodynamic drag instead approach each other due to the electrostatic interactions. Nevertheless, after initial approach, the droplets eventually repel due to the hydrodynamic forces. This clearly indicates that the a priori expectation of the stability of an emulsion in an electric field would require consideration of the relative time scales of electrostatic attraction, charge relaxation and flow time scale.
Moreover, the relative strength of the electrostatic and hydrodynamic forces, given by the ratio 
$F_e/F_h$ in (3.5), also assists in understanding the interaction dynamics of the droplets. For low 
$F_e/F_h$, hydrodynamic forces govern the motion of the droplets and the velocity of approach varies as 
$1/h^2$, whereas for higher 
$F_e/F_h$, electrostatic forces are dominant and the velocity of approach varies as 
$1/h^4$. We have shown that for certain parameters, where repulsion between the droplets is expected for 
$\sigma _r/\varepsilon _r > 1$, electrostatic forces can lead to the droplets coming together and coalescing. We believe that this numerical study would motivate some more experiments in regimes that show interesting dynamics of droplets such as the formation of a drop–drop doublet.
Funding
This research received no specific grant from any funding agency, commercial or not-for-profit sectors.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The corresponding author may be contacted for details of the Basilisk flow solver (open source: http://basilisk.fr/) implementation of the computational work.
Author contributions
G.T. derived the theory and S.K.D. performed the simulations. All authors contributed equally to analysing data and reaching conclusions, and to writing the paper.
Open access
