A three-dimensional small-deformation theory for electrohydrodynamics of dielectric drops

Abstract

Electrohydrodynamics of drops is a classic fluid mechanical problem where deformations and microscale flows are generated by application of an external electric field. In weak fields, electric stresses acting on the drop surface drive quadrupolar flows inside and outside and cause the drop to adopt a steady axisymmetric shape. This phenomenon is best explained by the leaky-dielectric model under the premise that a net surface charge is present at the interface while the bulk fluids are electroneutral. In the case of dielectric drops, increasing the electric field beyond a critical value can cause the drop to start rotating spontaneously and assume a steady tilted shape. This symmetry-breaking phenomenon, called Quincke rotation, arises due to the action of the interfacial electric torque countering the viscous torque on the drop, giving rise to steady rotation in sufficiently strong fields. Here, we present a small-deformation theory for the electrohydrodynamics of dielectric drops for the complete Melcher–Taylor leaky-dielectric model in three dimensions. Our theory is valid in the limits of strong capillary forces and highly viscous drops and is able to capture the transition to Quincke rotation. A coupled set of nonlinear ordinary differential equations for the induced dipole moments and shape functions are derived whose solution matches well with experimental results in the appropriate small-deformation regime. Retention of both the straining and rotational components of the flow in the governing equation for charge transport enables us to perform a linear stability analysis and derive a criterion for the applied electric field strength that must be overcome for the onset of Quincke rotation of a viscous drop.

1. Introduction

When applied to fluid media and interfaces, electric fields give rise to a wide variety of phenomena relevant for both fundamental research and industrial applications. The study of such phenomena has come to be known as electrohydrodynamics (EHD) and has a rich history in the fluid mechanics literature (Melcher & Taylor 1969; Saville 1997). The first study on the subject can be traced back to the work of Rayleigh (1882) on the stability of charged raindrops. Building on his pioneering work, Zeleny (1915, 1917) and Wilson & Taylor (1925) performed experiments and theoretical analyses on the stability of charged prolate-shaped droplets as well as uncharged soap bubbles in electric fields. A few decades later, Taylor (1964) performed an ingenious yet simple analysis showing that a conical interface between a conducting liquid and a dielectric medium can exist in equilibrium in an electric field but only when the cone has a semi-vertical angle of $49.3^\circ$ (Fernández de La Mora 2007; Collins et al. 2008, 2013). Under certain conditions, such Taylor cones on the tips of prolate drops can emit extremely fine jets. This EHD phenomenon forms the basis of electrospraying techniques used in mass spectrometry (Fenn et al. 1989). Early models for the deformation of drops assumed a balance of normal electric stresses with hydrodynamic pressure across the interface and were able to explain the steady prolate shapes seen in most experiments. However, other experiments (O'Konski & Harris 1957; Allan & Mason 1962) also reported steady oblate shapes, leading Taylor (1966) to posit the existence of a surface charge density at the drop surface giving rise to tangential interfacial stresses. These stresses then drive a quadrupolar fluid flow inside and outside the drop and are responsible for the reported oblate shapes. The leaky-dielectric model used in this analysis, subsequently formalized by Melcher & Taylor (1969), hypothesizes that the rate of accumulation of surface charge is balanced by Ohmic currents from the bulk and that surface charge is advected by the flow.

In a weak electric field, the drop assumes a steady axisymmetric shape that is either prolate or oblate depending on material properties, and we refer to this as the Taylor regime. However, on increasing the electric field strength, the system can become unstable and spontaneously adopt a non-axisymmetric shape with a tilted dipole moment and steady electro-rotation of the drop (Krause & Chandratreya 1998; Ha & Yang 2000b; Sato et al. 2006; Salipante & Vlahovska 2010). This spontaneous symmetry breaking, known as Quincke rotation, was first observed in solid spheres by Quincke (1896). Even though it is well explained by the Melcher–Taylor leaky-dielectric model proposed in 1969, a theoretical analysis of Quincke rotation was presented much later by Jones (1984). There has been a revival of interest in Quincke rotation of particles (Bricard et al. 2013, 2015; Das & Lauga 2019; Karani, Pradillo & Vlahovska 2019) and drops (Salipante & Vlahovska 2010, 2013; Rozynek, Bielas & Józefczak 2018), fuelled in part by an interest in designing self-propelled particles with controllable collective dynamics and suspensions or emulsions with tunable rheological properties (Pannacci, Lemaire & Lobry 2007). Readers are referred to Vlahovska (2019) and Papageorgiou (2019) for recent reviews on the subject.

While models for rigid particles are well developed (Jones 1984; Das & Saintillan 2013), the case of deformable droplets is significantly more complex due to the nonlinear coupling between deformations, fluid flow and interfacial charge dynamics. In his original analysis, Taylor (1966) performed a first-order small-deformation theory for an axisymmetric dielectric drop in an electric field, neglecting shape and charge transients as well as interfacial charge advection by the flow. Since his pioneering work, there have been several attempts to include various effects such as transient shape deformations (Esmaeeli & Sharifi 2011), transient charge relaxation and fluid acceleration (Lanauze, Walker & Khair 2013), or coupling with other fields such as gravity (Bandopadhyay et al. 2016; Yariv & Almog 2016). On the computational side, several simulation studies have solved the Melcher–Taylor leaky-dielectric model using boundary integral methods (Sherwood 1988; Lac & Homsy 2007; Lanauze, Walker & Khair 2015; Nganguia et al. 2016; Das & Saintillan 2017a) or grid-based methods (López-Herrera, Popinet & Herrada 2011; Hsu, Hu & Lai 2019; Theillard, Gibou & Saintillan 2019) and often match well with experiments (Ha & Yang 2000a,b; Salipante & Vlahovska 2010). It is also worth mentioning a related problem concerning the behaviour and breakup of conducting drops in electric fields (Dubash & Mestel 2007a,b; Karyappa, Deshmukh & Thaokar 2014).

Including charge convection in analytical theories is challenging, as it couples nonlinearly with the interfacial fluid velocity. However, charge convection can be significant in strong electric fields and is critical to include in any analysis of Quincke rotation. Feng (2002) performed an analysis in which he included the effect of both rotational and straining flows on charge transport. However, as his analysis was limited to two dimensions, the drop deformation and tilt angle in the Quincke regime were found not to depend on the viscosity ratio, similar to the Taylor regime, and consequently the critical electric field for the onset of rotation was the same as that for a solid cylindrical particle. Yariv & Frankel (2016) also focused on two dimensions but analysed the asymptotic limit of strong charge convection (large electric Reynolds number). They found that no fore–aft symmetric solutions exist in this limit, implying a transition to spontaneous rotation. Shutov (2002) and Shkadov & Shutov (2002) included charge convection for axisymmetric drops using small-deformation theory. However, they incorrectly neglected it at first order and only included it at the second order. He, Salipante & Vlahovska (2013) performed a three-dimensional analysis of the leaky-dielectric model but only included the effect of solid-body rotation in the charge convection. This assumption also results in the critical electric field being independent of viscosity as in the two-dimensional case. Tyatyushkin (2017) proposed a model for very viscous drops in three dimensions in which he focused on the effect of surface charge conduction and found a modified expression for the critical electric field for Quincke rotation. In our previous work (Das & Saintillan 2017b), we developed a small-deformation theory for the complete Melcher–Taylor leaky-dielectric model, including transient shape deformation, transient charge relaxation and nonlinear charge convection. However, the drop shape was assumed to remain axisymmetric, preventing the occurrence of Quincke rotation. Both experiments (Salipante & Vlahovska 2010) and numerical simulations (Das & Saintillan 2017a) have shown that the critical electric field for the onset of Quincke rotation of a deformable drop deviates from that of a solid sphere, which demands a theoretical analysis.

In this paper, we present a three-dimensional small-deformation theory for EHD of a drop using the complete Melcher–Taylor leaky-dielectric model. The theory is valid in the limits of strong capillary forces and highly viscous drops. The novelty of our work lies in the retention of the transient charge relaxation term and straining component of the flow in the surface charge evolution equation, which renders the governing equations nonlinear but allows for numerical solutions. Furthermore, inclusion of the latter enables us to perform a linear stability analysis to determine the critical electric field for the onset of rotation. The problem definition and its solution are presented in §§ 2 and 3, respectively. We show in § 4 that our theory recovers Quincke rotation of a solid sphere in the limit of infinite viscosity ratio. Results for deformable drops are compared with experiments and discussed in § 5. Conclusions are summarized in § 6.

2. Problem definition

Let us consider an uncharged neutrally buoyant dielectric droplet suspended in a weakly conducting fluid medium as shown in figure 1. The volumes occupied by the droplet and surrounding fluid are denoted by $V^-$ and $V^+$, respectively. In the absence of an electric field, the droplet remains spherical. However, under the influence of an applied field, the shape deviates from a sphere and undergoes small ellipsoidal deformations in the limit of high viscosity and high surface tension. The unit normal vector on the droplet–fluid interface ${S}$ is denoted by $\boldsymbol {n}$ and points into the suspending fluid medium. The electric field may point in any direction, and without loss of generality we choose a coordinate system in which it is aligned with the $y$-axis. Let $(\epsilon ^{\pm },\sigma ^{\pm },\mu ^{\pm })$ be the dielectric permittivities, electric conductivities and dynamic viscosities of the two fluids, respectively. Non-dimensionalization of the material properties yields three dimensionless groups, typically defined as

(2.1a–c)\begin{equation} R=\frac{\sigma^+}{\sigma^-},\quad Q= \frac{\epsilon^-}{\epsilon^+},\quad \lambda =\frac{\mu^-}{\mu^+}. \end{equation}

The low-drop-viscosity limit $\lambda \rightarrow 0$ describes a bubble, whereas $\lambda \rightarrow \infty$ describes a rigid particle. The product $RQ$ can also be interpreted as the ratio of the inner to outer charge relaxation times:

(2.2)\begin{equation} RQ=\frac{\tau^-}{\tau^+}\quad \mbox{where}\ \tau^\pm=\frac{\epsilon^\pm}{\sigma^\pm}. \end{equation}

Figure 1. Problem definition. (a) Prolate- and (b) oblate-shaped drops in the Taylor regime. Quadrupolar flow fields inside and outside the drop advect charges (depicted as $+$ and $-$ symbols) from equator to poles in the prolate case and from poles to equator in the oblate case. Here $V^{\pm }$ denote the exterior and interior domains, respectively, and $(\epsilon ^\pm , \sigma ^\pm ,\mu ^\pm )$ are the corresponding dielectric permittivities, electric conductivities and dynamic viscosities. (c) Tilted drop in the Quincke regime with both quadrupolar and rotational flow fields present. Here $L$ and $B$ denote the longest and shortest dimensions of the drop, while $\phi ^*$ is the angle between the direction of maximum deformation and the $x$-axis. The electric field is applied along the $y$-axis.

Under the assumptions of the Melcher–Taylor leaky-dielectric model (Melcher & Taylor 1969), all charges in the domain are concentrated on the drop–fluid interface, so that the electric potential in both fluid domains is governed by Laplace's equation:

(2.3)\begin{equation} \nabla^{2}\varphi^{\pm}(\boldsymbol{x})=0 \quad \mbox{for}\ \boldsymbol{x}\in V^{\pm}. \end{equation}

The electric potential and the tangential component of the local electric field are continuous across the interface:

(2.4a,b)\begin{equation} [\kern-1pt[ \varphi (\boldsymbol{x})]\kern-1pt] =0 \quad \mbox{and} \quad [\kern-1pt[ \boldsymbol{E}_t(\boldsymbol{x})]\kern-1pt] =\boldsymbol{0} \quad \mbox{for}\ \boldsymbol{x}\in S, \end{equation}

where $\boldsymbol {E}^{\pm }_{t}=({{\boldsymbol{\mathsf{I}}}}-\boldsymbol {nn})\boldsymbol {\cdot } \boldsymbol {E}^{\pm }$ and $\boldsymbol {E}^{\pm }=-\boldsymbol {\nabla } \varphi ^{\pm }$. We have introduced the notation $[\kern-1pt[ \,f(\boldsymbol {x})]\kern-1pt] \equiv f^+(\boldsymbol {x}) - f^-(\boldsymbol {x})$ for any field variable $f(\boldsymbol {x})$ defined on both sides of the interface. The normal component of the electric field $E_{n}^{\pm }=\boldsymbol {n}\boldsymbol {\cdot }\boldsymbol {E}^{\pm }$ undergoes a jump due to the mismatch in electrical properties between the two media (Landau, Lifshitz & Pitaevskiì 1984), resulting in a surface charge distribution $q(\boldsymbol {x})$ related to the normal displacement field by Gauss's law:

(2.5)\begin{equation} q(\boldsymbol{x})=[\kern-1pt[ \epsilon {E}_n(\boldsymbol{x})]\kern-1pt] \quad \mbox{for}\ \boldsymbol{x}\in S. \end{equation}

The surface charge density $q$ obeys a conservation law that prescribes a balance between temporal changes in charge, Ohmic currents from the bulk and charge advection by the fluid flow with velocity $\boldsymbol {v}(\boldsymbol {x})$ on the drop surface. Accordingly, it satisfies the conservation equation

(2.6)\begin{equation} \frac{\partial q}{\partial t} + [\kern-1pt[ \sigma {E}_n]\kern-1pt] +\boldsymbol{\nabla}_{s} \boldsymbol{\cdot} (q\boldsymbol{v})=0 \quad \mbox{for}\ \boldsymbol{x}\in S, \end{equation}

where $\boldsymbol {\nabla }_{s}\equiv ({{\boldsymbol{\mathsf{I}}}}-\boldsymbol {nn})\boldsymbol {\cdot }\boldsymbol {\nabla }$ is the surface gradient operator. Since the droplet and medium are leaky dielectrics, freely moving charged ions in the bulk are assumed to be present in negligible quantities and only surface charges are assumed to exist in the physical domain. A detailed derivation of the leaky-dielectric model from a more generic electrodiffusion model is provided in the recent work of Mori & Young (2018).

The fluid velocity field $\boldsymbol {v}^{\pm }(\boldsymbol {x})$ and corresponding pressure field $p^{\pm }(\boldsymbol {x})$ satisfy the Stokes equations in both fluid domains:

(2.7a,b)\begin{equation} -\mu^{\pm}\nabla^{2}\boldsymbol{v}^{\pm}+\boldsymbol{\nabla} p^{\pm}=\boldsymbol{0}\quad \mbox{and}\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}^{\pm}=0 \quad \mbox{for}\ \boldsymbol{x}\in V^{\pm}. \end{equation}

