3.1 Shape parametrization and expansion

In axisymmetric geometry, we parametrize the drop shape as a curve $\unicode[STIX]{x1D709}(r,\unicode[STIX]{x1D702})=0$ , where $r=|\boldsymbol{x}|$ is the distance from the drop centre and $\unicode[STIX]{x1D702}=\cos \unicode[STIX]{x1D703}$ is the cosine of the polar angle $\unicode[STIX]{x1D703}\in [0,\unicode[STIX]{x03C0}]$ measured from the field direction. For small deviations from sphericity, the drop shape is expanded on the basis of spherical harmonics as

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D709}(r,\unicode[STIX]{x1D702})=r-(1+\unicode[STIX]{x1D6FF}f_{1}+\unicode[STIX]{x1D6FF}^{2}f_{2})+O(\unicode[STIX]{x1D6FF}^{3}).\end{eqnarray}$$

The first- and second-order shape functions are linear combinations of Legendre polynomials ${\mathcal{L}}_{n}$ of order $n$ :

(3.2a ) $$\begin{eqnarray}\displaystyle f_{1} & = & \displaystyle f_{12}{\mathcal{L}}_{2}(\unicode[STIX]{x1D702}),\end{eqnarray}$$

(3.2b ) $$\begin{eqnarray}\displaystyle f_{2} & = & \displaystyle f_{20}+f_{22}{\mathcal{L}}_{2}(\unicode[STIX]{x1D702})+f_{24}{\mathcal{L}}_{4}(\unicode[STIX]{x1D702}),\end{eqnarray}$$

${\mathcal{L}}_{2}$

${\mathcal{L}}_{4}$

$b$

(3.3) $$\begin{eqnarray}\displaystyle \int _{0}^{\unicode[STIX]{x03C0}}{\mathcal{L}}_{i}(\unicode[STIX]{x1D702}){\mathcal{L}}_{j}(\unicode[STIX]{x1D702})\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}=\frac{2}{2i+1}\unicode[STIX]{x1D6FF}_{ij}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}_{ij}$ denotes the Kronecker delta, which will become useful later. The choice of Legendre functions in (3.2a )–(3.2b ) is a consequence of the quadratic nature of the Maxwell electric stresses acting on the fluid–drop interface, which in a uniform electric field only excite shape modes of order $2^{n}$ ( $n\in \mathbb{Z}^{+}$ ). In (3.2a )–(3.2b ) and in the rest of the paper, pairs of indices in coefficients of the form $f_{ij}$ refer to the order $i$ in the small-deformation expansion and to the order $j$ of the Legendre polynomial they multiply, respectively. In (3.2b ), the constant term $f_{20}$ is added to the second-order shape function $f_{2}$ to negate the perturbation in the drop volume due to the first-order shape function $f_{1}$ :

(3.4) $$\begin{eqnarray}2\unicode[STIX]{x03C0}\int _{0}^{\unicode[STIX]{x03C0}}\int _{0}^{r}\unicode[STIX]{x1D70C}^{2}\sin \unicode[STIX]{x1D703}\hspace{2.22198pt}\text{d}\unicode[STIX]{x1D70C}\hspace{2.22198pt}\text{d}\unicode[STIX]{x1D703}=\frac{4\unicode[STIX]{x03C0}}{3}+4\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}^{2}\left(\frac{f_{12}^{2}}{5}+f_{20}\right)+O(\unicode[STIX]{x1D6FF}^{3}).\end{eqnarray}$$

Requiring terms of order $\unicode[STIX]{x1D6FF}^{2}$ to vanish, we get $f_{20}=-f_{12}^{2}/5$ . The outward unit normal, tangent vector and curvature of the interface are also obtained as (Ajayi Reference Ajayi1978)

(3.5a,b ) $$\begin{eqnarray}\displaystyle & \boldsymbol{n}=\hat{\boldsymbol{e}}_{r}-\unicode[STIX]{x1D6FF}\,\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}f_{1}\,\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D703}}+O(\unicode[STIX]{x1D6FF}^{2}),\quad \boldsymbol{t}=\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D703}}+\unicode[STIX]{x1D6FF}\,\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}f_{1}\,\hat{\boldsymbol{e}}_{r}+O(\unicode[STIX]{x1D6FF}^{2}), & \displaystyle\end{eqnarray}$$

(3.5c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }\boldsymbol{n}=2-\unicode[STIX]{x1D6FF}L[f_{1}]-\unicode[STIX]{x1D6FF}^{2}\{L[f_{2}]-2f_{1}(L[f_{1}]-f_{1})\}+O(\unicode[STIX]{x1D6FF}^{3}), & & \displaystyle\end{eqnarray}$$

$L$

$L[f]=\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}\{(1-\unicode[STIX]{x1D702}^{2})\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}f\}+2f$

Using the above parametrization, the normal and tangential components of any vector $\boldsymbol{v}$ and second-order tensor $\unicode[STIX]{x1D64F}$ on the drop surface are related to their components in spherical coordinates by

(3.6a ) $$\begin{eqnarray}\displaystyle v^{n} & = & \displaystyle v_{0}^{r}+\unicode[STIX]{x1D6FF}(v_{1}^{r}+f_{1}\unicode[STIX]{x2202}_{r}v_{0}^{r}-\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}f_{1}v_{0}^{\unicode[STIX]{x1D703}})+O(\unicode[STIX]{x1D6FF}^{2}),\end{eqnarray}$$

(3.6b ) $$\begin{eqnarray}\displaystyle v^{t} & = & \displaystyle v_{0}^{\unicode[STIX]{x1D703}}+\unicode[STIX]{x1D6FF}(v_{1}^{\unicode[STIX]{x1D703}}+f_{1}\unicode[STIX]{x2202}_{r}v_{0}^{\unicode[STIX]{x1D703}}+\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}f_{1}v_{0}^{r})+O(\unicode[STIX]{x1D6FF}^{2}),\end{eqnarray}$$

(3.6c ) $$\begin{eqnarray}\displaystyle T^{nn} & = & \displaystyle T_{0}^{rr}+\unicode[STIX]{x1D6FF}(T_{1}^{rr}+f_{1}\unicode[STIX]{x2202}_{r}T_{0}^{rr}-2\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}f_{1}T_{0}^{r\unicode[STIX]{x1D703}})+O(\unicode[STIX]{x1D6FF}^{2}),\end{eqnarray}$$

(3.6d ) $$\begin{eqnarray}\displaystyle T^{nt} & = & \displaystyle T_{0}^{r\unicode[STIX]{x1D703}}+\unicode[STIX]{x1D6FF}[T_{1}^{r\unicode[STIX]{x1D703}}+f_{1}\unicode[STIX]{x2202}_{r}T_{0}^{r\unicode[STIX]{x1D703}}+\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}f_{1}(T_{0}^{rr}-T_{0}^{\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}})]+O(\unicode[STIX]{x1D6FF}^{2}),\end{eqnarray}$$

$r=1$

3.2 Electric problem

3.2.1 Spherical harmonic expansion

We first present the solution to the electric problem, which consists in solving (2.1)–(2.3) asymptotically. The electric potential outside and inside the drop can be expanded in powers of $\unicode[STIX]{x1D6FF}$ as

(3.7a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D711} & = & \displaystyle \unicode[STIX]{x1D711}_{e}(r,\unicode[STIX]{x1D703})+\unicode[STIX]{x1D711}_{0}(r,\unicode[STIX]{x1D703})+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}_{1}(r,\unicode[STIX]{x1D703})+O(\unicode[STIX]{x1D6FF}^{2}),\end{eqnarray}$$

(3.7b ) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D711}} & = & \displaystyle \unicode[STIX]{x1D711}_{e}(r,\unicode[STIX]{x1D703})+\bar{\unicode[STIX]{x1D711}}_{0}(r,\unicode[STIX]{x1D703})+\unicode[STIX]{x1D6FF}\bar{\unicode[STIX]{x1D711}}_{1}(r,\unicode[STIX]{x1D703})+O(\unicode[STIX]{x1D6FF}^{2}),\end{eqnarray}$$

which automatically satisfies the far-field boundary condition (2.3). We have yet to enforce continuity of the potential across the interface. To this end, we employ a domain perturbation approach in which all the boundary conditions are enforced approximately on the undeformed spherical surface $r=1$ . The potential on the interface is first expanded in the neighbourhood of $r=1$ using Taylor series:

(3.8a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D711} & = & \displaystyle \unicode[STIX]{x1D711}_{e}+\unicode[STIX]{x1D711}_{0}+\unicode[STIX]{x1D6FF}\left[\unicode[STIX]{x1D711}_{1}+f_{1}\unicode[STIX]{x2202}_{r}(\unicode[STIX]{x1D711}_{e}+\unicode[STIX]{x1D711}_{0})\right]+O(\unicode[STIX]{x1D6FF}^{2}),\end{eqnarray}$$

(3.8b ) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D711}} & = & \displaystyle \bar{\unicode[STIX]{x1D711}}_{e}+\bar{\unicode[STIX]{x1D711}}_{0}+\unicode[STIX]{x1D6FF}\left[\bar{\unicode[STIX]{x1D711}}_{1}+f_{1}\unicode[STIX]{x2202}_{r}(\unicode[STIX]{x1D711}_{e}+\bar{\unicode[STIX]{x1D711}}_{0})\right]+O(\unicode[STIX]{x1D6FF}^{2}).\end{eqnarray}$$

Applying continuity (2.2) and matching terms of zeroth and first order in $\unicode[STIX]{x1D6FF}$ provides two boundary conditions at $r=1$ :

(3.9a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D711}_{0} & = & \displaystyle \bar{\unicode[STIX]{x1D711}}_{0},\end{eqnarray}$$

(3.9b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D711}_{1}+f_{1}\unicode[STIX]{x2202}_{r}(\unicode[STIX]{x1D711}_{e}+\unicode[STIX]{x1D711}_{0}) & = & \displaystyle \bar{\unicode[STIX]{x1D711}}_{1}+f_{1}\unicode[STIX]{x2202}_{r}(\unicode[STIX]{x1D711}_{e}+\bar{\unicode[STIX]{x1D711}}_{0}).\end{eqnarray}$$

The zeroth-order problem, which is identical to the case of a sphere, is easily solved using decaying and growing spherical harmonics in terms of electric dipoles $P_{01}$ , $\bar{P}_{01}$ :

(3.10a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D711}_{0} & = & \displaystyle P_{01}r^{-2}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702}),\end{eqnarray}$$

(3.10b ) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D711}}_{0} & = & \displaystyle \bar{P}_{01}r{\mathcal{L}}_{1}(\unicode[STIX]{x1D702}),\end{eqnarray}$$

and we require that $\bar{P}_{01}=P_{01}$ to satisfy (3.9a ); solving for $P_{01}$ will require application of the charge conservation (2.5) as detailed below. After substitution into (3.9b ), we obtain a new first-order boundary condition:

(3.11) $$\begin{eqnarray}\unicode[STIX]{x1D711}_{1}-\bar{\unicode[STIX]{x1D711}}_{1}=3f_{1}P_{01}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})=3f_{12}P_{01}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702}){\mathcal{L}}_{2}(\unicode[STIX]{x1D702})={\textstyle \frac{3}{5}}f_{12}P_{01}[2{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})+3{\mathcal{L}}_{3}(\unicode[STIX]{x1D702})].\end{eqnarray}$$

The order of the polynomials appearing on the right-hand side suggests representing the first-order potentials in terms of both dipoles $P_{11}$ , $\bar{P}_{11}$ and octupoles $P_{13}$ , $\bar{P}_{13}$ :

(3.12a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D711}_{1} & = & \displaystyle P_{11}r^{-2}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})+P_{13}r^{-4}{\mathcal{L}}_{3}(\unicode[STIX]{x1D702}),\end{eqnarray}$$

(3.12b ) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D711}}_{1} & = & \displaystyle \bar{P}_{11}r{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})+\bar{P}_{13}r^{3}{\mathcal{L}}_{3}(\unicode[STIX]{x1D702}),\end{eqnarray}$$

and application of the boundary condition (3.11) yields the relations

(3.13a,b ) $$\begin{eqnarray}\bar{P}_{11}=P_{11}-{\textstyle \frac{6}{5}}f_{12}P_{01},\quad \bar{P}_{13}=P_{13}-{\textstyle \frac{9}{5}}f_{12}P_{01}.\end{eqnarray}$$

Having determined the electric potential, we can also obtain asymptotic expressions for the normal and tangential electric fields $E^{n}=-\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}$ and $E^{t}=-\boldsymbol{t}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}$ on the drop surface. Applying (3.6a ), we find

(3.14) $$\begin{eqnarray}E^{n}=E_{0}^{n}+\unicode[STIX]{x1D6FF}E_{1}^{n}+O(\unicode[STIX]{x1D6FF}^{2})=E_{01}^{n}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D6FF}[E_{11}^{n}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})+E_{13}^{n}{\mathcal{L}}_{3}(\unicode[STIX]{x1D702})]+O(\unicode[STIX]{x1D6FF}^{2}),\end{eqnarray}$$

with a similar expansion for $\bar{E}^{n}$ . Finally, the expansion for the tangential electric field, which is continuous across the interface, is obtained using (3.6b ) as

(3.15) $$\begin{eqnarray}E^{t}=E_{0}^{t}+\unicode[STIX]{x1D6FF}E_{1}^{t}+O(\unicode[STIX]{x1D6FF}^{2})=E_{00}^{t}\sin \unicode[STIX]{x1D703}+\unicode[STIX]{x1D6FF}[E_{10}^{t}+E_{12}^{t}{\mathcal{L}}_{2}(\unicode[STIX]{x1D702})]\sin \unicode[STIX]{x1D703}+O(\unicode[STIX]{x1D6FF}^{2}).\end{eqnarray}$$

The detailed expressions for the coefficients appearing in (3.14)–(3.15) are provided in appendix A.

3.2.2 Charge conservation and moment equations

To complete the solution of the electric problem, equations must be derived for the moments $P_{01}$ , $P_{11}$ and $P_{13}$ , which are time-dependent. These can be obtained as ordinary differential equations by application of the charge conservation (2.5). First, we expand the charge density in powers of $\unicode[STIX]{x1D6FF}$ as

(3.16) $$\begin{eqnarray}q=q_{0}+\unicode[STIX]{x1D6FF}q_{1}+O(\unicode[STIX]{x1D6FF}^{2})=q_{01}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D6FF}[q_{11}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})+q_{13}{\mathcal{L}}_{3}(\unicode[STIX]{x1D702})]+O(\unicode[STIX]{x1D6FF}^{2}).\end{eqnarray}$$

Similarly, we expand the jump in Ohmic currents $\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x27E6}\,\boldsymbol{j}\unicode[STIX]{x27E7}=\unicode[STIX]{x27E6}\,j\unicode[STIX]{x27E7}^{n}$ , scaled here by $\bar{\unicode[STIX]{x1D70E}}E_{0}$ , as