The fluid velocity is continuous across the interface $S$,

(2.8)\begin{equation} [\kern-1pt[ \boldsymbol{v}(\boldsymbol{x})]\kern-1pt] =\boldsymbol{0} \quad \mbox{for}\ \boldsymbol{x}\in S, \end{equation}

and the shape of the interface, parametrized by an implicit equation $\xi (\boldsymbol {x},t)=0$, deforms as a material surface,

(2.9)\begin{equation} \frac{\mathrm{D}\xi}{\mathrm{D}t}=\frac{\partial \xi}{\partial t}+ \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} \xi =0 \quad \mbox{for}\ \boldsymbol{x}\in S. \end{equation}

The dynamic boundary condition requires that the jumps in electric and hydrodynamic tractions across the interface balance surface tension forces:

(2.10)\begin{equation} [\kern-1pt[\,\boldsymbol{f}_{E}]\kern-1pt] + [\kern-1pt[\,\boldsymbol{f}_{H}]\kern-1pt] =\gamma (\boldsymbol{\nabla}_{s}\boldsymbol{\cdot}\boldsymbol{n})\boldsymbol{n} \quad \mbox{for}\ \boldsymbol{x}\in S.\end{equation}

Here, $\gamma$ is the constant surface tension and $\boldsymbol {\nabla }_{s} \boldsymbol {\cdot }\boldsymbol {n}=2\kappa _{m}$ is twice the mean surface curvature. The jumps in tractions are obtained in terms of the Maxwell stress tensor ${{\boldsymbol{\mathsf{T}}}}_{E}$ and hydrodynamic stress tensor ${{\boldsymbol{\mathsf{T}}}}_{H}$ as

(2.11)\begin{gather} [\kern-1pt[\,\boldsymbol{f}_{E}]\kern-1pt] = \boldsymbol{n}\boldsymbol{\cdot}[\kern-1pt[{{\boldsymbol{\mathsf{T}}}}_{E}]\kern-1pt]= \boldsymbol{n}\boldsymbol{\cdot}[\kern-1pt[ \epsilon (\boldsymbol{EE}-\tfrac{1}{2}E^{2}{{\boldsymbol{\mathsf{I}}}})]\kern-1pt], \end{gather}

(2.12)\begin{gather}[\kern-1pt[\,\boldsymbol{f}_{H}]\kern-1pt] =\boldsymbol{n}\boldsymbol{\cdot}[\kern-1pt[{{\boldsymbol{\mathsf{T}}}}_{H}]\kern-1pt] = \boldsymbol{n}\boldsymbol{\cdot} [\kern-1pt[-p{{\boldsymbol{\mathsf{I}}}}+\mu (\boldsymbol{\nabla} \boldsymbol{v}+\boldsymbol{\nabla}\boldsymbol{v}^{\textrm{T}} )]\kern-1pt]. \end{gather}

The jump in electric tractions can also be expressed as

(2.13)\begin{equation} [\kern-1pt[\,\boldsymbol{f}_{E}]\kern-1pt] =[\kern-1pt[ \epsilon E_{n} ]\kern-1pt] \boldsymbol{E}_{t}+\tfrac{1}{2}[\kern-1pt[ \epsilon(E_{n}^{2}-E_{t}^{2})]\kern-1pt] \boldsymbol{n} =q\boldsymbol{E}_{t}+ [\kern-1pt[ p_{E}]\kern-1pt] \boldsymbol{n}.\end{equation}

The first term on the right-hand side captures the tangential electric force acting on the interfacial charge, while the second term captures normal electric stresses and can be interpreted as a jump in electric pressure $p_E=\frac {1}{2}\epsilon (E_n^2-E_t^2)$ (Lac & Homsy 2007).

The characteristic length and time scales of the problem are the initial radius of the droplet $a_d$ and the Maxwell–Wagner relaxation time $\tau _{MW}$, which is the time scale for polarization of the drop surface upon application of the field (Das & Saintillan 2013):

(2.14)\begin{equation} \tau_{MW}=\frac{\epsilon^- +2\epsilon^+}{\sigma^- +2\sigma^+}. \end{equation}

Non-dimensionalization of the electric tractions with $\epsilon ^+ E_0^2$, hydrodynamic tractions with $\mu ^+ a/\tau _{MW}$ and surface tension forces with $\gamma /a$ gives rise to two additional dimensionless numbers. They consist of the electric capillary number $Ca_E$, denoting the ratio of electric to capillary forces, and the electric Mason number $Ma$, denoting the ratio of viscous to electric stresses:

(2.15a,b)\begin{equation} Ca_{E}=\frac{a_d\epsilon^+ E_{0}^{2}}{\gamma},\quad Ma = \frac{\mu^+}{\epsilon^+ \tau_{MW}E_{0}^2}. \end{equation}

The Mason number is inversely proportional to the electric Reynolds number $Re_E$ sometimes used in the literature (Salipante & Vlahovska 2010; Lanauze et al. 2015; Yariv & Frankel 2016) and defined as

(2.16)\begin{equation} Re_E=\frac{1}{Ma}\frac{1+2R}{R(Q+2)}. \end{equation}

Since we choose $\tau _{MW}$ as the time scale for the problem, the Mason number, $Ma$, appears in the stress balance equation (3.26). Alternatively, one can choose the EHD straining flow time, $\tau _{HD} = \mu ^+(1+\lambda )/(\epsilon ^+ E_0^2)$, as the relevant time scale for the problem, which would result in the Mason number, $Ma$, appearing in the charge conservation equation (3.36), as done in some previous studies (He et al. 2013).

3. Problem solution

We now present the solution to the governing equations introduced above using a small-deformation or domain perturbation theory (Taylor & Acrivos 1964; Matunobu 1966; Barthès-Biesel & Acrivos 1973; Rallison 1980). Note that the solutions are expressed in the laboratory frame of reference.

3.1. Shape of a slightly deformed drop

Let us parametrize the drop–fluid interface as $\xi \equiv r - [1+ \delta f(\theta ,\phi ,t)]=0$, where $(r,\theta ,\phi )$ are spherical coordinates with polar angle $\theta$ and azimuthal angle $\phi$. Departures from the spherical shape are assumed to be small, as evidenced by the small-deformation parameter $\delta$ to be defined more precisely later (see § 3.5). The shape $f$ of the drop is expanded on the basis of spherical harmonics:

(3.1)\begin{equation} f(\theta,\phi,t) = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} \mathcal{L}_{ml}(\cos \theta) [\,f_{ml}(t) \cos m\phi + \tilde{f}_{ml}(t)\sin m\phi], \end{equation}

where $\mathcal {L}_{ml}(\cos \theta )$ are the associated Legendre functions of order $m$ and degree $l$. Owing to the quadratic nature of electric stresses, the applied field induces perturbations in the drop's shape that are of even order (Das & Saintillan 2017b). Retaining terms to the leading order of deformation ($l=2$), we obtain the following expression for the shape:

(3.2)\begin{align} f(\theta,\phi,t)&= \tfrac{1}{2}(3\cos^2\theta - 1)\,f_{02}(t) - 3\cos\theta\sin\theta \, [\,f_{12}(t) \cos\phi + \tilde{f}_{12}(t) \sin\phi] \nonumber\\ &\quad +3\sin^2\theta \, [\,f_{22}(t)\cos 2\phi + \,\tilde{f}_{22}(t) \sin2\phi]. \end{align}

The time-dependent coefficients $\,f_{ml}$ and $\tilde {f}_{ml}$ are unknown and to be determined as part of the solution. Using (3.2), we can obtain the unit normal vector and mean curvature as

(3.3a,b)\begin{equation} \boldsymbol{n} = \left[\frac{\boldsymbol{\nabla} \xi}{|\boldsymbol{\nabla} \xi|}\right]_{r=1}=\boldsymbol{\hat{r}} - \delta \boldsymbol{\nabla} f + {O}(\delta^2),\quad \kappa_m = \tfrac{1}{2} \boldsymbol{\nabla}_s \boldsymbol{\cdot} \boldsymbol{n}. \end{equation}

Detailed expressions for $\boldsymbol {n}$ and $\kappa _m$ are provided in appendix A.

3.2. Electric problem

The solution for the electric potential in a spherical geometry is best written in terms of harmonic multipoles (Jackson 1998), with the leading-order term corresponding to a dipole in this problem. Higher-order multipoles can arise even at first order in deformation due to the nonlinear charge conservation equation, but we neglect them here. The leading order of these higher-order multipoles is ${O}(Ma^{-n} \lambda ^{-n})$ with $n>2$ and they become stronger with decreasing Mason number and viscosity ratio (see §§ 3.6 and 5.2 for details). Our theory is therefore valid for drop–fluid systems in which the straining component of the flow is relatively weak. It is also noted that including higher-order electric multipoles will in turn induce higher-order shape multipoles in (3.2).

The induced dipole moment on the drop is denoted by $\boldsymbol {P} = (P_x,P_y,P_z)$ and is non-dimensionalized by $a^3E_0$. The electric potential, non-dimensionalized by $aE_0$, at location $\boldsymbol {r}=r(\sin \theta \cos \phi ,\sin \theta \sin \phi ,\cos \phi )$ outside and inside the drop is then given by

(3.4a,b)\begin{equation} \varphi^+ = -\boldsymbol{\hat{e}}\boldsymbol{\cdot} \boldsymbol{r} + \frac{\boldsymbol{P}\boldsymbol{\cdot}\boldsymbol{r}}{r^3}, \quad \varphi^- = -\boldsymbol{\hat{e}}\boldsymbol{\cdot} \boldsymbol{r} + \boldsymbol{P}\boldsymbol{\cdot}\boldsymbol{r}, \end{equation}

respectively, where $\boldsymbol {\hat {e}}=\boldsymbol {E}_0/E_0$ denotes the direction of the external field. Equations (3.4a,b) satisfy the boundary conditions (2.4a,b). The normal component of the electric field is discontinuous, while its tangential components are continuous. Note that, in the leading-order perturbation analysis, the boundary conditions are enforced on the surface of the undeformed sphere $r=a$. Therefore, the leading-order electric field components are

(3.5a)\begin{gather} E^+_r = P^+_x\sin\theta\cos\phi + P^+_y \sin\theta\sin\phi + P^+_z \cos \theta, \end{gather}

(3.5b)\begin{gather}E^-_r = P^-_x \sin\theta\cos\phi + P^-_y \sin\theta\sin\phi + P_z^- \cos \theta, \end{gather}

(3.5c)\begin{gather}E_\theta = P^-_x\cos\theta \cos\phi + P^-_y\cos\theta\sin\phi - P^-_z\sin \theta, \end{gather}

(3.5d)\begin{gather}E_\phi = - P^-_x \sin\phi + P^-_y\cos\phi. \end{gather}

We have introduced the short-hand notation $\boldsymbol {P}^+=\boldsymbol {\hat {e}}+2\boldsymbol {P}$ and $\boldsymbol {P}^-=\boldsymbol {\hat {e}}-\boldsymbol {P}$ for the coefficients appearing in the combined electric field accounting for the applied field and induced dipole moment. Doing so allows us to keep the direction $\boldsymbol {\hat {e}}$ of the applied field arbitrary, only to be specified when dealing with the charge conservation equation in § 3.6. For example, if the electric field is applied in the $y$-direction, we substitute $P^+_y = 1 + 2P_y$, $P^-_y = 1 - P_y$, $P^+_{x,z} = 2P_{x,z}$ and $P^-_{x,z} = -P_{x,z}$.

The surface charge density, given by Gauss's law, and the jump in Ohmic current at the interface can both also be expressed in terms of surface harmonics as

(3.6)\begin{gather} q = E^+_r - Q E^-_r = q_{01} \cos \theta+ q_{11} \sin\theta\cos\phi + \tilde{q}_{11} \sin\theta\sin\phi, \end{gather}

(3.7)\begin{gather}\boldsymbol{n} \boldsymbol{\cdot} [\kern-1pt[\,\boldsymbol{j} ]\kern-1pt] = RE^+_r - E^-_r = j_{n,01} \cos \theta + j_{n,11} \sin\theta\cos\phi + \tilde{j}_{n,11} \sin\theta\sin\phi, \end{gather}

where the coefficients $q_{ml}$ and $j_{n,ml}$ are linear functions of the dipole moment:

(3.8a–c)\begin{gather} q_{01} = P_z^+ - Q P_z^-, \quad q_{11} = P_x^+ - Q P_x^-, \quad \tilde{q}_{11} = P_y^+ - Q P_y^-, \end{gather}

(3.9a–c)\begin{gather}j_{n,01} = P_z^+ - Q P_z^-, \quad j_{n,11} = P_x^+ - Q P_x^-, \quad \tilde{j}_{n,11} = P_y^+ - Q P_y^-. \end{gather}

The radial, polar and azimuthal components of the jump in electric tractions in (2.13) can similarly be written as

(3.10a)\begin{align} [\kern-1pt[ p_E ]\kern-1pt] &= [\kern-1pt[ p_E ]\kern-1pt]_{00} + [\kern-1pt[ p_E ]\kern-1pt]_{02} \cos^2\theta + ([\kern-1pt[ p_E ]\kern-1pt]_{12} \cos \phi + [\kern-1pt[ \widetilde{p_E} ]\kern-1pt]_{12} \sin \phi ) \sin 2\theta \nonumber\\ &\quad + ([\kern-1pt[ p_E ]\kern-1pt]_{22} \cos 2\phi + [\kern-1pt[ \widetilde{p_E} ]\kern-1pt]_{22} \sin2\phi ) \sin^2\theta, \end{align}

(3.10b)\begin{align} qE_\theta &= qE_{\theta ,02} \sin 2\theta+qE_{\theta ,10} \cos \phi + \widetilde{q E}_{\theta ,10} \sin \phi \nonumber\\ &\quad + (qE_{\theta ,12} \cos\phi + \widetilde{q E}_{\theta ,12} \sin \phi) \cos 2\theta + (qE_{\theta ,22} \cos 2\phi + \widetilde{q E}_{\theta ,22} \sin 2\phi) \sin 2\theta \end{align}

and

(3.10c)\begin{align} qE_\phi &= qE_{\phi ,01} \sin \theta + (qE_{\phi ,11} \cos\phi + \widetilde{q E}_{\phi ,11} \sin \phi ) \cos \theta \nonumber\\ &\quad + (qE_{\phi ,21} \cos 2\phi + \widetilde{q E}_{\phi ,21} \sin 2\phi) \sin \theta. \end{align}

The unknown coefficients in this case are quadratic functions of the dipole moments and are provided in appendix C. It is instructive to notice that the radial, polar and azimuthal components of the tractions contain harmonic functions of degrees two, two and one, and one, respectively. We also note that the terms containing $qE_{\theta ,10}$, $\widetilde {q E}_{\theta ,10}$ and $qE_{\phi ,01}$ in the tangential tractions capture the net electric torque acting on the droplet.

3.3. Flow problem

We solve the Stokes equations (2.7a,b) in spherical coordinates using Lamb's general solution (Happel & Brenner 1965; Kim & Karrila 2013). The pressure $p$ satisfies Laplace's equation and constitutes the particular solution for the Stokes momentum equation, while the homogeneous solution consists of the velocity potential $\varPhi$ and toroidal flow field $\chi$. Inside the drop, only growing harmonics are retained,

(3.11a)\begin{gather} p_l = r^l \sum_{m=0}^l \mathcal{L}_{ml} (\cos\theta)[a_{ml}\cos m\phi + \tilde{a}_{ml}\sin m\phi], \end{gather}

(3.11b)\begin{gather}\varPhi_l = r^l \sum_{m=0}^l \mathcal{L}_{ml} (\cos\theta)[b_{ml}\cos m\phi + \tilde{b}_{ml}\sin m\phi], \end{gather}

(3.11c)\begin{gather}\chi_l = r^l \sum_{m=0}^l \mathcal{L}_{ml} (\cos\theta)[c_{ml}\cos m\phi + \tilde{c}_{ml}\sin m\phi], \end{gather}

while decaying harmonics are used outside,

(3.12a)\begin{gather} p_{-l-1} = r^{-l-1} \sum_{m=0}^l \mathcal{L}_{ml} (\cos\theta)[A_{ml}\cos m\phi + \tilde{A}_{ml}\sin m\phi], \end{gather}

(3.12b)\begin{gather}\varPhi_{-l-1} = r^{-l-1} \sum_{m=0}^l \mathcal{L}_{ml} (\cos\theta)[B_{ml}\cos m\phi + \tilde{B}_{ml}\sin m\phi], \end{gather}

(3.12c)\begin{gather}\chi_{-l-1} = r^{-l-1} \sum_{m=0}^l \mathcal{L}_{ml} (\cos\theta)[C_{ml}\cos m\phi + \tilde{C}_{ml}\sin m\phi]. \end{gather}

The velocity fields inside and outside the drop are written as

(3.13)\begin{equation} \begin{aligned} \boldsymbol{v}^- & = \sum_{l=1}^\infty \left[\frac{(l+3)r^2\boldsymbol{\nabla}p_l - 2lp_l\boldsymbol{r}}{2\lambda(l+1)(2l+3)} + \boldsymbol{\nabla}\varPhi_l + \boldsymbol{\nabla}\chi_l \times \boldsymbol{r} \right], \\ \boldsymbol{v}^+ & = \sum_{l=1}^\infty \left[\frac{-(l-2)r^2\boldsymbol{\nabla}p_{-l-1} + 2(l+1)p_{-l-1}\boldsymbol{r}}{2 l(2l-1)} \right] + \sum_{l=0}^\infty \left[\boldsymbol{\nabla}\varPhi_{-l-1} + \boldsymbol{\nabla}\chi_{-l-1} \times \boldsymbol{r} \right]. \end{aligned} \end{equation}

and the corresponding pressure fields are $p^- = \sum _{l=0}^{\infty } p_l$ and $p^+ = \sum _{l=-1}^{\infty } p_{-l-1}$. The velocity and pressure fields are then used to obtain the hydrodynamic tractions inside and outside,

(3.14)\begin{equation} \boldsymbol{f}_H^\pm = -\boldsymbol{\hat{r}}p^\pm + \mu^\pm \boldsymbol{\hat{r}} \boldsymbol{\cdot} (\boldsymbol{\nabla} \boldsymbol{v}^\pm + \boldsymbol{\nabla} \boldsymbol{v}^{\pm \textrm{T}}) + {O}(\delta), \end{equation}

which must balance electric stresses as well as capillary forces. Inspecting the radial and tangential components of these various tractions, we can deduce the harmonics that need to be retained:

(3.15a)\begin{gather} \boldsymbol{f}_H^- =-\boldsymbol{\hat{r}} p_0^- + \tfrac{8}{21} r\boldsymbol{\nabla} p_2 - \tfrac{19}{21} \boldsymbol{\hat{r}} p_2 + 2\lambda r^{-1} \boldsymbol{\nabla}\varPhi_2, \end{gather}

(3.15b)\begin{gather}\boldsymbol{f}_H^+ =-\boldsymbol{\hat{r}} p_0^+ + \tfrac{1}{2} r \boldsymbol{\nabla} p_{-3} - \tfrac{3}{2} \boldsymbol{\hat{r}}p_{-3} - 8 r^{-1} \boldsymbol{\nabla}\varPhi_{-3} - 3 \boldsymbol{\nabla}\chi_{-2} \times \boldsymbol{\hat{r}}. \end{gather}

Using the identity $\boldsymbol {r} \boldsymbol {\cdot } \boldsymbol {\nabla } p_l = lp_l$, we find the radial, polar and azimuthal components of the jump in hydrodynamic tractions to be

(3.16a)\begin{gather} \boldsymbol{\hat{r}} \boldsymbol{\cdot} [\kern-1pt[\,\boldsymbol{f}_H ]\kern-1pt] = - [\kern-1pt[ p_0 ]\kern-1pt] - 3 p_{-3} + 24 \varPhi_{-3} + \tfrac{1}{7} p_2 - 4\lambda \varPhi_2, \end{gather}

(3.16b)\begin{gather}\boldsymbol{\hat{\theta}} \boldsymbol{\cdot} [\kern-1pt[\,\boldsymbol{f}_H ]\kern-1pt] = \tfrac{1}{2} \boldsymbol{\hat{\theta}} \boldsymbol{\cdot} \boldsymbol{\nabla} p_{-3} - (8\boldsymbol{\hat{\theta}} \boldsymbol{\cdot} \boldsymbol{\nabla}\varPhi_{-3} + 3 \boldsymbol{\hat{\phi}} \boldsymbol{\cdot} \boldsymbol{\nabla}\chi_{-2}) - \tfrac{8}{21} \boldsymbol{\hat{\theta}} \boldsymbol{\cdot} \boldsymbol{\nabla} p_2 - 2 \lambda \boldsymbol{\hat{\theta}} \boldsymbol{\cdot} \boldsymbol{\nabla}\varPhi_2, \end{gather}

(3.16c)\begin{gather}\boldsymbol{\hat{\phi}} \boldsymbol{\cdot} [\kern-1pt[\,\boldsymbol{f}_H ]\kern-1pt] = \tfrac{1}{2} \boldsymbol{\hat{\phi}} \boldsymbol{\cdot} \boldsymbol{\nabla} p_{-3} - (8 \boldsymbol{\hat{\phi}} \boldsymbol{\cdot} \boldsymbol{\nabla}\varPhi_{-3} - 3 \boldsymbol{\hat{\theta}} \boldsymbol{\cdot} \boldsymbol{\nabla}\chi_{-2}) - \tfrac{8}{21} \boldsymbol{\hat{\phi}} \boldsymbol{\cdot} \boldsymbol{\nabla} p_2 - 2 \lambda \boldsymbol{\hat{\phi}} \boldsymbol{\cdot} \boldsymbol{\nabla}\varPhi_2. \end{gather}

The radial stresses are conveniently expressed in terms of spherical harmonics of second degree $D$, while for the tangential stresses it is convenient to assume $({{\boldsymbol{\mathsf{I}}}}-\boldsymbol {rr}) \boldsymbol {\cdot } [\kern-1pt[ \,\boldsymbol {f}_H ]\kern-1pt] =({{\boldsymbol{\mathsf{I}}}}-\boldsymbol {rr}) \boldsymbol {\cdot } (\boldsymbol {\nabla } G + \boldsymbol {\nabla } F \times \boldsymbol {\hat {r}})$, where $G$ and $F$ are spherical harmonics of second and first degree, respectively. The hydrodynamic stresses are then compactly written as

(3.17a)\begin{align} \boldsymbol{\hat{r}} \boldsymbol{\cdot} [\kern-1pt[\,\boldsymbol{f}_H ]\kern-1pt] &= -[\kern-1pt[ p_0 ]\kern-1pt] + [\tfrac{1}{2}(3\cos^2\theta - 1)D_{02} - 3\cos\theta\sin\theta\, (D_{12} \cos\phi + \tilde{D}_{12} \sin\phi) \nonumber\\ &\quad +3\sin^2\theta\,(D_{22} (\cos 2\phi) + \tilde{D}_{22} \sin(2\phi))], \end{align}

(3.17b)\begin{align} \boldsymbol{\hat{\theta}} \boldsymbol{\cdot} [\kern-1pt[\,\boldsymbol{f}_H ]\kern-1pt] &= 3\sin 2\theta\,(-0.5G_{02} + G_{22} \cos 2\phi + \tilde{G}_{22} \sin 2\phi) \nonumber\\ &\quad -3\cos 2\theta \,(G_{12} \cos\phi + \tilde{G}_{12} \sin\phi) + F_{11}\sin\phi - \tilde{F}_{11} \cos\phi \end{align}

and

(3.17c)\begin{align} \boldsymbol{\hat{\phi}} \boldsymbol{\cdot} [\kern-1pt[\,\boldsymbol{f}_H ]\kern-1pt] &= 3\cos \theta \, (G_{12} \sin \phi - \tilde{G}_{12} \cos \phi) - 6\sin \theta \,(G_{22} \sin 2\phi - \tilde{G}_{22} \cos 2\phi) ] \nonumber\\ &\quad + F_{01}\sin\theta + \cos\theta \, (F_{11} \cos\phi + \tilde{F}_{11}\sin\phi), \end{align}

where the stress coefficients are given by linear combinations of the flow coefficients:

(3.18a)\begin{gather} D_{ml} = -3 A_{ml} + 24B_{ml} + \tfrac{1}{7}a_{ml} - 4 \lambda b_{ml}, \end{gather}

(3.18b)\begin{gather}G_{ml} = \tfrac{1}{2} A_{ml} - 8B_{ml} - \tfrac{8}{21} a_{ml} - 2 \lambda b_{ml}, \end{gather}

(3.18c)\begin{gather}F_{ml} = -3 C_{ml}. \end{gather}

At this order of approximation, (3.13) for the velocity fields becomes

(3.19a,b)\begin{equation} \boldsymbol{v}^- = \boldsymbol{\nabla}\chi_1 \times \boldsymbol{r} + \frac{5 r^2 \boldsymbol{\nabla} p_2 - 4p_2 \boldsymbol{r}}{42\lambda} + \boldsymbol{\nabla}\varPhi_2,\quad \boldsymbol{v}^+ = \boldsymbol{\nabla}\chi_{-2} \times \boldsymbol{r} + \tfrac{1}{2}p_{-3} \boldsymbol{r} + \boldsymbol{\nabla}\varPhi_{-3} . \end{equation}

3.4. Kinematic boundary condition

The no-slip and no-penetration boundary conditions of (2.8) and (2.9) are

(3.20a–c)\begin{equation} v_\theta^+ = v_\theta^-, \quad v_\phi^+ = v_\phi^-,\quad v_r^+ - \delta \boldsymbol{v}^+ \boldsymbol{\cdot} \boldsymbol{\nabla} f = v_r^- - \delta \boldsymbol{v}^- \boldsymbol{\cdot} \boldsymbol{\nabla} f= \delta \dot{f}, \end{equation}

which provide relationships between the flow coefficients inside and outside the droplet. First, continuity of the polar and azimuthal components of the velocity yields the two relations

(3.21a,b)\begin{equation} \frac{5}{42 \lambda}a_{ml} + b_{ml} = B_{ml} \quad \mathrm{and} \quad c_{ml}=C_{ml}. \end{equation}

The no-penetration boundary condition (3.20c) has two terms that are of ${O}(\delta )$, which implies that $v_r^\pm$ must be of order ${O}(\delta )$ as well. Since $v_r^\pm \sim \lambda ^{-1}$, we require $\lambda ^{-1} = {O}(\delta )$, making the theory most accurate for viscous droplets that undergo small deformations due to high capillary forces. The radial components of the velocity inside and outside the drop are

(3.22a,b)\begin{equation} v_r^- = \frac{1}{7 \lambda} r p_2 + \frac{2}{r} \varPhi_2 \quad \mathrm{and} \quad v_r^+ = \frac{1}{2} rp_{-3} - \frac{3}{r}\varPhi_{-3}. \end{equation}

The no-penetration boundary condition in the radial direction gives us relations involving the temporal derivatives of the shape and the flow coefficients,

(3.23)\begin{equation} \frac{1}{7\lambda}a_{ml} + 2b_{ml} = \frac{1}{2} A_{ml} - 3B_{ml} = \delta (\dot{f}_{ml} + [\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]_{ml} ), \end{equation}

where we have introduced the following notation for the nonlinear product $\boldsymbol {v} \boldsymbol {\cdot } \boldsymbol {\nabla } f$:

(3.24)\begin{align} [\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f] &= \tfrac{1}{2}(3\cos^2\theta - 1)[\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]_{02} - 3\cos\theta\sin\theta\,([\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]_{12} \cos\phi \nonumber\\ &\quad + \widetilde{[\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]}_{12} \sin\phi) + 3\sin^2\theta \,([\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]_{22} \cos 2\phi + \widetilde{[\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]}_{22} \sin2\phi). \end{align}

The coefficients $[\boldsymbol {v} \boldsymbol {\cdot } \boldsymbol {\nabla } f]_{ml}$ are quadratic functions of $C_{ml}$ and $\,f_{ml}$ and are provided in appendix C. Note that $\boldsymbol {v} \boldsymbol {\cdot } \boldsymbol {\nabla } f$ only involves rigid-body rotation, as it is independent of viscosity, so that the leading-order term in (3.20c) is ${O}(\delta )$. Physically, it simply represents rotation of the droplet shape with the angular velocity $C$. We choose to express all the flow coefficients in terms of coefficients $B$ and $C$ (through $[\boldsymbol {v} \boldsymbol {\cdot } \boldsymbol {\nabla } f]$) as well as the shape transient $\dot {f}$:

(3.25a)\begin{gather} A_{ml} = 2 \{\delta (\,\dot{f}_{ml} + [\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]_{ml} ) + 3B_{ml} \}, \end{gather}

(3.25b)\begin{gather}a_{ml} = \tfrac{21}{2} \lambda \{ -\delta (\,\dot{f}_{ml} + [\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]_{ml} ) + 2B_{ml} \}, \end{gather}

(3.25c)\begin{gather}b_{ml} = \tfrac{1}{4} \{5\delta (\,\dot{f}_{ml} + [\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]_{ml} ) - 6B_{ml} \}. \end{gather}

Flows associated with coefficients $B$ and $C$ represent straining and rigid-body rotational flows, respectively.

3.5. Dynamic boundary condition

The dimensionless stress balance on the drop interface reads

(3.26)\begin{equation} [\kern-1pt[ p_E ]\kern-1pt] \boldsymbol{n} + q (E_\theta \boldsymbol{\hat{\theta}} + E_\phi \boldsymbol{\hat{\phi}} )+ Ma [\kern-1pt[\,\boldsymbol{f}_H ]\kern-1pt] = \frac{2}{Ca_E} \kappa_m \boldsymbol{n}. \end{equation}

We retain terms up to order ${O}(\delta )$ in this equation. In the radial direction, this reads

(3.27)\begin{equation} [\kern-1pt[ p_E ]\kern-1pt]+ Ma\, \boldsymbol{\hat{r}} \boldsymbol{\cdot} [\kern-1pt[\,\boldsymbol{f}_H ]\kern-1pt] = \frac{2}{Ca_E}(1+2\delta f). \end{equation}

Using orthogonality of spherical harmonics, we obtain one set of relations between the dipole moments, flow and shape coefficients:

(3.28)\begin{equation} \begin{aligned} & [\kern-1pt[ p_E ]\kern-1pt]_{00} - Ma ([\kern-1pt[ p_0 ]\kern-1pt]+ \tfrac{1}{2} D_{02} ) = 2Ca_E^{-1} (1 - \delta \,f_{02}),\quad [\kern-1pt[ p_E ]\kern-1pt]_{02} + \tfrac{3}{2} Ma D_{02} = 6 Ca_E^{-1} \delta \,f_{02}, \\ & [\kern-1pt[ p_E ]\kern-1pt]_{12} - \tfrac{3}{2} Ma D_{12} = -6 Ca_E^{-1} \delta \,f_{12},\quad [\kern-1pt[ \widetilde{p_E} ]\kern-1pt]_{12} - \tfrac{3}{2} Ma \tilde{D}_{12} = -6 Ca_E^{-1} \delta \tilde{f}_{12}, \\ & [\kern-1pt[ p_E ]\kern-1pt]_{22} + 3 Ma D_{22} = 12 Ca_E^{-1} \delta \,f_{22}, \quad [\kern-1pt[ \widetilde{p_E} ]\kern-1pt]_{22} + 3 Ma \tilde{D}_{22} = 12 Ca_E^{-1} \delta \tilde{f}_{22}, \end{aligned} \end{equation}

where the hydrodynamic radial stress coefficients $D$ are expressed in terms of the flow and shape coefficients using (3.18) and (3.25) as

(3.29)\begin{equation} D_{ml} = 3(2 + 3\lambda) B_{ml} - \tfrac{1}{2}\delta (12 + 13\lambda) (\,\dot{f}_{ml} + [\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]_{ml} ). \end{equation}

The ${O}(1)$ term on the right-hand side of the first equation in (3.28) is balanced by $[\kern-1pt[ p_0]\kern-1pt]$ and captures the Young–Laplace pressure jump due to surface tension in the undeformed state.

A second set of relations is obtained from the polar stress boundary condition

(3.30)\begin{equation} q E_\theta + Ma\, \boldsymbol{\hat{\theta}} \boldsymbol{\cdot} [\kern-1pt[\,\boldsymbol{f}_H ]\kern-1pt] = 0, \end{equation}

which, after applying orthogonality, yields

(3.31)\begin{equation} \begin{aligned} & \tilde{F}_{11} = \frac{1}{Ma} qE_{\theta ,10} , \quad F_{11} = - \frac{1}{Ma} \widetilde{q E}_{\theta ,10}, \\ & G_{02} = \frac{2}{3Ma} qE_{\theta ,02}, \quad G_{12} = \frac{1}{3Ma} qE_{\theta ,12} , \quad \tilde{G}_{12} = \frac{1}{3Ma} \widetilde{q E}_{\theta ,12}, \\ & G_{22} = - \frac{1}{3Ma} qE_{\theta ,22} , \quad \tilde{G}_{22} = - \frac{1}{3Ma} \widetilde{q E}_{\theta, 22} , \end{aligned} \end{equation}

where $G$ and $F$, similar to $D$, are expressed in terms of the flow and shape coefficients as

(3.32a,b)\begin{equation} F_{ml} = -3C_{ml}, \quad G_{ml} = -5(1+\lambda)B_{ml} + \tfrac{1}{2}\delta (2+3\lambda) (\,\dot{f}_{ml} + [\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]_{ml} ). \end{equation}

Balancing stress in the azimuthal direction,

(3.33)\begin{equation} q E_\phi + Ma\, \boldsymbol{\hat{\phi}} \boldsymbol{\cdot} [\kern-1pt[\,\boldsymbol{f}_H ]\kern-1pt] = 0, \end{equation}

gives us the same set of relations as in (3.31), with one additional relation,

(3.34)\begin{equation} F_{01} = \frac{1}{Ma} qE_{\phi ,01}. \end{equation}

The driving force for the fluid velocity are the tangential electric stresses, which induce hydrodynamic tractions that scale as ${O}(Ma^{-1})$. The magnitude of the resulting flow therefore is such that both electric and hydrodynamic radial tractions in (3.28) are of order ${O}(1)$. Balancing these tractions with surface tension forces requires us to choose $\delta \propto Ca^{-1}$. We choose to define the small-deformation parameter $\delta$ as

(3.35)\begin{equation} \delta=\frac{3Ca_E}{4(1+2R)^2}, \end{equation}

to be consistent with previous small-deformation theories (Taylor 1966). Physical quantities required to parametrize the drop shape, such as the extent of ellipsoidal deformations or tilt angle, are independent of this choice of $\delta$. For given values of the dipole moments $P_{x,y,z}$, the relations provided in this section, (3.28), (3.31) and (3.34), completely solve the EHD flow problem.

3.6. Charge conservation equation

We now proceed to derive a dipole moment evolution equation in each direction starting from the charge conservation equation. In dimensionless form, it reads

(3.36)\begin{equation} \frac{\partial q}{\partial t}+\left(\frac{Q+2}{1+2R}\right)\boldsymbol{n}\boldsymbol{\cdot} [\kern-1pt[\,\boldsymbol{j}]\kern-1pt] +\boldsymbol{\nabla}_{s}\boldsymbol{\cdot} (q\boldsymbol{v})=0, \end{equation}

where the surface charge density and jump in electric current have been scaled by $\epsilon E_0$ and $\sigma ^- E_0$, respectively. Taking moments of the charge conservation equation yields evolution equations for the charge coefficients:

(3.37a)\begin{gather} \dot{q}_{01} + \left(\frac{Q+2}{1+2R}\right) j_{n,01} + \frac{3}{4{\rm \pi}} \int_{0}^{\rm \pi} \int_0^{2{\rm \pi}} [\boldsymbol{\nabla}_{s}\boldsymbol{\cdot} (q\boldsymbol{v})] [\cos \theta] \sin\theta \,\mathrm{d}\phi \,\mathrm{d}\theta = 0, \end{gather}

(3.37b)\begin{gather}\dot{q}_{11} + \left(\frac{Q+2}{1+2R}\right) j_{n,11} + \frac{3}{4{\rm \pi}} \int_{0}^{\rm \pi} \int_0^{2{\rm \pi}} [\boldsymbol{\nabla}_{s}\boldsymbol{\cdot} (q\boldsymbol{v})] [\sin\theta\cos\phi] \sin\theta \,\mathrm{d}\phi \,\mathrm{d}\theta = 0, \end{gather}

(3.37c)\begin{gather}\dot{\tilde{q}}_{11} + \left(\frac{Q+2}{1+2R}\right) \tilde{j}_{n,11} + \frac{3}{4{\rm \pi}} \int_{0}^{\rm \pi} \int_0^{2{\rm \pi}} [\boldsymbol{\nabla}_{s}\boldsymbol{\cdot} (q\boldsymbol{v})] [\sin\theta\sin\phi] \sin\theta \,\mathrm{d}\phi \,\mathrm{d}\theta = 0. \end{gather}

In view of (3.8a–c) and (3.9a–c), these can also be regarded as the governing equations for the components of the dipole moment $\boldsymbol {P}$. The interfacial velocity is simply found from (3.19a,b) as $\boldsymbol {v}_s = \boldsymbol {v}^-|_{r=1} = \boldsymbol {v}^+|_{r=1}$. After substituting the expression for $\boldsymbol {v}_s$ along with kinematic boundary conditions and applying orthogonality of spherical harmonics, we obtain the final set of governing equations:

(3.38a)\begin{align} & \dot{q}_{01} + \left(\frac{Q+2}{1+2R}\right) j_{n,01}+ \delta \Big\{\frac{4}{5}\left(\,\dot{f}_{02} + [\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]_{02}\right) q_{01} \nonumber\\ &\quad -\frac{6}{5}\Big[\Big( \dot{\tilde{f}}_{12} + \widetilde{[\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]}_{12}\Big) q_{11} + \left(\,\dot{f}_{12} + [\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]_{12}\right)\tilde{q}_{11}\Big]\Big\} \nonumber\\ &\quad +\tilde{q}_{11}C_{11} - q_{11}\tilde{C}_{11} - \frac{6}{5} q_{01}B_{02} + \frac{9}{5}\Big(q_{11} B_{12} + \tilde{q}_{11}\tilde{B}_{12}\Big) = 0, \end{align}