(3.17) $$\begin{eqnarray}\unicode[STIX]{x27E6}\,j\unicode[STIX]{x27E7}^{n}=\unicode[STIX]{x27E6}\,j\unicode[STIX]{x27E7}_{0}^{n}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x27E6}\,j\unicode[STIX]{x27E7}_{1}^{n}+O(\unicode[STIX]{x1D6FF}^{2})=\unicode[STIX]{x27E6}\,j\unicode[STIX]{x27E7}_{01}^{n}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D6FF}\{\unicode[STIX]{x27E6}\,j\unicode[STIX]{x27E7}_{11}^{n}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})+\unicode[STIX]{x27E6}\,j\unicode[STIX]{x27E7}_{13}^{n}{\mathcal{L}}_{3}(\unicode[STIX]{x1D702})\}+O(\unicode[STIX]{x1D6FF}^{2}).\end{eqnarray}$$

Expansion coefficients for both $q$ and $\unicode[STIX]{x27E6}\,j\unicode[STIX]{x27E7}^{n}$ are provided in appendix A. Finally, we formally expand the charge convection term in (2.5) as

(3.18) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v}) & = & \displaystyle [\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v})]_{0}+\unicode[STIX]{x1D6FF}[\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v})]_{1}+O(\unicode[STIX]{x1D6FF}^{2}),\nonumber\\ \displaystyle & = & \displaystyle [\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v})]_{01}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D6FF}\{[\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v})]_{11}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})+[\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v})]_{13}{\mathcal{L}}_{3}(\unicode[STIX]{x1D702})\}+O(\unicode[STIX]{x1D6FF}^{2}),\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where we have introduced the Legendre coefficients

(3.19) $$\begin{eqnarray}\displaystyle [\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v})]_{ij}=\frac{2j+1}{2}\int _{0}^{\unicode[STIX]{x03C0}}[\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v})]_{i}{\mathcal{L}}_{j}(\unicode[STIX]{x1D702})\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}. & & \displaystyle\end{eqnarray}$$

Detailed expressions for these coefficients require knowledge of the interfacial velocity $\boldsymbol{v}$ , whose calculation is presented in § 3.3. Note that the zeroth-order charge convection term arising from the nonlinear product of the charge density $q_{0}$ and interfacial velocity $\boldsymbol{v}_{0}$ also involves an additional term of the form $[\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v})]_{03}{\mathcal{L}}_{3}(\unicode[STIX]{x1D702})$ , which has been neglected here. To consider its effect, one would need to include a zeroth-order octupole $P_{03}$ in (3.10), which in turn would generate charge convection terms of orders $1,3,5,7,\ldots ,$ thus requiring additional higher-order odd multipoles. These multipoles becomes stronger with increasing electric Reynolds number or decreasing Mason number. Our theory is therefore valid in the limit of high Mason number, i.e. for drop–fluid systems in which charge convection is relatively weak.

Substituting the expansions (3.16)–(3.18) into the charge conservation (2.5), matching powers of $\unicode[STIX]{x1D6FF}$ , and applying orthogonality of Legendre polynomials leads to a set of relaxation equations for the charge coefficients. In dimensionless form, these read

(3.20) $$\begin{eqnarray}\dot{q}_{ij}+\frac{Q+2}{1+2R}\,\unicode[STIX]{x27E6}\,j\unicode[STIX]{x27E7}_{ij}^{n}+[\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v})]_{ij}=0,\end{eqnarray}$$

where the dot in the first term denotes differentiation with respect to time. If we further express $q_{ij}$ and $\unicode[STIX]{x27E6}\,j\unicode[STIX]{x27E7}_{ij}^{n}$ in terms of $P_{01}$ , $P_{11}$ and $P_{13}$ using (A 2) and (A 3), we arrive at a set of hierarchical differential equations for the dipole and octupole moments:

(3.21) $$\begin{eqnarray}\displaystyle {\dot{P}}_{01}+P_{01}=\frac{1-R}{1+2R}-\frac{1}{Q+2}[\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v})]_{01}, & & \displaystyle\end{eqnarray}$$

(3.22) $$\begin{eqnarray}\displaystyle {\dot{P}}_{11}+P_{11} & = & \displaystyle \frac{\text{d}}{\text{d}t}\left[\frac{6}{5}f_{12}\left(P_{01}\frac{1+2Q}{2+Q}+\frac{1-Q}{2+Q}\right)\right]+\frac{6}{5}f_{12}\left(P_{01}\frac{R+2}{2R+1}-\frac{1-R}{2R+1}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\frac{1}{Q+2}[\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v})]_{11},\end{eqnarray}$$

(3.23) $$\begin{eqnarray}\displaystyle {\dot{P}}_{13}+\frac{Q+2}{3Q+4}\frac{4R+3}{2R+1}P_{13} & = & \displaystyle \frac{\text{d}}{\text{d}t}\left[\frac{6}{5}f_{12}\left(P_{01}\frac{8+7Q}{8+6Q}-\frac{1-Q}{4+3Q}\right)\right]\nonumber\\ \displaystyle & & \displaystyle +\,\frac{6}{5}f_{12}\frac{Q+2}{3Q+4}\left(P_{01}\frac{8R+7}{4R+2}+\frac{1-R}{2R+1}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\frac{1}{3Q+4}[\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v})]_{13}.\end{eqnarray}$$

These coupled ordinary differential equations constitute the main result of this section. The external forcing in these equations is encapsulated in the first term on the right-hand side of (3.21), which describes the effect of the applied electric field on the leading-order dipole moment. It is also clear from (3.21) that charge convection cannot be neglected even at zeroth order, as was previously done in the theories of Shutov (Reference Shutov2002) and Shkadov & Shutov (Reference Shkadov and Shutov2002). Solving (3.21)–(3.23) requires the Legendre coefficients of the charge convection term as well as the first-order shape coefficient $f_{12}$ . These unknowns will be determined below after we solve for the fluid flow, which affects both interfacial charge convection and droplet deformation.

3.3 Flow problem: streamfunction formulation

We now turn to the calculation of the fluid flow outside and inside the drop. Upon application of the field, electric stresses develop at the interface, leading to deformations and flow. Since the flow is axisymmetric, we use a Stokes streamfunction $\unicode[STIX]{x1D6F9}(r,\unicode[STIX]{x1D703})$ to determine the fluid velocity, which has components

(3.24a,b ) $$\begin{eqnarray}v^{r}=\frac{1}{r^{2}\sin \unicode[STIX]{x1D703}}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}\unicode[STIX]{x1D6F9},\quad v^{\unicode[STIX]{x1D703}}=-\frac{1}{r\sin \unicode[STIX]{x1D703}}\unicode[STIX]{x2202}_{r}\unicode[STIX]{x1D6F9},\end{eqnarray}$$

in spherical coordinates. The streamfunction satisfies the biharmonic equation  $\unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D6F9}=0$ , the general solutions to which outside and inside the drop are (Kim & Karrila Reference Kim and Karrila2013):

(3.25a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F9}=\mathop{\sum }_{n=2}^{\infty }(A_{n}r^{-n+1}+B_{n}r^{-n+3}){\mathcal{G}}_{n}(\unicode[STIX]{x1D702}),\quad \bar{\unicode[STIX]{x1D6F9}}=\mathop{\sum }_{n=2}^{\infty }(\bar{A}_{n}r^{n}+\bar{B}_{n}r^{n+2}){\mathcal{G}}_{n}(\unicode[STIX]{x1D702}), & & \displaystyle\end{eqnarray}$$

where ${\mathcal{G}}_{n}(\unicode[STIX]{x1D702})$ are Gegenbauer functions of degree $-1/2$ of the first kind (Abramowitz & Stegun Reference Abramowitz and Stegun1972). They are related to Legendre polynomials and are regular everywhere in $-1\leqslant \unicode[STIX]{x1D702}\leqslant 1$ :