(3.38b)\begin{align} & \dot{q}_{11} + \left(\frac{Q+2}{1+2R}\right) j_{n,11} + \delta \Big\{-\frac{2}{5}\left(\,\dot{f}_{02} + [\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]_{02}\right)q_{11} - \frac{6}{5} \left( \,\dot{f}_{12} + [\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]_{12}\right) q_{01} \nonumber\\ &\quad +\frac{12}{5} \Big[ \Big( \dot{f}_{22} + [\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]_{22}\Big) q_{11} + \Big( \dot{\tilde{f}}_{22} + \widetilde{[\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]}_{22} \Big) \tilde{q}_{11}\Big] \Big\} \nonumber\\ &\quad + q_{01} \Big(\tilde{C}_{11} + \frac{9}{5} B_{12}\Big) + \frac{3}{5} q_{11} B_{02} + \tilde{q}_{11} C_{01} - \frac{18}{5} \Big(q_{11} B_{22} + \tilde{q}_{11} \tilde{B}_{22}\Big) = 0 \end{align}

and

(3.38c)\begin{align} & \dot{\tilde{q}}_{11} + \left(\frac{Q+2}{1+2R}\right) \tilde{j}_{n,11} + \delta \Big\{-\frac{2}{5} \Big( \dot{f}_{02} + [\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]_{02} \Big) \tilde{q}_{11} -\frac{6}{5} \Big( \dot{\tilde{f}}_{12} + \widetilde{[\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]}_{12} \Big) q_{01} \nonumber\\ &\quad \frac{12}{5} \Big[\Big( \dot{f}_{22} + [\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]_{22} \Big) \tilde{q}_{11} - \Big( \dot{\tilde{f}}_{22} + \widetilde{[\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f]}_{22}\Big) q_{11}\Big] \Big\} \nonumber\\ &\quad -q_{01} \Big(C_{11} - \frac{9}{5} \tilde{B}_{12}\Big) + \frac{3}{5} \tilde{q}_{11} B_{02} - q_{11}C_{01} + \frac{18}{5} \Big(\tilde{q}_{11} B_{22} - q_{11} \tilde{B}_{22}\Big) = 0. \end{align}

The initial condition for the charge moments is zero charge, i.e. $q_{ml}(0) = 0$. Note that, in principle, the nonlinear term in (3.36) may excite higher-order charge multipoles beyond $q_{01}$, $q_{11}$ and $\tilde {q}_{11}$. These higher-order terms scale as ${O}(Ma^{-n} \lambda ^{-n})$ with $n>2$ and are neglected here, consistent with the truncation of the multipole expansion (3.4a,b) after the dipole term. Our theory therefore captures the leading-order ${O}(Ma^{-1} \lambda ^{-1})$ effect of straining flow in the charge convection term.

4. Recovering Quincke rotation of a solid sphere

We discuss the limit of $\lambda \rightarrow \infty$, i.e. Quincke rotation of a solid sphere, by considering only the rigid-body rotational flow in the charge convection term. The charge conservation equation in the high-viscosity limit reduces to

(4.1)\begin{equation} \begin{aligned} & \dot{q}_{01} + \left(\frac{Q+2}{1+2R}\right) j_{n,01} + \tilde{q}_{11}C_{11} - q_{11}\tilde{C}_{11} = 0, \\ & \dot{q}_{11} + \left(\frac{Q+2}{1+2R}\right) j_{n,11} + q_{01} \tilde{C}_{11} + \tilde{q}_{11} C_{01} = 0, \\ & \dot{\tilde{q}}_{11} + \left(\frac{Q+2}{1+2R}\right) \tilde{j}_{n,11} - q_{01} C_{11} - q_{11} C_{01} = 0, \end{aligned} \end{equation}

where the solid-body rotational flows are given by the torque balance equations,

(4.2a–c)\begin{equation} C_{01} = \frac{1}{3Ma} qE_{\phi ,01}, \quad C_{11} = - \frac{1}{3Ma} \widetilde{q E}_{\theta ,10}, \quad \tilde{C}_{11} = \frac{1}{3Ma} qE_{\theta ,10}. \end{equation}

In order to proceed, we choose to apply an electric field in the $y$-direction, which implies $P^+_y = 1+ 2P_y$, $P^-_y = 1 - P_y$, $P^+_{x,z} = 2P_{x,z}$ and $P^-_{x,z} = P_{x,z}$. This choice of field direction makes it easier for a direct comparison with the results of He et al. (2013). Substituting these relations in the polar and azimuthal stress components, we obtain an expression for the rotational flow fields in terms of the dipole moments:

(4.3a–c)\begin{equation} C_{01} = \frac{P_x}{2Ma}, \quad C_{11} = \frac{P_z}{2Ma}, \quad \tilde{C}_{11}= 0. \end{equation}

Substituting the dipole moments in the charge and jump in electric current provides coupled equations for the dipole components:

(4.4)\begin{equation} \begin{aligned} & \dot{P}_y + P_y - \frac{P^2_x}{2Ma}- \frac{P^2_z}{2Ma} = \frac{1-R}{1+2R}, \\ & \dot{P}_x + P_x \left[ 1 + \frac{1}{2Ma} \left( P_y -P_y^0 \right) \right] = 0, \\ & \dot{P}_z + P_z \left[ 1 + \frac{1}{2Ma} \left( P_y-P_y^0 \right) \right] = 0, \end{aligned} \end{equation}

where $P_y^0=P_y(0) = -(1-Q)/(Q+2)$ denotes the initial condition for the dipole moment in the absence of surface charge. With no loss of generality, we assume a dipole moment in the $(x,y)$ plane. At steady state, the system of equations (4.4) admits two solutions corresponding to no rotation, and steady Quincke rotation:

(4.5a,b)\begin{gather} P_x = 0, \quad P_y = \frac{1-R}{1+2R} \quad (\text{no rotation}), \end{gather}

(4.6a,b)\begin{gather}P_x = \sqrt{ 2Ma\left[ \frac{3(RQ-1)}{(Q+2)(1+2R)} - 2Ma \right] },\quad P_y = P_y^0-2Ma \quad (\text{Quincke rotation}). \end{gather}

The solid-body angular velocity is then easily found as

(4.7)\begin{equation} C_{01} = \sqrt{ \frac{3(RQ-1)}{2Ma(Q+2)(1+2R)} - 1 } = \sqrt{ \frac{E_0^2}{E_{c,s}^2} - 1 }, \end{equation}

where $E_{c,s}$ is the critical electric field for Quincke rotation of a solid sphere,

(4.8a–c)\begin{equation} E_{c,s}=\sqrt{\frac{2\mu^+}{\epsilon^+ \tau_{MW}(\bar{\epsilon}-\bar{\sigma})}} \quad \mathrm{and}\quad \bar{\epsilon}=\frac{\epsilon^- -\epsilon^+}{\epsilon^- + 2\epsilon^+}, \quad \bar{\sigma}=\frac{\sigma^- -\sigma^+}{\sigma^- + 2\sigma^+}. \end{equation}

Equations (4.4)–(4.8a–c) match classic results for Quincke electro-rotation of rigid spheres (Jones 1984; Das & Saintillan 2013) and highlight a supercritical pitchfork bifurcation in which the axisymmetric equilibrium state with no rotation becomes unstable and breaks symmetry for field strengths exceeding $E_{c,s}$. The system then evolves towards steady rotation around an arbitrary axis normal to the field and is characterized by a tilted dipole moment.

5. Results and discussion

We now compare our theoretical results with the existing experimental data of Salipante & Vlahovska (2010). Most previous studies have quantified departures from sphericity using Taylor's deformation parameter $\mathcal {D}_T=(r^\parallel -r^\perp )/(r^\parallel +r^\perp )$, where $r^\parallel$ and $r^\perp$ denote the lengths of the drop in the directions parallel and perpendicular to the electric field, respectively. The sign of $\mathcal {D}_T$ distinguishes between oblate ($\mathcal {D}_T<0$) and prolate ($\mathcal {D}_T>0$) shapes. However, its definition is ambiguous when the drop is tilted at an angle with respect to the electric field. Following He et al. (2013), we define a new parameter $\mathcal {D}_Q$ to quantify drop deformation,

(5.1)\begin{equation} \mathcal{D}_Q = \frac{L - B}{L + B}, \end{equation}

where $L$ and $B$ denote the lengths of the longest and shortest axes of the drop, respectively (see figure 1). We also introduce the tilt angle $\phi ^*$ as the angle between the longest axis of the drop and the plane normal to the applied field, such that $\phi ^*=0$ in the Taylor regime and $\phi ^*>0$ in the Quincke regime. It is worth noting that, unlike Taylor's deformation parameter, $D_Q$ is always positive. In the Taylor regime of axisymmetric deformations, their magnitudes are the same: $\mathcal {D}_Q = |\mathcal {D}_T|$. As the drop deformation and tilt angle are measured in the $(x,y)$ plane, the steady shape deformation simplifies to

(5.2)\begin{equation} f(\theta={\rm \pi}/2,\phi) = -\tfrac{1}{2}\,f_{02} + 3\,f_{22} \cos 2\phi + 3\tilde{f}_{22} \sin2\phi. \end{equation}

The tilt angle $\phi ^*$ then reads

(5.3)\begin{equation} \phi^* = \frac{1}{2} \tan^{-1}\left( \frac{ \tilde{f}_{22}}{\,f_{22}} \right). \end{equation}

The lengths of the longest and shortest axes of the drop are

(5.4a)\begin{gather} L(\phi=\phi^*) = 2 - \delta (\,f_{02} -6 \,f_{22} \cos 2\phi^* -6\tilde{f}_{22} \sin2\phi^* ), \end{gather}

(5.4b)\begin{gather}B(\phi={\rm \pi}/2 + \phi^*) = 2 - \delta (\,f_{02} +6 \,f_{22} \cos 2\phi^* +6\tilde{f}_{22} \sin2\phi^* ), \end{gather}

respectively, yielding the following expression for the deformation parameter in terms of the shape coefficients and tilt angle:

(5.5)\begin{equation} \mathcal{D}_Q = 3 \delta (\,f_{22} \cos 2\phi^* + \tilde{f}_{22} \sin 2\phi^* ) + {O}(\delta^2). \end{equation}

In the calculations presented below, we use the material properties, drop sizes and electric field strengths listed in table 1, which match the experimental values of Salipante & Vlahovska (2010). The dimensionless parameters $R=36.6$ and $Q=0.57$ are held fixed, while the other three dimensionless parameters $\lambda$, $Ca_E$ and $Ma$ change as we vary the drop viscosity $\mu ^-$, initial radius $a_d$ and electric field strength $E_0$.

Table 1. Material properties used herein corresponding to the experiments of Salipante & Vlahovska (2010). Here $\epsilon _0=8.8542 \times 10^{-12}\ \text {F}\ \textrm {m}^{-1}$ denotes the permittivity of vacuum.

5.1. Drop dynamics ignoring transient charge relaxation and straining flow

In this section, we neglect transient charge relaxation as well as the straining component of the flow in the charge conservation equation. In this case, our theory reduces to the model previously proposed by He et al. (2013) exactly. As $P_z =0$, we immediately deduce ${[\kern-1pt[ p_E ]\kern-1pt] _{12} = [\kern-1pt[ \widetilde {p_E} ]\kern-1pt] _{12} = qE_{\theta ,10} = \widetilde {q E}_{\theta ,10} = qE_{\theta ,12} = \widetilde {q E}_{\theta ,12}=0}$. Using the radial and polar stress balance equations, we find the evolution equations for the relevant shape functions:

(5.6a)\begin{gather} \delta \dot{f}_{02} = H_{02} = - \frac{1}{Ma \varLambda_1} ( \varLambda_2 \,f_{02} - 2[\kern-1pt[ p_E ]\kern-1pt]_{02} +2\varLambda_3qE_{\theta ,02} ), \end{gather}

(5.6b)\begin{gather}\delta (\,\dot{f}_{22} + 2C_{01}\,\tilde{f}_{22}) = H_{22} = - \frac{1}{Ma \varLambda_1} (\varLambda_2 \,f_{22} - [\kern-1pt[ p_E ]\kern-1pt]_{22} - \varLambda_3qE_{\theta ,22} ), \end{gather}

(5.6c)\begin{gather}\delta (\,\dot{\tilde{f}}_{22} - 2C_{01}\,f_{22}) = \tilde{H}_{22} = - \frac{1}{Ma \varLambda_1} (\varLambda_2 \,\tilde{f}_{22} - [\kern-1pt[ \widetilde{p_E} ]\kern-1pt]_{22} - \varLambda_3 \widetilde{q E}_{\theta,22} ), \end{gather}

where

(5.7a–c)\begin{equation} \varLambda_1=\frac{3 (19\lambda+16)(2\lambda+3)}{10(1+\lambda)},\quad \varLambda_2=\frac{9}{(1+2R)^2}, \quad \varLambda_3=\frac{3(2+3\lambda)}{5(1+\lambda)}. \end{equation}

From these expressions and substituting $C_{01}=P_x/2Ma$ from the azimuthal stress balance (3.34), we find the steady-state value of the shape coefficients $\,f_{22}$ and ${\tilde {f}_{22}}$ required for calculating $\phi ^*$ and $\mathcal {D}_Q$:

(5.8)\begin{gather} f_{22} = S [\varLambda_2 (\varLambda_3 qE_{\theta ,22} + [\kern-1pt[ p_E ]\kern-1pt]_{22} ) - \varLambda_1 \delta P_x (\varLambda_3 \widetilde{q E}_{\theta, 22} + [\kern-1pt[ \widetilde{p_E} ]\kern-1pt]_{22} ) ], \end{gather}

(5.9)\begin{gather}\,\tilde{f}_{22} = S [\varLambda_2 (\varLambda_3 \widetilde{q E}_{\theta, 22} + [\kern-1pt[ \widetilde{p_E} ]\kern-1pt]_{22} ) + \varLambda_1 \delta P_x (\varLambda_3 qE_{\theta ,22} + [\kern-1pt[ p_E ]\kern-1pt]_{22} ) ], \end{gather}

where $S= [\varLambda _2^2 + (\varLambda _1\delta P_x)^2]^{-1}$. Neglecting the straining part in the charge convection term results in the electric problem being the same as that for a solid spherical particle. Consequently, the critical electric field for the onset of Quincke rotation of a drop remains the same as that of a solid sphere. The dipole moments $P_{x,y}$ required in the above expressions are therefore simply given by the solutions obtained in (4.6a,b) for Quincke rotation of a solid sphere. The tilt angle and deformation for two drops of different sizes and viscosity ratios $(a_d,\lambda )=(0.9\ \textrm {mm},\ 14.1)$ and $(a_d,\lambda )=(0.7\ \textrm {mm},\ 1.41)$ are plotted in figure 2 (dashed lines). These results are discussed in more detail in the next section.

Figure 2. Effect of straining flows. (a) Tilt and (b) deformation parameter for two different datasets matching the experiments of Salipante & Vlahovska (2010). Red and blue colours correspond to drop sizes and viscosity ratios $(a_d,\lambda )=(0.9\ \textrm {mm},\ 14.1)$ and $(a_d,\lambda )=(0.7\ \textrm {mm},\ 1.41)$, respectively. Markers are experimental results, while solid and dashed lines are predictions from theory with both straining and rotational flows ($S+R$) and only rotational flow ($R$), respectively. Inclusion of straining flows is observed to increase the critical electric field for the onset of rotation, in good agreement with experimental results. (c) Flow assessment parameter $\zeta$ shows the relative magnitude of straining and rotational flows as a function of electric field strength for the same two drops in the Taylor and Quincke regimes. Movies showing the flow fields and drop shapes for these two cases are provided in the supplementary material accompanying this article.

5.2. Drop dynamics with transient charge relaxation and straining flow

Next, we look at the full model including the effects of transient charge relaxation and straining flow in the charge conservation equation. The tangential components of the interfacial velocity are of a similar form as the radial components in (5.6):

(5.10a)\begin{gather} B_{02} = - \frac{1}{Ma \varLambda_4} (\varLambda_2 \,f_{02} - 2[\kern-1pt[ p_E ]\kern-1pt]_{02} + 2\varLambda_5 qE_{\theta ,02} ), \end{gather}

(5.10b)\begin{gather}B_{22} = - \frac{1}{Ma \varLambda_4} (\varLambda_2 \,f_{22} - [\kern-1pt[ p_E ]\kern-1pt]_{22} - \varLambda_5 qE_{\theta ,22} ), \end{gather}

(5.10c)\begin{gather}\tilde{B}_{22} = - \frac{1}{Ma \varLambda_4} (\varLambda_2\, \tilde{f}_{22} - [\kern-1pt[ \widetilde{p_E} ]\kern-1pt]_{22} - \varLambda_5 \widetilde{qE}_{\theta ,22} ), \end{gather}

where

(5.11a,b)\begin{equation} \varLambda_4 = \frac{3(19\lambda+16)(2\lambda+3)}{2+3\lambda},\quad \varLambda_5 = \frac{12+13\lambda}{2+3\lambda}. \end{equation}

It is instructive to note that these straining flows scale as ${O}(Ma^{-1}\lambda ^{-1})$ and in principle induce higher-order multipoles neglected in this work. Using the charge conservation equation (3.38), we also obtain the dipole moment evolution equations for a droplet subject to an electric field in the $y$-direction,

(5.12a)\begin{align} & \dot{P}_x + P_x+ \left(-\frac{2}{5} H_{02} + \frac{12}{5} H_{22}+\frac{3}{5}B_{02}- \frac{18}{5}B_{22}\right)P_x \nonumber\\ &\quad +\left(\frac{12}{5} \tilde{H}_{22}+C_{01}-\frac{18}{5}\tilde{B}_{22}\right) (P_y -P_y^0 ) = 0, \end{align}

(5.12b)\begin{align} & \dot{P}_y + P_y + \left(\frac{12}{5} \tilde{H}_{22}-C_{01}-\frac{18}{5}\tilde{B}_{22} \right)P_x \nonumber\\ &\quad +\left(-\frac{2}{5} H_{02} - \frac{12}{5} H_{22}+\frac{3}{5}B_{02}+\frac{18}{5}B_{22}\right) (P_y -P_y^0 ) = \frac{1-R}{1+2R}. \end{align}