(3.26) $$\begin{eqnarray}\displaystyle {\mathcal{G}}_{n}(\unicode[STIX]{x1D702})=\frac{{\mathcal{L}}_{n-2}(\unicode[STIX]{x1D702})-{\mathcal{L}}_{n}(\unicode[STIX]{x1D702})}{2n-1},~~~n\geqslant 2. & & \displaystyle\end{eqnarray}$$

The first two functions are defined as ${\mathcal{G}}_{0}(\unicode[STIX]{x1D702})=1$ and ${\mathcal{G}}_{1}(\unicode[STIX]{x1D702})=-\unicode[STIX]{x1D702}$ , and we also note the property: ${\mathcal{G}}_{n}^{\prime }(\unicode[STIX]{x1D702})=-{\mathcal{L}}_{n-1}(\unicode[STIX]{x1D702})$ .

Following the same methodology as for the electric problem, we seek solutions as expansions in powers of $\unicode[STIX]{x1D6FF}$ . As will become evident in § 3.5 when performing the stress balance on the interface, the zeroth- and first-order electric stresses acting on the interface at most induce fluid motions of the form

(3.27a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F9} & = & \displaystyle \unicode[STIX]{x1D6F9}_{03}{\mathcal{G}}_{3}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D6FF}[\unicode[STIX]{x1D6F9}_{13}{\mathcal{G}}_{3}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D6F9}_{15}{\mathcal{G}}_{5}(\unicode[STIX]{x1D702})]+O(\unicode[STIX]{x1D6FF}^{2}),\end{eqnarray}$$

(3.27b ) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D6F9}} & = & \displaystyle \bar{\unicode[STIX]{x1D6F9}}_{03}{\mathcal{G}}_{3}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D6FF}[\bar{\unicode[STIX]{x1D6F9}}_{13}{\mathcal{G}}_{3}(\unicode[STIX]{x1D702})+\bar{\unicode[STIX]{x1D6F9}}_{15}{\mathcal{G}}_{5}(\unicode[STIX]{x1D702})]+O(\unicode[STIX]{x1D6FF}^{2}),\end{eqnarray}$$

(3.28a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F9}_{03}=A_{03}r^{-2}+B_{03},\quad \bar{\unicode[STIX]{x1D6F9}}_{03}=\bar{A}_{03}r^{3}+\bar{B}_{03}r^{5}, & & \displaystyle\end{eqnarray}$$

(3.28c,d ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F9}_{13}=A_{13}r^{-2}+B_{13},\quad \bar{\unicode[STIX]{x1D6F9}}_{13}=\bar{A}_{13}r^{3}+\bar{B}_{13}r^{5}, & & \displaystyle\end{eqnarray}$$

(3.28e,f ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F9}_{15}=A_{15}r^{-4}+B_{15}r^{-2},\quad \bar{\unicode[STIX]{x1D6F9}}_{15}=\bar{A}_{15}r^{5}+\bar{B}_{15}r^{7}. & & \displaystyle\end{eqnarray}$$

3.4 Kinematic boundary condition

The kinematic boundary condition relates the shape deformation to the fluid velocity so as to satisfy the no-slip and no-penetration boundary conditions at the interface. The streamfunction $\unicode[STIX]{x1D6F9}$ can be used to determine the normal and tangential components of the fluid velocity on the drop surface, which are obtained by combining (3.6a )–(3.6b ) and (3.24) as

(3.29a ) $$\begin{eqnarray}\displaystyle v^{n} & = & \displaystyle v_{02}^{n}{\mathcal{L}}_{2}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D6FF}[v_{10}^{n}+v_{12}^{n}{\mathcal{L}}_{2}(\unicode[STIX]{x1D702})+v_{14}^{n}{\mathcal{L}}_{4}(\unicode[STIX]{x1D702})]+O(\unicode[STIX]{x1D6FF}^{2}),\end{eqnarray}$$

(3.29b ) $$\begin{eqnarray}\displaystyle \bar{v}^{n} & = & \displaystyle \bar{v}_{02}^{n}{\mathcal{L}}_{2}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D6FF}[\bar{v}_{10}^{n}+\bar{v}_{12}^{n}{\mathcal{L}}_{2}(\unicode[STIX]{x1D702})+\bar{v}_{14}^{n}{\mathcal{L}}_{4}(\unicode[STIX]{x1D702})]+O(\unicode[STIX]{x1D6FF}^{2}),\end{eqnarray}$$

(3.29c ) $$\begin{eqnarray}\displaystyle v^{t} & = & \displaystyle v_{01}^{t}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})\sin \unicode[STIX]{x1D703}+\unicode[STIX]{x1D6FF}[v_{11}^{t}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})+v_{13}^{t}{\mathcal{L}}_{3}(\unicode[STIX]{x1D702})]\sin \unicode[STIX]{x1D703}+O(\unicode[STIX]{x1D6FF}^{2}),\end{eqnarray}$$

(3.29d ) $$\begin{eqnarray}\displaystyle \bar{v}^{t} & = & \displaystyle \bar{v}_{01}^{t}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})\sin \unicode[STIX]{x1D703}+\unicode[STIX]{x1D6FF}[\bar{v}_{11}^{t}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})+\bar{v}_{13}^{t}{\mathcal{L}}_{3}(\unicode[STIX]{x1D702})]\sin \unicode[STIX]{x1D703}+O(\unicode[STIX]{x1D6FF}^{2}).\end{eqnarray}$$

$v^{n}=\bar{v}^{n}=-\dot{\unicode[STIX]{x1D709}}$

(3.30a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle v_{02}^{n}=\bar{v}_{02}^{n}=\unicode[STIX]{x1D6FF}{\dot{f}}_{12},\quad v_{10}^{n}=\bar{v}_{10}^{n}=\unicode[STIX]{x1D6FF}{\dot{f}}_{20}, & \displaystyle\end{eqnarray}$$

(3.30c,d ) $$\begin{eqnarray}\displaystyle & \displaystyle v_{12}^{n}=\bar{v}_{12}^{n}=\unicode[STIX]{x1D6FF}{\dot{f}}_{22},\quad v_{14}^{n}=\bar{v}_{14}^{n}=\unicode[STIX]{x1D6FF}{\dot{f}}_{24}. & \displaystyle\end{eqnarray}$$

Similarly, the no-slip boundary condition $v^{t}=\bar{v}^{t}$ dictates that

(3.31a-c ) $$\begin{eqnarray}v_{01}^{t}=\bar{v}_{01}^{t},\quad v_{11}^{t}=\bar{v}_{11}^{t},\quad v_{13}^{t}=\bar{v}_{13}^{t}.\end{eqnarray}$$

The matching of orders in (3.30) might seem surprising at first due to the presence of terms involving $\unicode[STIX]{x1D6FF}$ on the right-hand side. However, it is the only possible solution, as the leading-order term in $\dot{\unicode[STIX]{x1D709}}$ involves $\unicode[STIX]{x1D6FF}$ . This implies that temporal derivatives of the shape functions in fact scale as $\unicode[STIX]{x1D6FF}^{-1}$ , suggesting that the characteristic time scale for the shape transient is not the Maxwell–Wagner relaxation time used here for non-dimensionalization. This point will be made clearer in § 4.2.

3.5 Dynamic boundary condition

We now proceed to enforce the dynamic boundary condition of (2.7), which in dimensionless form reads

(3.32) $$\begin{eqnarray}\displaystyle \boldsymbol{n}\boldsymbol{\cdot }(\unicode[STIX]{x27E6}\boldsymbol{T}^{E}\unicode[STIX]{x27E7}+Ma\,\unicode[STIX]{x27E6}\boldsymbol{T}^{H}\unicode[STIX]{x27E7})=\frac{1}{Ca_{E}}(\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{n}, & & \displaystyle\end{eqnarray}$$

and requires us to evaluate electric and hydrodynamic stresses on the interface.

3.5.1 Electric and hydrodynamic stresses

As previously shown in (2.9), the jump in electric tractions can be decomposed into tangential and normal components, both of which involve quadratic products of expansions derived above. The tangential component $q\boldsymbol{E}^{t}=qE^{t}\boldsymbol{t}$ is continuous and is expanded as

(3.33) $$\begin{eqnarray}qE^{t}=[qE^{t}]_{01}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})\sin \unicode[STIX]{x1D703}+\unicode[STIX]{x1D6FF}\{[qE^{t}]_{11}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})+[qE^{t}]_{13}{\mathcal{L}}_{3}(\unicode[STIX]{x1D702})\}\sin \unicode[STIX]{x1D703}+O(\unicode[STIX]{x1D6FF}^{2}).\end{eqnarray}$$

Similarly, the expansion for the jump in electric pressure in (2.9) is found to be

(3.34) $$\begin{eqnarray}\unicode[STIX]{x27E6}\,p^{E}\unicode[STIX]{x27E7}=\unicode[STIX]{x27E6}\,p^{E}\unicode[STIX]{x27E7}_{00}+\unicode[STIX]{x27E6}\,p^{E}\unicode[STIX]{x27E7}_{02}{\mathcal{L}}_{2}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D6FF}\{\unicode[STIX]{x27E6}\,p^{E}\unicode[STIX]{x27E7}_{10}+\unicode[STIX]{x27E6}\,p^{E}\unicode[STIX]{x27E7}_{12}{\mathcal{L}}_{2}(\unicode[STIX]{x1D702})+\unicode[STIX]{x27E6}\,p^{E}\unicode[STIX]{x27E7}_{14}{\mathcal{L}}_{4}(\unicode[STIX]{x1D702})\}+O(\unicode[STIX]{x1D6FF}^{2}).\end{eqnarray}$$

The expressions for the coefficients are provided in appendix A.

The jump in hydrodynamic tractions is evaluated using (3.6c )–(3.6d ), in which the requisite components of the stress tensor in spherical coordinates are obtained from the velocity components as

(3.35a,b ) $$\begin{eqnarray}\displaystyle T^{H,rr} & = & \displaystyle -p^{H}+2\unicode[STIX]{x2202}_{r}v^{r},\quad \bar{T}^{H,rr}=-\bar{p}^{H}+2\unicode[STIX]{x1D706}\unicode[STIX]{x2202}_{r}\bar{v}^{r},\end{eqnarray}$$

(3.35c,d ) $$\begin{eqnarray}\displaystyle T^{H,r\unicode[STIX]{x1D703}} & = & \displaystyle r^{-1}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}v^{r}+r\unicode[STIX]{x2202}_{r}(v^{\unicode[STIX]{x1D703}}r^{-1}),\quad \bar{T}^{H,r\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D706}[r^{-1}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}\bar{v}^{r}+r\unicode[STIX]{x2202}_{r}(\bar{v}^{\unicode[STIX]{x1D703}}r^{-1})],\end{eqnarray}$$

(3.35e,f ) $$\begin{eqnarray}\displaystyle T^{H,\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}} & = & \displaystyle -p^{H}+2r^{-1}(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}v_{\unicode[STIX]{x1D703}}+v_{r}),\quad \bar{T}^{H,\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}}=-\bar{p}^{H}+2r^{-1}\unicode[STIX]{x1D706}(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}\bar{v}_{\unicode[STIX]{x1D703}}+\bar{v}_{r}).\end{eqnarray}$$

The diagonal stress components $T^{H,rr}$ and $T^{H,\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}}$ involve the fluid pressure $p^{H}$ , which can be obtained from the velocity by integration of the momentum equation. After some algebra, the jumps in hydrodynamic stresses induced by the zeroth- and first-order streamfunctions $\unicode[STIX]{x1D6F9}_{0}$ , $\unicode[STIX]{x1D6F9}_{1}$ , scaled with $\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70F}_{MW}$ , are found as

(3.36a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}^{nn} & = & \displaystyle \unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{0}^{nn}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{1}^{nn}+O(\unicode[STIX]{x1D6FF}^{2})\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{00}^{nn}+\unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{02}^{nn}{\mathcal{L}}_{2}(\unicode[STIX]{x1D702})\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FF}\{\unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{10}^{nn}+\unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{12}^{nn}{\mathcal{L}}_{2}(\unicode[STIX]{x1D702})+\unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{14}^{nn}{\mathcal{L}}_{4}(\unicode[STIX]{x1D702})\}+O(\unicode[STIX]{x1D6FF}^{2}),\end{eqnarray}$$

(3.36b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}^{nt} & = & \displaystyle \unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{0}^{nt}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{1}^{nt}+O(\unicode[STIX]{x1D6FF}^{2})\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{01}^{nt}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})\sin \unicode[STIX]{x1D703}+\unicode[STIX]{x1D6FF}\{\unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{11}^{nt}{\mathcal{L}}_{1}(\unicode[STIX]{x1D702})+\unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{13}^{nt}{\mathcal{L}}_{3}(\unicode[STIX]{x1D702})\}\sin \unicode[STIX]{x1D703}+O(\unicode[STIX]{x1D6FF}^{2}),\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where the various coefficients can all be expressed in terms of $B_{03}$ , $B_{13}$ , $B_{15}$ , as detailed in appendix B.

3.5.2 Stress balance

The electric and hydrodynamic traction jumps can now be substituted into the stress balance (3.32) to satisfy the dynamic boundary condition. In the normal direction, the stress balance requires:

(3.37a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x27E6}\,p^{E}\unicode[STIX]{x27E7}_{00}+Ma\,\unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{00}^{nn} & = & \displaystyle \frac{2}{Ca_{E}},\end{eqnarray}$$

(3.37b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x27E6}\,p^{E}\unicode[STIX]{x27E7}_{02}+Ma\,\unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{02}^{nn} & = & \displaystyle \frac{4}{Ca_{E}}\unicode[STIX]{x1D6FF}f_{12},\end{eqnarray}$$

(3.37c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x27E6}\,p^{E}\unicode[STIX]{x27E7}_{10}+Ma\,\unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{10}^{nn} & = & \displaystyle -\frac{2}{Ca_{E}}\unicode[STIX]{x1D6FF}f_{12}^{2},\end{eqnarray}$$

(3.37d ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x27E6}\,p^{E}\unicode[STIX]{x27E7}_{12}+Ma\,\unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{12}^{nn} & = & \displaystyle \frac{4}{Ca_{E}}\unicode[STIX]{x1D6FF}\left(f_{22}-{\textstyle \frac{5}{7}}f_{12}^{2}\right),\end{eqnarray}$$

(3.37e ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x27E6}\,p^{E}\unicode[STIX]{x27E7}_{14}+Ma\,\unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{14}^{nn} & = & \displaystyle \frac{18}{Ca_{E}}\unicode[STIX]{x1D6FF}\left(f_{24}-{\textstyle \frac{2}{7}}f_{12}^{2}\right).\end{eqnarray}$$

In the tangential direction, it yields

(3.38a ) $$\begin{eqnarray}\displaystyle [qE^{t}]_{01}+Ma\,\unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{01}^{nt} & = & \displaystyle 0,\end{eqnarray}$$

(3.38b ) $$\begin{eqnarray}\displaystyle [qE^{t}]_{11}+Ma\,\unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{11}^{nt} & = & \displaystyle 0,\end{eqnarray}$$