Equations (5.12) are subject to initial conditions $[P_x(0),P_y(0)]=(0,P_y^0)$ for the dipole, and we also impose ${\,f_{02}(0)=\,f_{22}(0)=\tilde {f}_{22}(0)=0}$ for an initially spherical drop. The solution to the full problem is then obtained by numerically integrating (5.6) and (5.12) in time using an explicit fourth-order Runge–Kutta scheme (Pozrikidis 1998) until steady-state values for the shape coefficients and dipole moments are reached. In the simulations, we typically introduce a small perturbation in the initial condition for $P_x$ at $t=0$, which either decays or amplifies depending on whether the drop is in the Taylor or Quincke regime, respectively. The MATLAB code used to integrate the equations has been made available to the reader and can be found in the supplementary material accompanying this article (available at https://doi.org/10.1017/jfm.2020.924).

In figure 2, we compare our theory with the experimental results of Salipante & Vlahovska (2010). Here, again, we consider two different systems with distinct drop sizes and viscosity ratios: $(a_d,\lambda )=(0.9\ \textrm {mm},\ 14.1)$ (in red) and $(a_d,\lambda )=(0.7\ \textrm {mm},\ 1.41)$ (in blue). The experiments clearly show that the critical electric field $E_{c,d}$ for the onset of Quincke rotation is higher for the less viscous drop. The complete theory including both straining and rotational flows (solid lines) in the charge convection equation accurately predicts the experimental results and captures the increase in the critical electric field. This effect is more pronounced for the less viscous drop, as the straining flow is stronger in that case, and indeed our model predicts $E_{c,d} = 1.081E_{c,s}$ for $\lambda =14.1$ and $E_{c,d} = 1.241E_{c,s}$ for $\lambda =1.41$ in terms of the critical field $E_{c,s}$ for a rigid sphere. On the other hand, the theory neglecting the straining flow (dashed lines) fails to predict any shift in the critical electric field. It is interesting to note that, as the electric field is increased into the Quincke regime, the rotational component of the flow becomes stronger and the difference between the predictions of the two theories therefore decreases. In figure 2(b), we plot the drop deformation $\mathcal {D}_Q$ as a function of electric field strength. In the Taylor regime, the drop keeps stretching and becomes more oblate with increasing field strength due to an increase in the straining flow. However, once the drop enters the Quincke regime, the rotational flow tends to compete with the straining flow, making the drop more spherical in shape. This effect seen in experiments is captured well by our theory.

In order to compare the relative strength of the straining and rotational components of the flow, we introduce a flow assessment parameter defined as

(5.13)\begin{equation} \zeta = \frac{\mathrm{tr}({\boldsymbol{\mathsf{S}}}^2) - \mathrm{tr}({\boldsymbol{\mathsf{W}}}^2)}{\mathrm{tr}({\boldsymbol{\mathsf{S}}}^2) + \mathrm{tr}({\boldsymbol{\mathsf{W}}}^2)}, \end{equation}

where ${\boldsymbol{\mathsf{S}}}$ and ${\boldsymbol{\mathsf{W}}}$ are the rate-of-strain and rate-of-rotation tensors, respectively (Das & Saintillan 2017a). With this definition, $\zeta =1$ represents purely straining flow while $\zeta =-1$ represents purely rotational flow. In order to determine the nature of the flow field inside the drop, we compute its value at the origin where it assumes a simple form in terms of the flow coefficients,

(5.14)\begin{equation} \zeta (\boldsymbol{r}=\boldsymbol{0})= \frac{9(b_{02}^2 + 8b_{22}^2 + 8\tilde{b}_{22}^2 ) + 2C_{01}^2}{9(b_{02}^2 + 8b_{22}^2 + 8\tilde{b}_{22}^2 ) - 2C_{01}^2}, \end{equation}

which we plot for the two different drops in figure 2(c). In the Taylor regime, $\zeta$ is identically 1, as rotational flows are absent, but it sharply drops in magnitude at the onset of Quincke rotation. The increase in critical electric field is seen in this plot as well if straining flow is included in the charge conservation equation (solid lines). In the Quincke regime, $\zeta \approx -1$ for the high-viscosity drop (red lines), as the rotational flow is much stronger than the straining flow. On the other hand, $\zeta$ assumes an intermediate value for the lower-viscosity drop (blue lines), as the straining flow is stronger in magnitude in that case. Movies showing the flow fields and drop shapes for these two cases are provided in the supplementary material.

Figure 3 compares the drop tilt angle predicted by the theory with the experiments of Salipante & Vlahovska (2010) for three different viscosity ratios, $\lambda =0.14, 1.41$ and 14.1, and a range of drop sizes, $a_d=0.50\text {--}2.45$ m. For a given viscosity ratio, our model predicts that the onset of Quincke rotation increases with increasing drop size, though this trend is not seen in experiments, as we further discuss below. Once Quincke rotation takes place, the smaller drops have a higher tilt angle compared to the larger ones, up to a certain electric field strength, after which this trend reverses. This is particularly visible for the lowest viscosity ratio of $\lambda =0.14$ in figure 3(c). The theory does best at capturing the onset of rotation in the case of high- and medium-viscosity drops in figures 3(a) and 3(b), respectively, but underpredicts the critical electric field for $\lambda =0.14$ in figure 3(c). This can be attributed to three reasons. Firstly, lower-viscosity drops undergo larger deformations so that the small-deformation theory becomes less accurate. Secondly, we recall from § 3.4 that the asymptotic order of our theory requires $\lambda ={O}(\delta ^{-1})$. Lastly, lower-viscosity drops have much higher straining flows, as the strength of the straining flow scales as ${\sim }\lambda ^{-1}$ for fixed Mason number $Ma$. These straining flows result in stronger charge convection, which can only be captured by including higher-order terms in ${O}(Ma^{-1}\lambda ^{-1})$ (see relevant discussion in § 3.6). However, once rotation has ensued, drop deformation decreases and straining flow is overcome by rotational flow, and, as a consequence, the drop tilt angle is very well captured by the theory for sufficiently large $E/E_{c,s}$ past the onset of rotation, regardless of the viscosity ratio. We also note that, while the theory developed here is valid in the weak charge convection regime, the results are found to match well with the experiments even when $Ma^{-1} \lambda ^{-1} \sim 1$, i.e. for the high- and medium-viscosity drops.

Figure 3. Tilt angle for drops of various sizes $a_d$ and viscosity ratios (a) $\lambda = 14.1$, (b) $\lambda = 1.41$ and (c) $\lambda = 0.14$. Lines are predictions from theory, which includes both straining and rotational components of the flow, as presented in § 5.2. Squares, diamonds, triangles and circles are experimental results of Salipante & Vlahovska (2010). The critical electric field for the onset of Quincke rotation increases with a decrease in viscosity ratio and is well captured by the theory.

5.3. Linear stability analysis for critical electric field

In this section, we perform a linear stability analysis on the governing equations (5.6) and (5.12) to determine the critical electric field for the onset of electro-rotation in the case of a drop. The equilibrium state, denoted with a superscript $e$, is the axisymmetric Taylor regime with no rotation and a dipole moment perfectly aligned with the field. For an applied field in the $y$-direction, we obtain $P_y^e$ by solving the cubic equation

(5.15)\begin{equation} P_{y}^e + \frac{2(Q+2)}{25Ma(1+\lambda)} (P_{y}^e-1 ) (P_{y}^e -P_y^0 )^2 = \frac{1-R}{1+2R}. \end{equation}

In the Taylor regime, the shape is axisymmetric. It is captured by coefficients $f^e_{02}$ and $f^e_{22}$, which can be found by using (5.6). We perturb the equilibrium base state as

(5.16)\begin{equation} \begin{pmatrix} P_x \\ P_y \\ f_{02} \\ f_{22} \\ \,\tilde{f}_{22} \end{pmatrix} = \begin{pmatrix} 0 \\ P^e_y \\ \,f_{02}^e \\ f_{22}^e \\ 0 \end{pmatrix}+\alpha \begin{pmatrix} P_x^\prime \\ P_y^\prime \\ f_{02}^\prime \\ f_{22}^\prime \\ \,\tilde{f}_{22}^\prime \end{pmatrix}, \end{equation}

where primed variables denote perturbations. In the limit of $\alpha \ll 1$, we linearize (5.6) and (5.12) and obtain an eigenvalue problem,

(5.17)\begin{equation} \frac{\mathrm{d}}{\mathrm{d} t} \begin{pmatrix} P_x^\prime \\ P_y^\prime \\ f_{02}^\prime \\ f_{22}^\prime \\ \,\tilde{f}_{22}^\prime \end{pmatrix} = \begin{pmatrix} \mathcal{J}_{11} & 0 & 0 & 0 & \mathcal{J}_{15}\\ 0 & \mathcal{J}_{22} & \mathcal{J}_{23} & \mathcal{J}_{24} & 0 \\ 0 & \mathcal{J}_{32} & \mathcal{J}_{33} & 0 & 0 \\ 0 & \mathcal{J}_{42} & 0 & \mathcal{J}_{44} & 0 \\ \mathcal{J}_{51} & 0 & 0 & 0 & \mathcal{J}_{55}\\ \end{pmatrix} \begin{pmatrix} P_x^\prime \\ P_y^\prime \\ f_{02}^\prime \\ f_{22}^\prime \\ \,\tilde{f}_{22}^\prime \end{pmatrix}, \end{equation}

where the exact expressions for the components of the constant Jacobian matrix $\mathcal {J}=\mathcal {J}(Ma,\lambda ,R,Q)$ are provided in appendix D. Inspection of (5.17) shows that the transverse dipole moment and shape perturbation $P_x^\prime$ and $\smash {\tilde {f}^\prime _{22}}$ depend on each other but are uncoupled from the other variables, a consequence of the quadratic nature of electric stresses. To analyse the stability, we therefore need only to consider the simpler system

(5.18)\begin{equation} \frac{\mathrm{d}}{\mathrm{d} t} \begin{pmatrix} P_x^\prime \\ \,\tilde{f}_{22}^\prime \end{pmatrix} = \begin{pmatrix} \mathcal{J}_{11} & \mathcal{J}_{15}\\ \mathcal{J}_{51} & \mathcal{J}_{55}\\ \end{pmatrix} \begin{pmatrix} P_x^\prime \\ \,\tilde{f}_{22}^\prime \end{pmatrix}. \end{equation}

Its eigenvalues are then determined numerically, and the critical value $E_c/E_{c,s}$ of the applied field at which at least one eigenvalue becomes positive is recorded.

In figure 4, we plot the variation of the critical electric field $E_c$ for the onset of Quincke rotation of drops for a wide range of drop sizes and viscosity ratios, comparing experiments with numerical solutions. Figure 4(a) shows the variation of the critical electric field with drop radius for various viscosity ratios $\lambda =0.14\text {--}14.1$. The critical field for a given drop size is observed to decrease with increase in viscosity ratio. However, across the range of viscosity ratios, it is observed in experiments that bigger drops have a lower critical electric field than smaller ones. This is shown more clearly in figure 4(b), where we plot the critical electric field as a function of viscosity ratio for various drop sizes. The observation that smaller drops have a higher critical electric field for Quincke rotation has also been reported in boundary element simulations (Das & Saintillan 2017a). However, it is not captured by the first-order small-deformation theory, which predicts the opposite trend. We speculate that higher-order contributions of the dipole moment may be responsible for the discrepancy. It has already been noted in previous studies that the critical electric field of an oblate spheroid (bigger and more deformed drop) is higher than that of a solid sphere (smaller and less deformed drop) due to a stronger induced dipole in the former case (Cēbers, Lemaire & Lobry 2000; Salipante & Vlahovska 2013). We also note that the hydrodynamic torque required for rotating an oblate spheroid is smaller compared to that of a sphere, which acts to decrease the critical electric field for the onset of rotation (Kim & Karrila 2013). A full analysis of these effects is beyond the scope of this paper. It is also noted inexplicably that, while the trend of $E_c$ as a function of the drop radius $a_d$ disagrees with experiments, its variation with the viscosity ratio $\lambda$ tends to agree well with experiments. Predictions of the linear stability analysis for the critical electric field give the same result as that obtained by numerically integrating the complete set of ordinary differential equations (ODEs) (5.6) and (5.12).

Figure 4. Critical electric field versus (a) drop radius and (b) viscosity ratio. The critical field $E_c$ for a drop is normalized by the critical value $E_{c,s}$ for a rigid spherical particle. Lines are predictions from the theory including straining flow, while markers are experimental results of Salipante & Vlahovska (2010). The critical field obtained from the linear stability analysis gives the same result as the numerical solutions; hence, it is not shown here for brevity.

6. Conclusions

We have developed a leading-order small-deformation theory for the EHD of a dielectric drop in an applied electric field using the complete Melcher–Taylor leaky-dielectric model in three dimensions. Our model improves upon past theories that had neglected the effect of straining flows or been limited to axisymmetric deformations. It is most accurate in the limits of strong capillary forces, highly viscous drops and weak but finite charge convection, and is able to capture the transition to Quincke rotation seen in experiments and simulations. A coupled set of nonlinear ODEs involving dipole moments, shape and flow field coefficients were derived. Numerical solutions of this systems of ODEs were compared with existing experiments and showed excellent agreement in the appropriate limits of high viscosity and small deformations.

The main novelty of our work is the retention of the nonlinear straining flow in the governing equation for charge transport, which allowed us to explore variations in the critical electric field for electro-rotation observed in both experiments (Salipante & Vlahovska 2010) and simulations (Das & Saintillan 2017a). Predictions of the theory for the drop deformation and tilt angle, and for the critical electric field, showed good quantitative agreement with experimental data in the appropriate limits, with the exception of the dependence of $E_c$ on drop size, which we hypothesize requires higher-order effects to be captured accurately. This may perhaps be achieved using a spheroidal drop theory allowing for large deformations as in the work of Zhang, Zahn & Lin (2013). While drop deformation increases monotonically with electric field strength in the Taylor regime due to the sole presence of straining flows, our theory shows that, once Quincke rotation occurs, the rotational flow dominates and tends to decrease deformations in agreement with experimental observations. Our theoretical model, in addition to furthering our fundamental understanding of the role of interfacial flows in droplet EHD, may also be useful for building models for emulsions of interacting drops (Varshney et al. 2016; Zabarankin 2020) and for understanding other instabilities and complex dynamics observed in related systems (Dommersnes et al. 2013; Ouriemi & Vlahovska 2014; Brosseau & Vlahovska 2017).