(3.38c ) $$\begin{eqnarray}\displaystyle [qE^{t}]_{13}+Ma\,\unicode[STIX]{x27E6}T^{H}\unicode[STIX]{x27E7}_{13}^{nt} & = & \displaystyle 0.\end{eqnarray}$$

The above balances now allow us to define more explicitly the value of the small-deformation parameter $\unicode[STIX]{x1D6FF}$ . The driving force for the flow is the tangential electric stress $q\boldsymbol{E}^{t}$ , which according to (3.38) induces hydrodynamic tractions scaling with $O(Ma^{-1})$ . The magnitude of the resulting flow therefore is such that all normal tractions, both electric and hydrodynamic, in (3.37) are of order $O(1)$ . Balancing these tractions with surface tension forces thus requires us to choose $\unicode[STIX]{x1D6FF}\propto Ca_{E}$ . For consistency with previous small-deformation theories, we define $\unicode[STIX]{x1D6FF}$ explicitly as

(3.39) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}=\frac{3\,Ca_{E}}{4(1+2R)^{2}}. & & \displaystyle\end{eqnarray}$$

In particular, the assumption of small deformation yields no restriction on the magnitude of the Mason number, which is only constrained by the approximation discussed in § 3.2.2.

3.6 Nonlinear charge convection

As a final calculation, we determine the Legendre coefficients of the nonlinear convection term in the charge convection (2.5). The convection term is straightforward to calculate after applying the identity

(3.40) $$\begin{eqnarray}\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v})=qv^{n}(\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }\boldsymbol{n})+\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v}^{t}),\end{eqnarray}$$

in which the expansions for $q$ , $v^{n}$ , $\boldsymbol{v}^{t}=v^{t}\boldsymbol{t}$ and $\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }\boldsymbol{n}$ can be substituted together with

(3.41) $$\begin{eqnarray}\unicode[STIX]{x1D735}_{s}=[\unicode[STIX]{x1D644}-\hat{\boldsymbol{e}}_{r}\hat{\boldsymbol{e}}_{r}+\unicode[STIX]{x1D6FF}(\hat{\boldsymbol{e}}_{r}\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D703}}+\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D703}}\hat{\boldsymbol{e}}_{r})]\boldsymbol{\cdot }\unicode[STIX]{x1D735}+O(\unicode[STIX]{x1D6FF}^{2}).\end{eqnarray}$$

All calculations done, the relevant Legendre coefficients appearing in (3.21)–(3.23) for the dipole and octupole moments are found to be

(3.42a ) $$\begin{eqnarray}\displaystyle [\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v})]_{01} & = & \displaystyle -{\textstyle \frac{2}{5}}q_{01}B_{03}+{\textstyle \frac{6}{5}}q_{01}\unicode[STIX]{x1D6FF}{\dot{f}}_{12},\end{eqnarray}$$

(3.42b ) $$\begin{eqnarray}\displaystyle [\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v})]_{11} & = & \displaystyle {\textstyle \frac{2}{5}}q_{01}A_{13}+{\textstyle \frac{2}{5}}q_{11}A_{03}-{\textstyle \frac{6}{35}}q_{13}A_{03}-{\textstyle \frac{54}{35}}q_{01}A_{03}f_{12}\nonumber\\ \displaystyle & & \displaystyle +\,{\textstyle \frac{4}{5}}q_{01}\unicode[STIX]{x1D6FF}{\dot{f}}_{22}+{\textstyle \frac{4}{5}}q_{11}\unicode[STIX]{x1D6FF}{\dot{f}}_{12}+{\textstyle \frac{18}{35}}q_{13}\unicode[STIX]{x1D6FF}{\dot{f}}_{12}+{\textstyle \frac{38}{35}}q_{01}\unicode[STIX]{x1D6FF}f_{12}{\dot{f}}_{12},\end{eqnarray}$$

(3.42c ) $$\begin{eqnarray}\displaystyle [\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(q\boldsymbol{v})]_{13} & = & \displaystyle {\textstyle \frac{8}{5}}q_{01}A_{13}+{\textstyle \frac{4}{3}}q_{01}A_{15}+{\textstyle \frac{2}{3}}q_{01}B_{15}+{\textstyle \frac{8}{5}}q_{11}A_{03}\nonumber\\ \displaystyle & & \displaystyle +\,{\textstyle \frac{4}{15}}q_{13}A_{03}-{\textstyle \frac{104}{15}}q_{01}A_{03}f_{12}+{\textstyle \frac{6}{5}}q_{01}\unicode[STIX]{x1D6FF}{\dot{f}}_{22}+{\textstyle \frac{8}{9}}q_{01}\unicode[STIX]{x1D6FF}{\dot{f}}_{24}\nonumber\\ \displaystyle & & \displaystyle +\,{\textstyle \frac{6}{5}}q_{11}\unicode[STIX]{x1D6FF}{\dot{f}}_{12}+{\textstyle \frac{8}{15}}q_{13}\unicode[STIX]{x1D6FF}{\dot{f}}_{12}-{\textstyle \frac{4}{5}}q_{01}\unicode[STIX]{x1D6FF}f_{12}{\dot{f}}_{12}.\end{eqnarray}$$

4.1 Taylor deformation parameter

For easy comparison with previous theories and experiments, we introduce Taylor’s deformation parameter ${\mathcal{D}}$ , defined as

(4.1) $$\begin{eqnarray}{\mathcal{D}}=\frac{r^{\Vert }-r^{\bot }}{r^{\Vert }+r^{\bot }},\end{eqnarray}$$

where $r^{\Vert }$ and $r^{\bot }$ denote the length of the drop in directions parallel and perpendicular to the electric field, respectively. The sign of ${\mathcal{D}}$ distinguishes between oblate ( ${\mathcal{D}}<0$ ) and prolate ( ${\mathcal{D}}>0$ ) shapes. For an axisymmetric drop, $r^{\Vert }$ and $r^{\bot }$ are reached at $\unicode[STIX]{x1D703}=0$ and $\unicode[STIX]{x03C0}/2$ , respectively:

(4.2a ) $$\begin{eqnarray}\displaystyle r^{\Vert } & = & \displaystyle r(0)=1+\unicode[STIX]{x1D6FF}f_{12}+\unicode[STIX]{x1D6FF}^{2}\left(f_{20}+f_{22}+f_{24}\right)+O(\unicode[STIX]{x1D6FF}^{3}),\end{eqnarray}$$

(4.2b ) $$\begin{eqnarray}\displaystyle r^{\bot } & = & \displaystyle r(\unicode[STIX]{x03C0}/2)=1-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D6FF}f_{12}+\unicode[STIX]{x1D6FF}^{2}\left(f_{20}-{\textstyle \frac{1}{2}}f_{22}+{\textstyle \frac{3}{8}}f_{24}\right)+O(\unicode[STIX]{x1D6FF}^{3}),\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle {\mathcal{D}}={\textstyle \frac{3}{4}}\left[\unicode[STIX]{x1D6FF}f_{12}+\unicode[STIX]{x1D6FF}^{2}\left(f_{22}+{\textstyle \frac{5}{12}}f_{24}-{\textstyle \frac{1}{4}}f_{12}^{2}\right)\right]+O(\unicode[STIX]{x1D6FF}^{3}). & & \displaystyle\end{eqnarray}$$

4.2 First-order theory

We first summarize the first-order theory, which allows us to compare our results with those of Taylor (Reference Taylor1966), Esmaeeli & Sharifi (Reference Esmaeeli and Sharifi2011) and Lanauze et al. (Reference Lanauze, Walker and Khair2013). The zeroth-order stress balance (3.37b ) and (3.38a ), together with the dipole relaxation (3.21), provide three coupled equations for the three unknowns $B_{03}$ , $f_{12}$ and $P_{01}$ . We first eliminate $B_{03}$ by combining (3.37a ) and (3.38a ), and after manipulations we arrive at a coupled system of first-order ordinary differential equations of the form

(4.4) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\left[\begin{array}{@{}c@{}}P_{01}\\ f_{12}\\ \end{array}\right]={\mathcal{F}}_{1}(P_{01},f_{12};Ca_{E},Ma,R,Q,\unicode[STIX]{x1D706}),\end{eqnarray}$$

where ${\mathcal{F}}_{1}$ is a nonlinear function whose explicit form is cumbersome and is omitted here for brevity. These equations can be integrated numerically in time subject to initial conditions. In all of the results shown below, we assume that the drop surface is initially spherical and does not carry any charge at $t=0$ , which provides the initial conditions:

(4.5a,b ) $$\begin{eqnarray}P_{01}(0)=\frac{Q-1}{Q+2},\quad f_{12}(0)=0.\end{eqnarray}$$

Equations (4.4) can easily be compared to previous first-order theories. First, neglecting charge convection decouples the dipole evolution equation from the fluid problem, yielding the simple relaxation equation

(4.6) $$\begin{eqnarray}\displaystyle {\dot{P}}_{01}+P_{01}=\frac{1-R}{1+2R}, & & \displaystyle\end{eqnarray}$$

the solution to which is:

(4.7) $$\begin{eqnarray}\displaystyle P_{01}=\frac{1-R}{1+2R}+\left(\frac{Q-1}{Q+2}-\frac{1-R}{1+2R}\right)\text{e}^{-t}. & & \displaystyle\end{eqnarray}$$

Substituting (4.7) into (4.4) then yields a simplified model which is similar to that of Lanauze et al. (Reference Lanauze, Walker and Khair2013) when the effect of fluid inertia is negligible. If we further neglect charge relaxation, we can easily solve for the transient deformation parameter as

(4.8) $$\begin{eqnarray}{\mathcal{D}}(t)={\mathcal{D}}_{T}(1-\text{e}^{-t/\unicode[STIX]{x1D70F}_{d}})\quad \text{where }\unicode[STIX]{x1D70F}_{d}=\frac{a\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D6FE}}\frac{(19\unicode[STIX]{x1D706}+16)(2\unicode[STIX]{x1D706}+3)}{40(\unicode[STIX]{x1D706}+1)},\end{eqnarray}$$

which matches the result of Esmaeeli & Sharifi (Reference Esmaeeli and Sharifi2011). In particular, the viscous-capillary time scale $\unicode[STIX]{x1D70F}_{d}$ emerges as the characteristic time scale for shape deformations, which also rationalizes the seeming contradiction in the matching of terms in the kinematic boundary of (3.30). Here, ${\mathcal{D}}_{T}$ is the steady first-order deformation parameter first obtained by Taylor (Reference Taylor1966) as

(4.9) $$\begin{eqnarray}{\mathcal{D}}_{T}=\frac{9}{16}\frac{\unicode[STIX]{x1D6F7}_{T}}{(1+2R)^{2}}Ca_{E}\end{eqnarray}$$

in terms of Taylor’s discriminating function $\unicode[STIX]{x1D6F7}_{T}$ :

(4.10) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{T}=(1-R)^{2}+R(1-RQ)\left[2+\frac{3}{5}\frac{2+3\unicode[STIX]{x1D706}}{1+\unicode[STIX]{x1D706}}\right].\end{eqnarray}$$

Note that (4.8) predicts an exponential relaxation towards the steady drop shape and therefore fails to capture the non-monotonic transient deformation observed in experiments and simulations (Lanauze et al. Reference Lanauze, Walker and Khair2015) and also predicted by the full solution of (4.4), as we discuss in § 5.

4.3 Second-order theory

The first-order theory can then be improved by solution of the second-order equations, which involve the additional unknowns $B_{13}$ , $B_{15}$ , $f_{22}$ , $f_{24}$ , $P_{11}$ and $P_{13}$ . These are provided by the first-order normal and tangential stress balances of (3.37c ), (3.37e ), (3.38b ) and (3.38c ), together with the moment evolution (3.22)–(3.23). The flow unknowns $B_{13}$ and $B_{15}$ can be eliminated by manipulating (3.37c ) and (3.38b ) for $B_{13}$ , and (3.37e ) and (3.38c ) for $B_{15}$ . When combined with the moment evolution equations, this yields a system of differential equations of the form

(4.11) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\left[\begin{array}{@{}c@{}}P_{11}\\ P_{13}\\ f_{22}\\ f_{24}\end{array}\right]={\mathcal{F}}_{2}(P_{11},P_{13},f_{22},f_{24};P_{01},f_{12};Ca_{E},Ma,R,Q,\unicode[STIX]{x1D706}),\end{eqnarray}$$

where ${\mathcal{F}}_{2}$ is another nonlinear function. Once again, these equations can be integrated in time numerically to obtain the multipole moments as well as shape functions entering Taylor’s deformation parameter of (4.3). The initial conditions for these variables in the case of an initially spherical and uncharged drop are

(4.12) $$\begin{eqnarray}P_{11}(0)=P_{13}(0)=f_{22}(0)=f_{24}(0)=0.\end{eqnarray}$$

If charge convection is neglected, (3.22)–(3.23) for the moments become uncoupled from the flow problem and only involve electric parameters. At steady state, the first-order multipole moments are then obtained as

(4.13a,b ) $$\begin{eqnarray}\displaystyle P_{11}=\frac{6}{5}f_{12}\left(\frac{1-R}{1+2R}\right)^{2},\quad P_{13}=\frac{9}{5}f_{12}\frac{1-R}{1+2R}, & & \displaystyle\end{eqnarray}$$

which matches (25) and (26) in the work of Ajayi (Reference Ajayi1978). The numerical codes solving systems (4.4) and (4.11) are available upon request.

5.1 Effect of transient charge relaxation and shape deformation

In this section, we first neglect nonlinear charge convection and focus on the effects of transient charge relaxation and transient shape deformation alone. Here we adopt the experimental values of system 1b. The drop deformation is plotted as a function of time in figure 2 for three distinct cases. In figure 2( $a$ ), both nonlinear charge convection and transient charge relaxation are neglected. In this case, the only time dependence enters through the temporal derivatives of the shape functions. We find that the drop shape becomes oblate ( ${\mathcal{D}}<0$ ), and our theoretical results asymptote at long times towards the steady-state predictions of Taylor (Reference Taylor1966) and Ajayi (Reference Ajayi1978) at first and second order, respectively. Both steady states, however, overpredict the drop deformation, and it is found, rather curiously, that the theory performs more poorly at second order than at first order; this was already noted by Ajayi (Reference Ajayi1978) and is a consequence of neglecting charge convection, as further discussed below. The transient is also poorly captured: the model predicts a monotonic increase of the drop deformation towards the oblate steady state and fails to capture the initial dynamics seen in experiments, where the drop first adopts a prolate shape before becoming oblate. Figure 2( $b$ ) shows the opposite situation, in which transient shape deformation is neglected but transient charge relaxation is included. In this case, the shape instantaneously adjusts to the charge distribution, which explains the immediate deformation to a prolate shape at $t=0$ as a result of the instantaneous polarization of the drop according to (4.5). The deformation subsequently relaxes monotonically towards its steady oblate value. However, accounting for both transient phenomena in figure 2( $c$ ) captures the transient dynamics correctly while still evolving towards the steady deformation values of Taylor (Reference Taylor1966) and Ajayi (Reference Ajayi1978) in the absence of charge convection. These results underscore the importance of including all transient effects if one wants to capture the correct shape dynamics.

5.2 Effect of nonlinear charge convection

We now turn to the full theoretical model, which includes transient charge and shape relaxation as well as nonlinear charge convection. As we show here, the main effect of charge convection is to reduce the strength of the interfacial velocity, thereby causing oblate drops to deform less but prolate drops to deform more, in agreement with computational studies (Feng Reference Feng1999; Lanauze et al. Reference Lanauze, Walker and Khair2015). We first consider the dynamics in a relatively weak electric field using the parameters of system 1a in figure 3. First, we note in figure 3( $a$ ) that the boundary element simulations perform best and capture both the transient and the steady state with very good accuracy. Our small-deformation theory with charge convection also captures the transient very well, but still slightly overpredicts the steady deformation parameter, albeit not as much as the models of Taylor (Reference Taylor1966) and Ajayi (Reference Ajayi1978). Interestingly, we find that while the second-order theory of Ajayi is worse than the first-order theory of Taylor in the absence of charge convection, such is not the case in our model, where including second-order terms is seen to improve the solution. The poor performance of Ajayi’s model is a direct consequence of the neglect of charge convection, which results in a stronger dipole moment and in turn leads to larger deformations. Charge convection by the flow, on the other hand, causes the transport of positive and negative charges from the poles towards the equator, thus effectively reducing the induced dipole. This point is evident in figure 3( $b$ ), showing the steady charge distribution on the drop surface, where we see that the second-order theory with charge convection best approximates the charge profile from boundary element simulations. This numerical charge profile, however, exhibits a sharper transition from negative to positive values at the equator.

Figure 3.

( $a$ ) Deformation parameter ${\mathcal{D}}$ as a function of time for the parameters of system 1a. ( $b$ ) Steady interfacial charge profile. The plots show experimental results of Lanauze et al. (Reference Lanauze, Walker and Khair2015), fully nonlinear boundary element simulations, first- and second-order small-deformation theory (with nonlinear charge convection), and the steady results of Taylor (Reference Taylor1966) and Ajayi (Reference Ajayi1978) that neglected charge convection.

Figure 4.

( $a$ ) Deformation parameter ${\mathcal{D}}$ as a function of time for the parameters of system 1b. ( $b$ ) Steady interfacial charge profile. For these parameter values, the charge distribution predicted by the boundary element simulation develops a discontinuity at the equator. See supplementary online materials at https://doi.org/10.1017/jfm.2016.704 for a movie showing the dynamics and flow field in this case.

Figure 5.

( $a$ ) Deformation parameter ${\mathcal{D}}$ as a function of time for the parameters of system 1c. ( $b$ ) Steady interfacial charge profile. The steady deformation values predicted by the models of Taylor (Reference Taylor1966) and Ajayi (Reference Ajayi1978) in this case are $-0.75$ and $-1.40$ , respectively, and out of the frame of panel ( $a$ ). For these parameter values, the charge discontinuity at the equator is so severe that the boundary element simulations blow up before reaching steady state; in this case, the charge profile shown in ( $b$ ) corresponds to a time shortly before the instability develops.

The effect of increasing field strength is shown in figure 4, corresponding to system 1b. Unsurprisingly, stronger fields cause larger drop deformations, which are not as easily captured by the theory. While the boundary element simulation matches the experimental data quite well, our nonlinear small-deformation theory captures the transient well but shows a significant departure at steady state. Nevertheless, the second-order theory still outperforms all previous theoretical models. The difficulty in capturing the steady state accurately can be understood by considering the charge profile in figure 4( $b$ ), where a sharp gradient is observed across the equator in the numerical data from boundary element simulations. This sharp gradient cannot be captured using only two Legendre functions as in the expansion of (3.16), which explains the discrepancy. The problem becomes yet more severe in stronger fields, as shown in figure 5 in the case of system 1c. There, an actual discontinuity appears in the charge profile, leading to the very poor performance of small-deformation theories and to numerical instabilities in the boundary element simulation, which blows up before reaching steady state. The formation of a charge shock in strong fields was first observed in the simulations of Lanauze et al. (Reference Lanauze, Walker and Khair2015), who also were not able to resolve it numerically using their boundary element algorithm based on spline interpolation. The boundary element method used here and described in appendix C solves the charge conservation equation using finite volumes, and yet is still unable to capture the discontinuity, suggesting that higher-order non-oscillating numerical schemes should be employed towards this purpose (LeVeque Reference Leveque2002).

Figure 6.

( $a$ ) Deformation parameter ${\mathcal{D}}$ as a function of time for the parameters of system 4, which correspond to a steady prolate shape. ( $b$ ) Steady interfacial charge profile.

The case of prolate deformations is illustrated in figure 6 using the parameters of system 4. In this case, the drop deformation increases monotonically with time. The steady deformation parameter obtained by simulations with $Ma=0.5$ is ${\mathcal{D}}=0.27$ , which slightly exceeds the value of ${\mathcal{D}}=0.22$ found by Lac & Homsy (Reference Lac and Homsy2007), who neglected charge convection ( $Ma\rightarrow \infty$ ); the experiments of Ha & Yang (Reference Ha and Yang2000), for which the value of $Ma$ is unknown, reported a deformation of ${\mathcal{D}}=0.25$ . Our small-deformation theory only provides a modest improvement at steady state over the predictions of Taylor (Reference Taylor1966) and Ajayi (Reference Ajayi1978), again confirming that nonlinear charge convection has a weaker effect for prolate drops. This again can be rationalized by considering the interfacial charge profile in figure 6( $b$ ): convection by the flow is seen to cause charge accumulation at the drop poles, and thus does not cause any discontinuity as in the oblate case. Instead, the charge profile remains relatively smooth and therefore can be reasonably well approximated using Legendre polynomials.

Figure 7.

Steady drop deformation ${\mathcal{D}}$ as a function of electric capillary number $Ca_{E}$ for the parameters of: ( $a$ ) system 2a, ( $b$ ) system 2b, and ( $c$ ) system 3. The various models are compared to the experimental measurements of Tsukada et al. (Reference Tsukada, Katayama, Ito and Hozawa1993) and Salipante & Vlahovska (Reference Salipante and Vlahovska2010).

As a final test, we compare our theoretical and numerical predictions for the steady drop shapes with the experimental results of Salipante & Vlahovska (Reference Salipante and Vlahovska2010) for systems 2a and 2b, with oblate drops, and Tsukada et al. (Reference Tsukada, Katayama, Ito and Hozawa1993) for system 3, with prolate drops, in figure 7. The experimental system 2 used two different drop sizes but identical material properties. At a given value of the electric capillary number $Ca_{E}$ , increasing drop size is equivalent to decreasing the electric field or increasing the Mason number $Ma$ , which reduces the effect of charge convection. Charge convection is therefore more significant in figure 7( $a$ ) for the smaller drop size, and indeed departures of our numerical and theoretical results from the small-deformation theories of Taylor (Reference Taylor1966) and Ajayi (Reference Ajayi1978) are more significant in this case. In both cases, our theoretical model performs quite well at predicting the steady drop shape, but still slightly overpredicts the experimental values, especially as $Ca_{E}$ increases; nonetheless the agreement is noticeably better than in previous models. The effect of charge convection is extremely weak in the case of prolate drops in system 3, as shown in figure 7( $c$ ). As a consequence, our first- and second-order theories are indistinguishable from Taylor’s and Ajayi’s results, and the deformations predicted by our axisymmetric boundary element method only slightly exceed the finite-element simulation values of Tsukada et al. (Reference Tsukada, Katayama, Ito and Hozawa1993), who neglected charge convection.