# A spectral boundary integral method for simulating electrohydrodynamic flows
in viscous drops

Mohammadhossein Firouznia  Department of Mechanical and Aerospace Engineering,
University of California San Diego, 9500 Gilman Drive, La Jolla, CA 92093, USA
Spencer H. Bryngelson  School of Computational Science & Engineering, Georgia
Institute of Technology, 756 West Peachtree Street NW, Atlanta, GA 30332, USA
David Saintillan  Corresponding author<EMAIL_ADDRESS>Department of Mechanical
and Aerospace Engineering, University of California San Diego, 9500 Gilman
Drive, La Jolla, CA 92093, USA

###### Abstract

A weakly conducting liquid droplet immersed in another leaky dielectric liquid
can exhibit rich dynamical behaviors under the effect of an applied electric
field. Depending on material properties and field strength, the nonlinear
coupling of interfacial charge transport and fluid flow can trigger
electrohydrodynamic instabilities that lead to shape deformations and complex
dynamics. We present a spectral boundary integral method to simulate droplet
electrohydrodynamics in a uniform electric field. All physical variables, such
as drop shape and interfacial charge density, are represented using spherical
harmonic expansions. In addition to its exponential accuracy, the spectral
representation affords a nondissipative dealiasing method required for
numerical stability. A comprehensive charge transport model, valid under a
wide range of electric field strengths, accounts for charge relaxation, Ohmic
conduction, and surface charge convection by the flow. A shape
reparametrization technique enables the exploration of significant droplet
deformation regimes. For low-viscosity drops, the convection by the flow
drives steep interfacial charge gradients near the drop equator. This
introduces numerical ringing artifacts we treat via a weighted spherical
harmonic expansion, resulting in solution convergence. The method and
simulations are validated against experimental data and analytical predictions
in the axisymmetric Taylor and Quincke electrorotation regimes.

##  I Introduction

A wide range of engineering applications involve liquid drops immersed in
another fluid while subject to an applied electric field. Some examples
include ink-jet printing [1], electrospraying [2], and microfluidic devices
and pumps [3]. These systems exhibit rich dynamics due to the electric field
and fluid flow coupling. When an interface separating two immiscible fluids is
subject to an otherwise uniform electric field, the electric field undergoes a
jump across the interface due to the mismatch in material properties. This
discontinuity in the electric field induces electric stresses that can deform
the interface and drive the fluid into motion.

We are interested in leaky dielectric liquids such as oils, which serve as
poor conductors. Unlike electrolyte solutions where diffuse Debye layers
affect the system’s dynamics, leaky dielectrics are characterized by the
absence of diffuse Debye layers [4]. The free charges instead concentrate on
the interfaces between different phases in the system. Consequently, the
electric field acting on the interfacial charge creates electric stresses
along the normal and tangential directions, which cause deformations and fluid
motion. Surface tension has a stabilizing effect in general, trying to restore
the equilibrium shape. Melcher and Taylor developed a framework for studying
electrohydrodynamic phenomena in leaky dielectric systems, known as the leaky
dielectric model (LDM) [5]. Central to their work is a charge conservation
model that describes a balance between Ohmic fluxes from the bulk, interfacial
charge convection, and finite charge relaxation. It was previously shown that
LDM can be derived asymptotically from electrokinetic models in the limit of
strong electric fields and thin Debye layers [6, 7].

This work focuses on the dynamics of a leaky dielectric drop immersed in
another dielectric fluid under a uniform DC electric field. This canonical
problem has been a long-standing research problem in electrohydrodynamics. In
his pioneering work, Taylor [8] formulated a small-deformation theory for an
isolated drop based on LDM and could predict oblate and prolate steady shapes
depending on the material properties. While Taylor’s theory shows good
agreement with experimental data in the limit of vanishing electric capillary
number $\mathrm{Ca}_{\mathrm{E}}$ (ratio of electric to capillary forces), the
discrepancy is significant at larger values of $\mathrm{Ca}_{\mathrm{E}}$.
Therefore, other researchers attempted to extend Taylor’s work by accounting
for second-order effects in $\mathrm{Ca}_{\mathrm{E}}$ [9], considering
spheroidal drops [10, 11], including inertial effects [12] and interfacial
charge convection [13, 14, 15, 16].

A variety of computational models have been developed to study drop dynamics
under strong electric fields at finite deformations, a problem untractable
using analytical theories. In the limit of negligible inertia, boundary
integral equations can be used to formulate and solve the coupled
electrohydrodynamic problem. Sherwood was the first to develop a boundary
element method for an axisymmetric drop in an equiviscous system and applied
it to capture breakup modes in prolate drops [17]. His original work was
subsequently extended to study drop pair interaction [18] and to cover a wider
range of fluid and electric parameters [19]. These earlier attempts used a
simplified boundary condition for the electric problem, which neglected
transient charge relaxation and interfacial charge convection by the flow.
These two effects have recently been shown to play a significant role in drop
dynamics and deformations [16]. Lanauze et al. [20] and Das and Saintillan
[21] recently addressed this problem and developed axisymmetric and three-
dimensional boundary element methods based on the full Melcher-Taylor LDM. The
effect of charge convection was specifically addressed in [21], where it was
shown to be responsible for Quincke electrorotation. These methods, however,
were found to lack accuracy and stability in the regime of strong electric
fields. Other numerical approaches have been used to study drop
electrohydrodynamics, including immersed boundary [22], level set [23, 24],
and finite element methods [25, 26, 27]. More recently, finite element
simulations [28, 29] were also used to investigate electrohydrodynamic
instabilities such as tip and equatorial streaming in drops under strong
electric fields. These latter techniques all include finite fluid inertia and,
with few exceptions [24], do not treat the drop surface as a sharp interface.

Improved accuracy within the boundary integral framework can be achieved using
spectral methods, which rely on expansions of the shape and interfacial
variables based on spherical harmonics. Such methods were recently developed
to simulate electrohydrodynamics of lipid vesicles [30] and also extended to
the case of individual drops and drop pairs [31, 32, 33]. These studies,
however, all neglected charge relaxation and charge convection and were thus
restricted to weak electric fields. Accurately capturing charge convection is
especially challenging as it nonlinearly couples fluid flow and charge
transport on a deformed interface. It can result in spurious aliasing errors
with negative consequences for accuracy and stability. This work addresses
this challenge and presents a spectral boundary integral method for the
electrohydrodynamics of deformable liquid drops based on the complete
Melcher–Taylor LDM. Interfacial charge convection is rigorously accounted for,
and dealiasing and reparametrization techniques are implemented to improve
accuracy and stability and enable long-time simulations.

The paper is organized as follows. We define the problem and discuss the
governing equations and boundary conditions in Sec. II.1, along with their
non-dimensionalization in Sec. II.2. Sec. II.3 presents the integral form of
the governing equations and boundary conditions used in developing the
boundary integral method. We discuss different aspects of the numerical method
in Sec. III: the spectral representation of all variables in terms of
spherical harmonics is discussed in III.1, followed by details of the
dealiasing method in Sec. III.2. Next, in Sec. III.3, we summarize the
numerical integration methods used in this study and correction methods to
ensure charge neutrality and incompressibility in Sec. III.4. As explained in
Sec. III.5, we also use a reparametrization method to improve the numerical
stability in simulations where the drop undergoes significant deformations. We
test and validate our computational model by applying it to a wide range of
dynamical behaviors, such as the axisymmetric Taylor regime under weak
electric fields in Sec. IV.1, and Quincke electrorotation under stronger
electric fields in Sec. IV.2. We also investigate the dynamics of low-
viscosity drops in Sec. IV.3, where charge convection plays an important role.
Finally, we discuss our conclusions and possible extensions of our work in
Sec. V.

##  II Problem definition

###  II.1 Governing equations

We consider a neutrally buoyant drop of a fluid occupying volume $V^{-}$
immersed in an infinite body of another fluid $V^{+}$ while subject to a
uniform electric field $\bm{E}_{\infty}=E_{\infty}\,\hat{\bm{e}}_{z}$ as
depicted schematically in Fig. 1. The interface $D$ separates the two fluid
media, and the surface unit normal $\bm{n}(\bm{x})$ is pointed towards the
suspending fluid. Initially, the drop is uncharged and spherical with radius
$r_{0}$. The material properties, namely the dielectric permittivities,
electric conductivities, and dynamic viscosities, are denoted by
$(\epsilon^{\pm},\,\sigma^{\pm},\,\mu^{\pm})$ inside and outside the drop,
respectively. Under the Taylor–Melcher leaky dielectric model [8], any net
charge in the system appears on the interface $D$, and the bulk of the fluids
remain electroneutral. Therefore, the electric potential is harmonic in the
bulk:

Figure 1:  Problem definition: a leaky dielectric drop with
$(\sigma^{-},\,\epsilon^{-},\,\mu^{-})$ is suspended in another leaky
dielectric fluid with $(\sigma^{+},\,\epsilon^{+},\,\mu^{+})$ and subject to
an external electric field $\bm{E}_{\infty}$. The drop deforms and diverges
from its initially spherical shape.

$\nabla^{2}\varphi^{\pm}(\bm{x})=0,\qquad\bm{x}\in V^{\pm}.$ (1)

Far away from the interface, the electric field $\bm{E}=-\bm{\nabla}\varphi$
tends to the applied electric field:

$\bm{E}^{+}\rightarrow\bm{E}_{\infty}=E_{\infty}\,\hat{\bm{e}}_{z},\qquad\text{as}~{}~{}|\bm{x}|\rightarrow\pm\infty.$
(2)

While the tangential component of the electric field is continuous across the
interface, its normal component undergoes a jump due to the mismatch in
material properties:

$\bm{n}\times\llbracket\bm{E}\rrbracket=\bm{0},\qquad\bm{x}\in D.$ (3)

We define the operator
$\llbracket\mathcal{F}\rrbracket\coloneqq\mathcal{F}^{+}-\mathcal{F}^{-}$ as
the jump in any variable $\mathcal{F}$ across the interface $D$. A surface
charge density develops at the interface and follows Gauss’s law,

$q(\bm{x})=\bm{n}\bm{\cdot}\llbracket\epsilon\bm{E}\rrbracket,\qquad\bm{x}\in
D.$ (4)

The surface charge evolves due to Ohmic currents from the bulk and convective
currents on the interface. Consequently, it satisfies the conservation
equation:

$\partial_{t}q+\bm{n}\bm{\cdot}\llbracket\sigma\bm{E}\rrbracket+\bm{\nabla}_{s}\bm{\cdot}(q\bm{u})=0,\qquad\bm{x}\in
D,$ (5)

where $\bm{\nabla}_{s}=({\bm{I}}-\bm{n}\bm{n})\bm{\cdot}\bm{\nabla}$ is the
surface gradient operator and $\bm{u}$ is the fluid velocity.

Neglecting the effect of inertia and gravity, the velocity and pressure fields
satisfy the Stokes and continuity equations:

${\mu^{\pm}}\nabla^{2}\bm{u}^{\pm}-\bm{\nabla}p^{\pm}=\mathbf{0},\qquad\bm{\nabla}\bm{\cdot}\bm{u}^{\pm}=0,\qquad~{}\bm{x}\in
V^{\pm}.$ (6)

The velocity vector is continuous across the interface and vanishes far from
it as

$\displaystyle\llbracket\bm{u}(\bm{x})\rrbracket=\bm{0},$
$\displaystyle\qquad\bm{x}\in D,$ (7)
$\displaystyle\bm{u}^{+}(\bm{x})\rightarrow\bm{0},$
$\displaystyle\qquad\text{as}~{}~{}|\bm{x}|\rightarrow\infty.$ (8)

The balance of interfacial forces requires that the jump in hydrodynamic and
electric tractions across the interface balance capillary forces:

$\llbracket\bm{f}^{\mathrm{H}}\rrbracket+\llbracket\bm{f}^{\mathrm{E}}\rrbracket=\gamma(\bm{\nabla}_{s}\bm{\cdot}\bm{n})\bm{n},\qquad\bm{x}\in
D.$ (9)

We neglect Marangoni effects due to variations in surface tension,
$\bm{\nabla}_{s}\gamma=\bm{0}$. Hydrodynamic and electric tractions are
expressed in terms of the Newtonian and Maxwell stress tensors, respectively:

$\displaystyle\bm{f}^{\mathrm{H}}=\bm{n}\bm{\cdot}\bm{T}^{\mathrm{H}},$
$\displaystyle\qquad\bm{T}^{\mathrm{H}}=-p\bm{I}+\mu\big{(}\bm{\nabla}\bm{u}+{\bm{\nabla}\bm{u}}^{T}\big{)},$
(10) $\displaystyle\bm{f}^{\mathrm{E}}=\bm{n}\bm{\cdot}\bm{T}^{\mathrm{E}},$
$\displaystyle\qquad\bm{T}^{\mathrm{E}}=\epsilon\left(\bm{E}\bm{E}-\tfrac{1}{2}E^{2}\bm{I}\right).$
(11)

The jump in electric tractions can be decomposed into tangential and normal
components as

$\bm{f}^{\mathrm{E}}=\llbracket\epsilon
E^{n}\rrbracket\bm{E}^{t}+\frac{1}{2}\llbracket\epsilon({E^{n}}^{2}-{E^{t}}^{2})\rrbracket\bm{n}=q\bm{E}^{t}+\llbracket
p^{\mathrm{E}}\rrbracket\bm{n},$ (12)

where $p^{\mathrm{E}}=\epsilon({E^{n}}^{2}-{E^{t}}^{2})/2$ is the electric
pressure [19]. The first term on the right-hand side represents the tangential
electric stresses in leaky dielectrics, and it vanishes when both fluids are
either perfect dielectrics or perfect conductors.

###  II.2 Non-dimensionalization

For the system described above, a dimensional analysis yields five
dimensionless groups, three of which characterize the mismatch of material
properties in the drop and the suspending fluid:

$\mathrm{R}=\frac{\sigma^{+}}{\sigma^{-}},\qquad\mathrm{Q}=\frac{\epsilon^{-}}{\epsilon^{+}},\qquad\mathrm{\lambda}=\frac{\mu^{-}}{\mu^{+}}.$
(13)

The limits of $\mathrm{\lambda}\rightarrow 0$ and $\infty$ correspond to a
bubble and a rigid particle, respectively. The remaining dimensionless groups
describe the system’s dynamics and can be obtained by comparing the
characteristic time scales in the problem. First, note that the response of
each fluid phase to Ohmic conduction is characterized by the charge relaxation
time:

$\tau_{c}^{\pm}=\frac{\epsilon^{\pm}}{\sigma^{\pm}}.$ (14)

The product $\mathrm{RQ}=\tau^{-}/\tau^{+}$ is the ratio of charge relaxation
times in two fluids and plays an important role in the dynamics of the drop
[21]. The polarization time for a rigid sphere under an applied electric field
is the Maxwell–Wagner relaxation time

$\tau_{\mathrm{MW}}=\frac{\epsilon^{-}+2\epsilon^{+}}{\sigma^{-}+2\sigma^{+}}=\mathrm{\frac{R(Q+2)}{1+2R}}\,\tau_{c}^{+},$
(15)

which provides an approximate timescale for polarization of the drop. The
accumulation of free charges on the interface creates electric forces that
drive the fluid into motion on the electrohydrodynamic time scale

$\tau_{\mathrm{EHD}}=\frac{\mu^{+}}{\epsilon^{+}{E^{2}_{\infty}}}.$ (16)

Deformations away from the equilibrium spherical shape relax under the effect
of surface tension on the capillary time scale

$\tau_{\mathrm{\gamma}}=\frac{\mu^{+}\,r_{0}}{\gamma}.$ (17)

By taking the ratios of these time scales, the two remaining dimensionless
groups can be defined as

$\mathrm{Ca}_{\mathrm{E}}=\frac{\tau_{\gamma}}{\tau_{\mathrm{EHD}}}=\frac{\epsilon{E^{2}_{\infty}}\,r_{0}}{\gamma}.\quad\mathrm{Ma}=\frac{\tau_{\mathrm{EHD}}}{\tau_{\mathrm{MW}}}=\frac{\mu^{+}}{{\tau_{\mathrm{MW}}{\epsilon^{+}}E^{2}_{\infty}}}.$
(18)

The electric capillary number $\mathrm{Ca}_{\mathrm{E}}$ compares electric
forces versus capillary forces, while the Mason number $\mathrm{Ma}$
characterizes the importance of charge conduction against surface charge
convection. Alternatively, we can construct $\mathrm{Ca}_{\mathrm{MW}}$ to be
independent of the electric field

$\mathrm{Ca}_{\mathrm{MW}}=\frac{\tau_{\mathrm{\gamma}}}{\tau_{\mathrm{MW}}}=\frac{\mu^{+}(1+\lambda)\,r_{0}}{\gamma\tau_{\mathrm{MW}}}=(1+\lambda)\,\mathrm{Ca}_{\mathrm{E}}\,\mathrm{Ma}.$
(19)

For a given set of material properties, changing $\mathrm{Ca}_{\mathrm{MW}}$
corresponds to varying the drop radius $r_{0}$.

We scale the governing equations and boundary conditions using length scale
$r_{0}$, time scale $\tau_{\mathrm{MW}}$, pressure scale
$\epsilon^{+}E^{2}_{\infty}$, and the characteristic electric potential
$E_{\infty}r_{0}$. In the remainder of this manuscript, all governing
equations and boundary conditions are dimensionless, and results are presented
in terms of the corresponding dimensionless variables.

###  II.3 Boundary integral formulation

The electric problem is formulated in integral form based on the solution to
Laplace’s equation as [17, 18, 19]

$\varphi(\bm{x}_{0})=-\bm{x}_{0}\bm{\cdot}\bm{E}_{\infty}-\int_{D}\bm{n}\bm{\cdot}\llbracket\nabla\varphi(\bm{x})\rrbracket\,\mathcal{G}\left(\bm{x}_{0};\bm{x}\right)\,\mathrm{d}s(\bm{x}),\qquad\text{for}~{}\bm{x}_{0}\in
V^{\pm},D,$ (20)

where the evaluation point $\bm{x}_{0}$ can be anywhere in space, and $\bm{x}$
denotes the integration point on the interface. The free-space Green’s
function for Laplace’s equation, $\mathcal{G}\left(\bm{x}_{0};\bm{x}\right)$,
captures the electric potential due to a point charge in an unbounded domain
as

$\mathcal{G}\left(\bm{x}_{0};\bm{x}\right)=\dfrac{1}{4\pi
r},\quad\text{where}~{}~{}\bm{r}=\bm{x}_{0}-\bm{x},~{}r=|\bm{r}|.$ (21)

Taking the gradient of Eq. (20) with respect to $\bm{x}_{0}$ and using Gauss’s
law (4), we derive an integral equation for the jump in the normal electric
field as a function of the surface charge distribution:

$\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{D}\llbracket
E^{n}(\bm{x})\rrbracket[\bm{n}(\bm{x}_{0})\bm{\cdot}\bm{\nabla}_{0}\mathcal{G}]\,\mathrm{d}s(\bm{x})-\mathrm{\dfrac{1+Q}{2(1-Q)}}\llbracket
E^{n}(\bm{x}_{0})\rrbracket=E^{n}_{\infty}(\bm{x}_{0})-\dfrac{q(\bm{x}_{0})}{1-\mathrm{Q}},\qquad\text{for}~{}\bm{x}_{0}\in
D.$ (22)

For a given charge distribution $q(\bm{x})$, Eq. (22) can determine
$\llbracket E^{n}(\bm{x})\rrbracket$, from which $E^{{n}{+}}$ and $E^{{n}{-}}$
follow as

$E^{{n}{+}}(\bm{x})=\dfrac{q(\bm{x})-\mathrm{Q}\llbracket
E^{n}(\bm{x})\rrbracket}{1-\mathrm{Q}},\qquad
E^{{n}{-}}(\bm{x})=\dfrac{q(\bm{x})-\llbracket
E^{n}(\bm{x})\rrbracket}{1-\mathrm{Q}}.$ (23)

The tangential electric field
$\bm{E}^{t}(\bm{x})=-\bm{\nabla}_{s}\varphi(\bm{x})$ can be computed by
differentiation of the electric potential (20) along the tangential direction.
The interfacial jump in electric tractions
$\llbracket\bm{f}^{\mathrm{E}}\rrbracket$ follows from Eq. (12) based on the
tangential and normal electric fields calculated at every point on the
interface. This can be used to determine the jump in hydrodynamic tractions
$\llbracket\bm{f}^{\mathrm{H}}\rrbracket$ using the dynamic boundary condition
(9) as

$\llbracket\bm{f}^{\mathrm{H}}\rrbracket=-\llbracket\bm{f}^{\mathrm{E}}\rrbracket+\mathrm{Ca}_{\mathrm{E}}^{-1}{(\bm{\nabla}_{{s}}\bm{\cdot}\bm{n})\bm{n}},$
(24)

which enters the calculation of the velocity field, as explained next.

The flow problem is also recast into a boundary integral form as [34, 35]

$\displaystyle\bm{u}(\bm{x}_{0})=$
$\displaystyle-\dfrac{1}{4\pi\,\mathrm{{Ma}\,(1+\lambda)}}\int_{D}\llbracket\bm{f}^{\mathrm{H}}(\bm{x})\rrbracket\bm{\cdot}\bm{G}(\bm{x}_{0};\bm{x})\,\mathrm{d}s(\bm{x})$
(25)
$\displaystyle+\mathrm{\dfrac{1-\lambda}{4\pi(1+\lambda)}}\,\,\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{D}\bm{u}(\bm{x})\bm{\cdot}\bm{T}(\bm{x}_{0};\bm{x})\bm{\cdot}\bm{n}(\bm{x})\,\mathrm{d}s(\bm{x}),\quad\text{for}~{}\bm{x}_{0}\in
D.$

Here, $\bm{G}$ is the free-space Green’s function for the Stokeslet or flow
due to a unit point force in an unbounded domain, and $\bm{T}$ is the
corresponding stress tensor:

$\bm{G}(\bm{x}_{0};\bm{x})=\dfrac{\bm{I}}{r}+\dfrac{\bm{r}\bm{r}}{r^{3}},\quad\bm{T}(\bm{x}_{0};\bm{x})=6\dfrac{\bm{r}\bm{r}\bm{r}}{r^{5}}.$
(26)

Note that the integral equations (20), (22) and (25) exhibit singular
behaviors of different orders as $\bm{x}$ approaches $\bm{x}_{0}$. This is due
to the singularity of the Green’s function for Laplace’s (21) and Stokes
equations (26). Conventional quadrature schemes have poor accuracy in the
presence of singular integrands and may not converge by increasing the level
of discretization. Therefore, accurate numerical integration requires special
treatment of the singularities, which we will discuss in Sec. III.3.

##  III Numerical methods

We solve Eqs. (22) and (25) by building upon a spectral boundary integral
method introduced by Zhao  _et al._ [36]. This method was also previously
applied to study vesicle dynamics under shear and extensional flows [37, 38]
as well as flows of confined red blood cells [39, 40, 41]. All variables,
including the interfacial shape, velocity, and charge, are represented using
truncated series of spherical harmonic expansions as discussed in Sec. III.1.
Nonlinear operations, geometrical quantities (such as mean curvature), and
spatial derivatives are computed accurately using a nondissipative dealiasing
method discussed in Sec. III.2. The boundary integrals are computed using a
quadrature scheme for the surface collocation points, with a special treatment
for the singular integrands, as summarized in Sec. III.3. Any changes to the
volume and net charge of the drop due to numerical errors are corrected as
explained in Sec. III.4 to ensure numerical stability over long simulation
times. In addition, we adopt a reparameterization technique in Sec. III.5 to
minimize the high-frequency component of the interfacial shape, which improves
the stability of the numerical method for cases with significant
deformations.The code associated with the methods and simulations in this
article is publicly available at https://github.com/mfirouzn/EHD_Drop_3D.

At $t=0$, the drop is uncharged with a spherical shape. The numerical
algorithm used in this study follows that of Das and Saintillan [21],
Firouznia and Saintillan [42], and Firouznia  _et al._ [43]. We perform the
following steps at every time iteration:

  1. 1.

Given the current charge distribution $q(\bm{x})$ and shape of the interface,
compute $\llbracket E^{n}(\bm{x})\rrbracket$ by numerically inverting (22)
using GMRES [44]. From $\llbracket E^{n}(\bm{x})\rrbracket$, we obtain
$E^{{n}+}(\bm{x})$ and $E^{{n}-}(\bm{x})$ via (23).

  2. 2.

Determine the potential $\varphi$ along the interface by evaluating (20).

  3. 3.

Differentiate the surface potential numerically along the interface in order
to obtain the tangential electric field
$\bm{E}^{\mathrm{t}}=-\bm{\nabla}_{{s}}\varphi$.

  4. 4.

Knowing both components of the electric field, determine the jump in the
electric traction $\llbracket\bm{f}^{\mathrm{E}}\rrbracket$ and use it to
obtain $\llbracket\bm{f}^{\mathrm{H}}\rrbracket$ using (24).

  5. 5.

Solve for the interfacial velocity using the Stokes boundary integral equation
(25).

  6. 6.

Compute $\partial_{t}q$ via (5) and update the charge distribution using a
second-order Runge-Kutta scheme.

  7. 7.

Update the position of the interface by advecting the grid with the normal
component of the interfacial velocity:
$\partial_{t}{\bm{x}}=(\bm{u}\bm{\cdot}\bm{n})\bm{n}$.

  8. 8.

Apply corrections to the shape and charge distribution to ensure
incompressibility and charge neutrality, as discussed in Sec. III.4.

  9. 9.

Reparametrize the interfacial shape following the method discussed in Sec.
III.5 to minimize high-frequency components in the spherical harmonic
expansion.

###  III.1 Surface representation

The shape of the drop is assumed to be smooth and of spherical topology.
Therefore, the surface is parameterized by a truncated series of spherical
harmonic expansion from a rectangular domain
$\mathbb{S}^{2}=\\{(\theta,\,\phi)|\,\theta\in(0,\,\pi),\,\phi\in[0,\,2\pi)\\}$
to $\mathbb{R}^{3}$:

$\bm{x}(\theta,\phi)=\sum_{n=0}^{N-1}\sum_{m=0}^{n}\,{\bar{P}}_{n}^{m}(\cos\theta)\left(\bm{a}_{nm}\,\cos
m\phi\,+\bm{b}_{nm}\,\sin m\phi\right),$ (27)

where $\theta$ and $\phi$ are the latitude and longitude angles, and
$\bm{s}=\\{\bm{a}_{nm},\,\bm{b}_{nm}\\}$ are the coefficients of the expansion
in a compact form [45]. The representation above yields $N^{2}$ spherical
harmonic modes per each component of $\bm{x}$ (total of $3N^{2}$ modes). The
normalized associated Legendre polynomials of degree $n$
$(n=0,\,1,\,2,\,\dots)$ and order $m$ $(m\leq n)$ are defined as

${\bar{P}}_{n}^{m}(\eta)=\frac{1}{2^{n}\,n!}\sqrt{\frac{(2n+1)\,(n-m)!}{2\,(n+m)!}}\,(1-\eta^{2})^{m/2}\,\frac{\mathrm{d}^{n+m}}{\mathrm{d}x^{\,n+m}}(\eta^{2}-1)^{n},$
(28)

and satisfy the orthogonality condition

$\int_{-1}^{1}\bar{P}^{m}_{n}(\eta)\,\bar{P}^{m}_{n^{\prime}}(\eta)\,\mathrm{d}\eta=\delta_{nn^{\prime}}.$
(29)

Similarly, the surface charge distribution is represented as

$q(\theta,\phi)=\sum_{n=0}^{N-1}\sum_{m=0}^{n}\,{\bar{P}}_{n}^{m}(\cos\theta)\left(\tilde{a}_{nm}\,\cos
m\phi\,+\tilde{b}_{nm}\,\sin m\phi\right).$ (30)

The rectangular domain $\mathbb{S}^{2}$ is discretized based on the roots of
the Legendre polynomial $P_{N}(\cos\theta)$ along $\theta$, and uniformly
along $\phi$. Forward and backward transformations are performed using the
SPHEREPACK library [46, 47]. Partial derivatives of a given distribution can
be computed using recurrence relations for the derivatives of the associated
Legendre polynomials [46, 48]. Besides the spectral accuracy, the spherical
harmonic representation allows for nondissipative dealiasing, which improves
the numerical stability of the simulations [36].

A local coordinate system is constructed at every point
$\bm{x}(\theta,\,\phi)$ on the surface of the drop, using two tangent vectors
$\bm{a}_{1,2}$ and the unit normal $\bm{a}_{3}$:

$\centering\bm{a}_{1}=\partial_{\theta}\bm{x},\qquad\bm{a}_{2}=\partial_{\phi}\bm{x},\qquad\bm{a}_{3}=\bm{n}=\dfrac{\bm{a}_{1}\times\bm{a}_{2}}{|\bm{a}_{1}\times\bm{a}_{2}|}.\@add@centering$
(31)

Consequently, the first and second fundamental forms of the drop surface have
the following components:

$L_{ij}=\bm{a}_{i}\bm{\cdot}\bm{a}_{j},\quad\text{and}\quad
B_{ij}=\bm{a}_{i,j}\bm{\cdot}\bm{n},\quad(i,j=1,2),$ (32)

which will be used in the subsequent derivations. Given a scalar function
$f(\theta,\,\phi)$ on the surface $D$ defined by the parameterization
introduced in (27), the surface gradient $\bm{\nabla}_{s}f$ is

$\bm{\nabla}_{s}\,f=\left(\dfrac{L_{22}\,\bm{a}_{1}-L_{12}\,\bm{a}_{2}}{W^{2}}\right)\partial_{\theta}f+\,\left(\dfrac{L_{11}\,\bm{a}_{2}-L_{12}\,\bm{a}_{1}}{W^{2}}\right)\partial_{\phi}f,$
(33)

where $W=(\det{\bm{L}})^{1/2}$ is the area element. Similarly, the surface
divergence of a vector field $\bm{v}(\theta,\,\phi)$ can be expressed as

$\bm{\nabla}_{s}\bm{\cdot}\bm{v}=\left(\dfrac{L_{22}\,\bm{a}_{1}-L_{12}\,\bm{a}_{2}}{W^{2}}\right)\bm{\cdot}\partial_{\theta}\bm{v}+\,\left(\dfrac{L_{11}\,\bm{a}_{2}-L_{12}\,\bm{a}_{1}}{W^{2}}\right)\bm{\cdot}\partial_{\phi}\bm{v}.$
(34)

In this study, we consider a drop of an incompressible fluid. In the absence
of Marangoni effects, the capillary stress is a function of the mean curvature
$H$:

$H=\frac{1}{2}\operatorname{Tr}(\bm{L}^{-1}\bm{B})=\frac{1}{2}\left(\dfrac{L_{22}\,B_{11}-2L_{12}\,B_{12}+L_{11}\,B_{22}}{W^{2}}\right).$
(35)

###  III.2 Aliasing errors

Samples of different functions may become indistinguishable on a discrete grid
by so-called aliasing [49]. This means, for instance, that a high-frequency
spherical harmonic basis function may be aliased to lower frequencies and
cause numerical instability in the simulations. Nonlinear operations and
differentiation broaden the deformation and charge distribution spectra on a
drop. As a result, the energy is moved to frequencies not resolved by the grid
resolution and may alias to the resolved frequencies [49]. Therefore, it is
physically consistent with the discretization level to remove the
corresponding energy from the solution.

Nonlinear operations that are susceptible to aliasing in our simulations
include the calculation of mean curvature, electrical stress, convective
charge flux, and quadrature. We use a mesh with finer resolution $M>N$ in the
simulations to prevent aliasing errors and filter the solution following any
nonlinear manipulation [36]. Our numerical experiments show that $M/N=2$ is
sufficient for lower-order and polynomial nonlinearities. For higher-order and
non-polynomial nonlinear manipulations such as those incurred in calculating
mean curvature and charge convection, we use an adaptive algorithm that
follows Rahimian  _et al._ [48]. Given a function $\bm{f}^{M}$ sampled over an
$M$-grid ($2M^{2}$ points in physical space and $M^{2}$ spherical harmonic
modes), one can interpolate the distribution on a finer $P$-grid
$(P/M=u_{f}>1)$ by upsampling:

$\bm{s}^{P}=\\{\bm{a}_{nm},\,\bm{b}_{nm}\\}^{P}=\begin{dcases}~{}\\{\bm{a}_{nm},\,\bm{b}_{nm}\\}^{M},&\quad\text{for
}~{}~{}n\leq M~{}~{}\text{and}~{}~{}0\leq m\leq n,\\\ ~{}0,&\quad\text{for
}~{}~{}n>M.\end{dcases}$

The upsampling factor $u_{f}$ can be determined based on the desired
tolerance. Fig. 2 shows how the error in the mean curvature decays as a
function of the finer grid resolution $P$ for oblate and prolate spheroids
representing typical drop shapes subject to an electric field.

Figure 2:  Error in the mean curvature as a function of the number of modes
used in the upsampled $P-$space. Two ellipsoids with $l/b=5$ (prolate) and
$l/b=0.2$ (oblate) are investigated, representing the typical shapes of a drop
due to EHD flows with large deformations. $\bm{H}^{*}$ is the mean curvature
computed on the finest upsampled grid and is chosen as a surrogate for the
exact curvature values.

###  III.3 Numerical integration

This section summarizes the numerical integration scheme used in this study,
which follows that used by Zhao  _et al._ [36]. The boundary integrals in
(20), (22), and (25) can be written in the general form

$I(\bm{x}_{0})=\int_{D}K(\bm{x},\bm{x}_{0})\,f(\bm{x})\,\mathrm{d}s(\bm{x})=\int_{\mathbb{S}^{2}}K(\bm{x}(\theta,\,\phi),\bm{x}_{0})\,f(\bm{x}(\theta,\phi))\,J(\bm{x}(\theta,\phi))\,\sin\theta\,\mathrm{d}\theta\,\mathrm{d}\phi,$
(36)

where $f(\theta,\phi)$ is a smooth scalar function over surface $D$, and $K$
is the kernel containing the Green’s function for Stokes or Laplace’s
equations.

The approximate spacing of a mesh with $N\times 2N$ points is
$h=\sqrt{A/2N^{2}}$, where $A$ is the drop surface area. Due to the singular
behavior of $K$ as $\bm{x}_{0}\rightarrow\bm{x}$, the drop surface $D$ is
divided into two regions where the integral $I$ has different behaviors.
Following the method of floating partition of unity [50, 51, 36], a local
polar patch is considered centered around $\bm{x}_{0}$. For any point $\bm{x}$
on $D$, $\rho(\bm{x},\bm{x}_{0})$ is defined as the distance along the great
circle that connects $\bm{x}$ to $\bm{x}_{0}$ on $\mathbb{S}^{2}$. In the next
step, a mask function is defined based on this coordinate as

$\eta(\rho)=\begin{dcases}~{}\exp{\left(\dfrac{2e^{-1/t}}{t-1}\right)},&\quad\text{for
}~{}~{}t=\rho/\rho_{1}<1,\\\ ~{}0&\quad\text{for
}~{}~{}\rho\geq\rho_{1},\end{dcases}$

where $\rho_{1}$ is the cutoff radius. This representation allows for accurate
calculation of the integral (36) as $\eta$ is a smooth function.

Next, the surface integral (36) is split into two parts:

$I=I_{1}+I_{2}=\int_{D}K(\bm{x},\bm{x}_{0})\,\eta(\rho(\bm{x},\bm{x}_{0}))\,f(\bm{x})\,\mathrm{d}s(\bm{x})+\int_{D}K(\bm{x},\bm{x}_{0})\,[1-\eta(\rho(\bm{x},\bm{x}_{0}))]\,f(\bm{x})\,\mathrm{d}s(\bm{x}).$
(37)

The integrand of $I_{1}$ has support only inside the patch.

Using a transformation to the local polar coordinate system $(\rho,\varphi)$,
we write:

$I_{1}=\int_{0}^{2\pi}\int_{D}^{\rho_{1}}K(\bm{x},\bm{x}_{0})\,\eta(\rho)f(\rho,\varphi)\sin{\rho}\,\mathrm{d}\rho\,\mathrm{d}\varphi,$
(38)

where the integrand is finite and periodic in $\varphi$. The integral (38) is
computed using Gauss quadrature along $\rho$ from $0$ to $\rho_{1}$, with a
uniform mesh in $\varphi$ from $0$ to $2\pi$. We set the patch radius on the
reference sphere to $\rho_{1}=\pi/\sqrt{N}$, which means the patch radius is
$O(h^{1/2})$ in $\mathbb{S}^{3}$. Inside the patch, $\sqrt{N}$ points are
considered along $\rho$, and $2\sqrt{N}$ points along $\varphi$. We note that
the quadrature points inside the patch do not coincide with the surface mesh
points. Therefore, bi-cubic spline interpolation is used to evaluate the
coordinates and other functions at quadrature points of $I_{1}$. The error of
singular integration with the mentioned choice of patch size is $O(h^{3})$
[51]. Higher-order accuracy can be achieved using a larger patch size at the
expense of computational cost.

The second part of integral (37) is $I_{2}$, which has a smooth integrand.
Therefore, it is computed accurately as:

$I_{2}\approx\sum_{i=1}^{N}\sum_{j=1}^{2N}K(\bm{x}_{ij},\bm{x}_{0})\,\eta(\rho(\bm{x}_{ij},\bm{x}_{0}))\,f_{ij}\,J_{ij}\,\mathit{w}_{ij},$
(39)

where $\bm{x}_{ij}=\bm{x}(\theta_{i},\,\phi_{j})$ are the quadrature points,
$w_{ij}$ are the corresponding weights and $J_{ij}$ is the Jacobian of the
transformation from $\mathbb{S}^{2}$ to $\mathbb{R}^{3}$. The quadrature (39)
converges exponentially with the mesh size $h$, using quadrature points with
Gaussian and uniform distributions along the $\theta$ and $\phi$ directions,
respectively. Further details on the numerical integration scheme can be found
in Zhao  _et al._ [36].

###  III.4 Incompressibility and charge neutrality

We expect no change in the drop volume, as both fluid phases are
incompressible. However, small changes in the volume occur due to numerical
errors, which we correct by adjusting the shape of the drop along the normal
direction [36]. In the reported simulations, the magnitude of these
adjustments is smaller than $10^{-8}r_{0}$ at every time step, where $r_{0}$
is the initial drop radius. Similarly, the net charge (less than
$~{}10^{-17}\varepsilon E_{0}$) is subtracted from the surface charge
distribution at every time step to ensure charge neutrality throughout the
simulations. The mentioned corrections are especially important for long
simulation times.

###  III.5 Reparametrization

The fluid-fluid interface of the drop evolves and deforms. During this
process, there is no physical mechanism to inhibit in-plane distortions of the
grid (i.e., depletion, aggregation, and skewness) since there is no bending
rigidity, in-plane shear resistance, or surface inextensibility, High-
frequency components of the spherical harmonic expansions thus grow,
exacerbating aliasing errors and resulting in numerical instability. In
addition to correcting aliasing errors, as explained in Sec. III.2, one must
develop a reparameterization strategy that ensures stable and accurate
simulations over long simulation times. Here, we use an algorithm that
minimizes the high-frequency components in the spherical harmonic expansion of
the surface parametrization. This method was introduced by Veerapenani et al.
[52, 48] and improved by Sorgentone and Tornberg [53].

Consider an implicit representation of the interface as a smooth function
$F:\mathbb{R}^{3}\mapsto\mathbb{R}$ such that $F(\bm{x})=0$ for all $\bm{x}\in
D$ where $D$ is the drop surface. The unit normal vector can be expressed as
$\bm{n}=\bm{\nabla}F/|\bm{\nabla}F|$ at every point on the surface. We define
a quality metric $E:\mathcal{X}\mapsto\mathbb{R}$ where $\mathcal{X}$ is the
space of smooth functions defined over $D$. Now, we can view the
reparameterization problem as a minimization of $E(\bm{x})$ subject to the
constraint $F(\bm{x})=0$ as

$\min_{\bm{x}\in D}\\{E(\bm{x})\\}~{}\text{subject to}~{}F(\bm{x})=0.$ (40)

It can be shown that the solution to the problem above is

$(\bm{I}-\bm{n}(\bm{x})\bm{n}(\bm{x}))\bm{\cdot}\bm{\nabla}E(\bm{x})=\bm{0},\quad\text{and}~{}~{}F(\bm{x})=0.$
(41)

The choice of quality metric $E(\bm{x})$ is not unique and could vary
depending on the reparametrization strategy. [52] proposed an equivalent of
the following quality metric:

$E(\bm{x})=\sum_{n=0}^{N-1}\sum_{m=0}^{n}\alpha_{nm}\left(|\bm{a}_{nm}|^{2}_{2}+|\bm{b}_{nm}|^{2}_{2}\right),$
(42)

where $|\bm{y}|_{2}$ is the $L^{2}$ norm of vector $\bm{y}$, and $\alpha_{nm}$
is the weight for $(n,m)$-th spherical harmonic. $E$, as defined in (42), can
be viewed as the spectral energy of the parametrization. We aim to minimize
the high-frequency part of $\bm{x}$. Therefore, the weights $\alpha_{nm}$ must
be small for low frequencies and larger for high frequencies. Here, we use the
following perfect low-pass filter

$\alpha_{nm}=\begin{dcases}1,&\quad\text{for }~{}~{}n>N_{\text{cutoff}},\\\
~{}0,&\quad\text{for }~{}~{}n\leq N_{\text{cutoff}}.\end{dcases}$

that was used [52, 53] to simulate vesicles and surfactant-laden drops and
improved long-time simulations. $N_{\text{cutoff}}$ is the cutoff frequency
and is chosen adaptively based on the energy spectrum as

$N_{\text{cutoff}}=\min\\{k\in\mathbb{N},~{}~{}1\leq k\leq
N-1~{}~{}|~{}~{}E_{k}/E_{1}\leq P_{\text{cutoff}}\\},$ (43)

where
$E_{k}=\sum_{n=k}^{N-1}\sum_{m=0}^{n}\alpha_{nm}(|\bm{a}_{nm}|^{2}_{2}+|\bm{b}_{nm}|^{2}_{2})$
and $P_{\text{cutoff}}$ determines the fraction of modes we penalize
($P_{\text{cutoff}}=0.2$ in our simulations) [53].

We solve (41) by marching in pseudo-time $\tau$ along the tangential direction
at every point:

$\partial_{\tau}\bm{x}+(\bm{I}-\bm{n}(\bm{x})\bm{n}(\bm{x}))\bm{\cdot}\bm{\nabla}E(\bm{x})=\bm{0},\quad\mathrm{with}\,\,\,\,\bm{x}(\tau=0)=\bm{x}_{0}.$
(44)

However, following (44) the volume is not necessarily conserved. In addition,
the charge distribution is distorted and hence not spectrally accurate. We use
the method of Sorgentone and Tornberg [53] to address these challenges. We
project the linear pseudo-velocity $\partial_{\tau}\bm{x}$ along the
latitudinal and longitudinal directions to obtain the angular pseudo-velocity
$(\partial_{\tau}\theta,\,\partial_{\tau}\phi)$. This strategy comprises
updates of the angular coordinates $(\theta,\,\phi)$ for every point and uses
the original spherical harmonic expansions for the shape and charge
distribution to interpolate the updated distributions. Using the chain rule,

$\partial_{\tau}\bm{x}=\frac{\partial\bm{x}}{\partial\theta}\,\frac{\partial\theta}{\partial\tau}+\frac{\partial\bm{x}}{\partial\phi}\,\frac{\partial\phi}{\partial\tau}=\bm{a}_{1}\,\partial_{\tau}\theta+\bm{a}_{2}\,\partial_{\tau}\phi.$
(45)

Next, we project $\partial_{\tau}\bm{x}$ along the tangential directions by
taking inner products with $\bm{a}_{1}$ and $\bm{a}_{2}$:

$\displaystyle\partial_{\tau}\bm{x}\bm{\cdot}\bm{a}_{1}=~{}~{}|\bm{a}_{1}|^{2}\,~{}~{}\partial_{\tau}\theta+\bm{a}_{1}\bm{\cdot}\bm{a}_{2}\,~{}\partial_{\tau}\phi,$
(46)
$\displaystyle\partial_{\tau}\bm{x}\bm{\cdot}\bm{a}_{2}=\bm{a}_{1}\bm{\cdot}\bm{a}_{2}\,~{}\partial_{\tau}\theta+~{}~{}|\bm{a}_{2}|^{2}\,~{}~{}\partial_{\tau}\phi.$
(47)

The solution to (46) and (47) is

$\partial_{\tau}\bm{\theta}=\bm{L}^{-1}\bm{w},$ (48)

where
$\partial_{\tau}\bm{\theta}=\\{\partial_{\tau}\theta,\,\partial_{\tau}\phi\\}^{T}$
and
$\bm{w}=\\{\partial_{\tau}\bm{x}\bm{\cdot}\bm{a}_{1},\,\partial_{\tau}\bm{x}\bm{\cdot}\bm{a}_{2}\\}^{T}$,
and $\bm{L}$ is the surface metric defined in (32). We update the angular
coordinates of every grid point by marching explicitly in pseudo-time to
obtain new angles $(\theta^{*},\,\phi^{*})$. Using the new angular
coordinates, we interpolate the position and surface charge at every point
based on the spherical harmonic expansions of the grid before
reparametrization (at $\tau=0$):

$\displaystyle\bm{x}^{*}(\theta^{*},\phi^{*})=\sum_{n=0}^{N-1}\sum_{m=0}^{n}\,{\bar{P}}_{n}^{m}(\cos\theta^{*})\left(\bm{a}^{0}_{nm}\,\cos
m\phi^{*}\,+\bm{b}^{0}_{nm}\,\sin m\phi^{*}\right),$ (49) $\displaystyle
q^{*}(\theta^{*},\phi^{*})=\sum_{n=0}^{N-1}\sum_{m=0}^{n}\,{\bar{P}}_{n}^{m}(\cos\theta^{*})\left(\tilde{a}^{0}_{nm}\,\cos
m\phi^{*}\,+\tilde{b}^{0}_{nm}\,\sin m\phi^{*}\right).$ (50)

Finally, a forward (FSHT) and backward (BSHT) spherical harmonic
transformation of $\bm{x}^{*}$ and $q^{*}$ yields distributions on a standard
grid that is spaced uniformly along $\phi$ and Gaussian along $\theta$:

$\displaystyle\bm{x}^{*}\,\xrightarrow{\text{FSHT}}\,\\{\bm{a}_{nm},\,\bm{b}_{nm}\\}\xrightarrow{\text{BSHT}}\,\bm{x},$
(51) $\displaystyle
q^{*}\,\xrightarrow{\text{FSHT}}\,\\{\tilde{a}_{nm},\,\tilde{b}_{nm}\\}\xrightarrow{\text{BSHT}}\,q.$
(52)

##  IV Numerical results

In this section, we compare the results of our numerical method with existing
theoretical solutions and previous computational and experimental studies on
electrohydrodynamic flows in drops. Table 1 lists all systems considered in
this study, including their dimensional and non-dimensional parameters.

Table 1:  Physical systems studied using numerical simulations, with their dimensional and non-dimensional parameters. $\mu^{+}={0.69}\,{\mathrm{Pa\,s}}$, $\gamma={4.5}\,\mathrm{mN\,m^{-1}}$, and $\epsilon^{0}={8.8542e-12}\,\mathrm{F\,m^{-1}}$ denotes the permittivity of vacuum.  system | $\mathrm{R}$ | $\mathrm{Q}$ | $\mathrm{\lambda}$ | $\mathrm{Ca}_{\mathrm{E}}$ | $\mathrm{Ma}$ | $\sigma^{+}~{}(\mathrm{S\,m^{-1}})$ | $\epsilon^{+}/\epsilon^{0}$ | $r_{0}~{}(\mathrm{mm})$
---|---|---|---|---|---|---|---|---
S1 | $0.57$ | $36.59$ | $1.41$ | $0.05-0.75$ | $0.21-3.25~{}$ | $4.5\times 10^{-11}$ | $5.3$ | $0.7$
S2 | $0.57$ | $36.59$ | $1.41$ | $0.05-0.75$ | $0.65-9.75~{}$ | $4.5\times 10^{-11}$ | $5.3$ | $2.1$
S3 | $0.57$ | $36.59$ | $14.12$ | $0.01-1.44$ | $0.28-5.0~{}$ | $4.5\times 10^{-11}$ | $5.3$ | $0.25,~{}0.75,~{}1.25,~{}1.75$
S4 | $0.60$ | $10.56$ | $0.07$ | $0.15-1.0$ | $0.43-2.86$ | $3.8\times 10^{-11}$ | $4.7$ | $1.9$

###  IV.1 Axisymmetric regime

Under weak electric fields ($\mathrm{Ca}_{\mathrm{E}}\ll 1$), the drop adopts
a steady axisymmetric shape and charge distribution, as shown in Fig. 3.
Following the pioneering work of Taylor [8], various analytical models have
been proposed to improve the predictions of the steady shape of a drop subject
to weak electric fields and small deformations [34, 16]. We characterize the
deviation from the spherical shape using Taylor’s deformation parameter
$\mathcal{D}$ defined as

$\mathcal{D}=\dfrac{l-b}{l+b},$ (53)

where $l$ and $b$ are the lengths of the major axes of the drop parallel and
perpendicular to the applied electric field. Figure 4 shows the steady-state
deformation $\mathcal{D}$ from our simulations compared to the predictions of
different small-deformation theories (SDT), prior numerical simulations, as
well as experimental data. As shown in Fig. 4, our simulation results closely
match the experimental data and past predictions. In agreement with past
studies [20, 16, 21], interfacial charge convection by the fluid flow plays a
significant role in the dynamics and tends to weaken drop deformations.

Figure 3: Axisymmetric regime: evolution of a drop of system S2 subject to a
uniform electric field with $\mathrm{Ca_{E}}=0.2$ ($N=8,~{}M=3N$). Colors show
the interfacial charge density $q$, while arrows show the interfacial fluid
velocity. Also see the corresponding video in the Supplemental Material.

Figure 4:  Steady drop deformation as a function of electric capillary number
$\mathrm{Ca}_{\mathrm{E}}$ for systems S1 in $(a)$ and S2 in $(b)$. Results
from our spectral boundary integral method (SBIM) are compared with the
experimental data (EXP) of Salipante and Vlahovska [54], along with
axisymmetric and three-dimensional boundary integral simulations (BIM) by Das
and Saintillan [21, 16] and various small-deformation theories (SDT) [8, 34,
16]. We use $N=8$ and $M=3N$ in the simulations shown here.

###  IV.2 Quincke regime

Upon increasing the intensity of the electric field, the drop can undergo a
symmetry-breaking bifurcation to a dynamical regime characterized by a
rotational component to the EHD flow, as shown in Fig. 5 [55, 56, 54]. When
$\mathrm{RQ}>1$, the induced flow is from the poles to the equator, and the
induced electric dipole is anti-parallel to the applied electric field. This
condition is unfavorable for stability as perturbations to the system result
in a destabilizing electric torque that can drive spontaneous rotation beyond
a critical electric field strength $E_{\mathrm{c}}$. Quincke first observed a
similar phenomenon in solid spheres subject to uniform electric fields [57].
In solid spheres, the threshold for Quincke electrorotation is given by [58]

$E_{\mathrm{c}}=\sqrt{\frac{2\mu^{+}}{3\epsilon^{+}\,\tau_{\mathrm{MW}}\,(\bar{\epsilon}-\bar{\sigma})}},\quad\mathrm{where}\quad\bar{\sigma}=\frac{\sigma^{-}-\sigma^{+}}{\sigma^{-}+2\sigma^{+}},\quad\bar{\epsilon}=\frac{\epsilon^{-}-\epsilon^{+}}{\epsilon^{-}+2\epsilon^{+}}.$
(54)

Figure 5:  Quincke regime: spontaneous rotation and tilt of a drop of system
S3 subject to an applied electric field with $E_{\infty}/E_{\mathrm{c}}=1$ and
$\mathrm{Ca}_{\mathrm{MW}}=1.32$ ($N=10,~{}M=3N$). The tilt angle $\alpha$ is
measured from the initial orientation of the drop as depicted. See the
Supplemental Material for the corresponding video showing the evolution of the
drop in the Quincke regime.

Figure 5 shows the velocity field and charge distribution in a leaky
dielectric drop during Quincke electrorotation. Following the onset of
instability, the drop tilts away from its initial orientation. The tilt angle
$\alpha$ and deformation oscillate until they reach their steady-state values.
Figure 6 shows the steady-state tilt angle and deformation as a function of
the applied electric field $E_{\infty}/E_{\mathrm{c}}$ for a drop of the
system S3. For a given set of material properties, different
$\mathrm{Ca}_{\mathrm{MW}}$ correspond to different drop sizes, and the effect
of capillary forces is stronger for smaller drops (small
$\mathrm{Ca}_{\mathrm{MW}}$). Therefore, we expect the drop tilt angle
$\alpha$ to approach that of a solid sphere $\beta$ [58, 54]:

$\beta=\frac{\pi}{2}-\arctan\left(\frac{E^{2}_{\infty}}{E^{2}_{\mathrm{c}}}-1\right)^{-\frac{1}{2}},\qquad\text{for}~{}~{}{E_{\infty}}\geq{E_{\mathrm{c}}},$
(55)

for sufficiently large $\mathrm{\lambda}$ and small
$\mathrm{Ca}_{\mathrm{MW}}$. This is verified in Fig. 6 (a), where the green
curve for $\mathrm{Ca}_{\mathrm{MW}}=0.44$ closely follows the black line
corresponding to $\beta$.

Figure 6:  Steady tilt angle $\alpha$ (a), and deformation $\mathcal{D}$ (b)
as functions of the applied electric field $E_{\infty}/E_{\mathrm{c}}$ for
different values of $\mathrm{Ca}_{\mathrm{MW}}$ (different $r_{0}$) for system
S3. The solid black line in (a) shows the tilt angle $\beta$ of a rigid
sphere. Filled markers in (b) show cases where Quincke electrorotation occurs.
We used $N=10$ and $M=3N$ in all simulations.

###  IV.3 Low-viscosity drops

When a leaky dielectric drop with $\mathrm{RQ}<1$ is subject to a uniform
electric field, it compresses along the field direction. The induced
quadrupolar flow is from the poles to the equator at the interface. Under weak
electric fields, the characteristic interfacial velocity is:

$u_{\mathrm{T}}=\dfrac{9}{10}\mathrm{\dfrac{R(RQ-1)}{(2R+1)^{2}\,(1+\lambda)}}\,r_{0}\tau^{-1}_{\mathrm{EHD}}$
(56)

based on Taylor’s classic solution [8]. According to Eq. (56), we expect the
EHD flow to be stronger in low-viscosity drops ($\lambda<1$) as
$u_{\mathrm{T}}\propto(1+\lambda)^{-1}$. As we increase the intensity of the
applied electric field, low-viscosity drops may undergo different types of EHD
instabilities, such as dimpling, equatorial streaming, or Quincke rotation,
depending on their material properties [59]. In other experiments, colloidal
particles adsorbed on a drop interface were observed to accumulate at the
equator and form a belt, which broke into vortices of particles [60, 61]. In
this section, we study the dynamics of an oblate drop with a small viscosity
ratio $\lambda<0.1$, represented by S4, which is similar to the systems
studied in [61].

Our results show that charge convection significantly affects the system’s
dynamics in this regime. As a result of the nonlinear coupling between the
flow and charge dynamics, strong charge gradients build up around the equator,
as shown in Fig. 7. The nonlinear steepening disappears when we switch off the
convective term in the charge conservation equation (5). Therefore, we can
conclude that the charge convection is responsible for the nonlinear behavior
mentioned. This can be verified in Fig. 7(c) by comparing the charge profiles
for simulations with and without convection. The convergent EHD flow sweeps
charges of opposite signs towards the equator, which is the stagnation line of
the flow where the electric field is locally tangent to the interface. It has
been shown that in systems where an interface is subject to a tangential
electric field and a converging flow, strong charge gradients develop around
the stagnation line due to the effect of charge convection [43].

The emergence of sharp local features in the charge profile poses a
fundamental challenge for our numerical method, as spectral methods are most
efficient for smooth field variables. It is known that large gradients or
discontinuities in a function result in Gibbs phenomenon or ringing artifacts
in its spherical harmonic representation [62]. Following the emergence of
ringing artifacts, spurious oscillations cannot be contained locally as
spherical harmonics are global basis functions. This phenomenon occurs due to
the inability of the finite spherical harmonic expansions to properly
represent the infinite (or very large) derivatives in the sharp regions of a
discontinuous (or nearly discontinuous) function.

Figure 7:  Charge dynamics in a low-viscosity drop of S4 with
$(\mathrm{Ca}_{\mathrm{E}},\,\mathrm{Ma})=(0.3,\,1.43)$: steady-state velocity
and charge density fields in the presence of charge convection (a), and with
no convection (b), and the corresponding charge profiles along the direction
of the applied electric field (c). We used $N=8,~{}M=3N$ in both simulations.

Figure 8:  WSH convergence study on a drop of S4 with
$(\mathrm{Ca_{E}},\,\mathrm{Ma})=(1.0,\,0.43)$: (a) Error in the steady-state
charge density distribution as a function of mesh resolution $N$ for a fixed
value of $\delta=10^{-3}$ ($N^{*}=42$). Open markers correspond to cases with
ringing artifacts. (b) Error in the steady-state charge density distribution
as a function of relaxation factor $\delta$ ($N^{*}=80$,
$\delta/\delta_{0}=0.75$ and $\delta_{0}=10^{-4}$ for
$\smash{q^{*}_{\mathrm{WSH}}}$). (c) Interfacial charge profiles corresponding
to (b). See the simulation video for $N^{*}=80$, $\delta/\delta_{0}=0.75$ and
$\delta_{0}=10^{-4}$ in the Supplemental Material.

For relatively small electric capillary numbers $\mathrm{Ca_{E}}<0.2$, we can
suppress ringing artifacts by using a finer mesh corresponding to a larger
number of spherical harmonics modes, a more restricted time step size, and,
therefore, a higher computational cost. However, due to the stronger effect of
charge convection at larger values of $\mathrm{Ca_{E}}$, refining the
resolution is insufficient to remove the ringing artifacts when
$\mathrm{Ca_{E}}$ is increased.

To resolve this problem, we propose a weighted spherical harmonic expansion
(WSH) for the charge density distribution, where the contribution from high-
frequency components is exponentially relaxed [63]. For a given charge density
distribution $q(\theta,\phi)$, the spherical harmonic expansion (30) gives
coefficients $\tilde{s}=\\{\tilde{a}_{nm},\,\tilde{b}_{nm}\\}$ where
$n=0,\,1,\,2,\,\dots,\,N-1$ and $m\leq n$. Accordingly, we define
$\text{WSH}\\!:q\mapsto q_{\mathrm{WSH}}$ as:

$q_{\mathrm{WSH}}(\theta,\phi)=\sum_{n=0}^{N-1}\sum_{m=0}^{n}\,e^{-n(n+1)\delta}\left(\tilde{a}_{nm}\,\cos
m\phi\,+\tilde{b}_{nm}\,\sin m\phi\right)\,{\bar{P}}_{n}^{m}(\cos\theta),$
(57)

where $\delta$ is the relaxation factor. Based on (57), we note that WSH does
not affect the mean, and $q_{\mathrm{WSH}}$ converges to $q$ in the limit of
$\delta\rightarrow 0$. In Fig. 8(a), we demonstrate that increasing the mesh
resolution at a fixed value of $\delta$ eliminates the ringing artifacts. In
addition, Fig. 8(b) shows that the steady-state charge profile at a fixed
resolution converges as the relaxation factor $\delta$ decreases. We combine
both approaches and increase $N$ while decreasing $\delta$ to refine the
numerical solution while avoiding the emergence of ringing artifacts. The
numerical solution with the finest resolution, $N^{*}=\text{max}(N)$, is then
chosen as a surrogate to the exact solution in Fig. 8(a,b). To compare the
charge profiles computed on different mesh resolutions, we upsample them to a
fine grid with $N^{u}>N^{*}$. Convergence of the charge profile is illustrated
in Fig. 8(c), where the coupling of the flow and charge dynamics is found to
result in the formation of strong gradients in the charge density profile near
the equator, consistent with previous studies [16]. We note that our
simulations for system S4 under electric capillary numbers of up to
$\mathrm{Ca_{E}}\sim O(1)$ reach a steady axisymmetric profile and do not
exhibit any electrohydrodynamic instability, suggesting that the equatorial
vortices observed by Ouriemi and Vlahovska [61] could be a consequence of the
colloidal particles present at the interface in their experiments.

##  V Conclusions

We developed a spectral boundary integral method for simulating
electrohydrodynamic flows in leaky dielectric viscous drops. The drop surface,
charge density, and all other variables are represented using truncated series
of spherical harmonic expansions. In addition to the excellent accuracy of the
spectral representation, it enables us to develop a nondissipative dealiasing
method required for numerical stability. The charge transport is modeled using
a conservation equation for Ohmic conduction from the bulk and surface charge
convection by the flow and finite charge relaxation. We employ a
reparametrization technique which allows us to explore regimes where drops
undergo significant deformations. Our results closely match the existing
experimental data and analytical predictions in the axisymmetric Taylor and
Quincke electrorotation regimes.The code used in this article is publicly
available at https://github.com/mfirouzn/EHD_Drop_3D.

Moreover, our simulations confirmed that the dynamics of low-viscosity drops
are strongly influenced by interfacial charge convection. In this regime, the
interplay between the flow and charge dynamics results in steep gradients in
the interfacial charge density. The development of these sharp features in the
charge density profile results in the formation of ringing artifacts due to
the Gibbs phenomenon. Increasing the mesh resolution would eliminate the Gibbs
phenomenon under relatively small electric capillary numbers
$\mathrm{Ca_{E}}>0.2$. However, the effect of nonlinear steepening by the flow
is sufficiently strong at larger $\mathrm{Ca_{E}}$ that increasing the
resolution does not avoid ringing. To solve this problem, we introduced a
weighted spherical harmonic transformation that serves as a relaxation method
to exponentially damp high-wave-number coefficients. The convergence of our
numerical solution with decreasing level of relaxation shows that it provides
a close approximation to the true solution in this regime. Further
characterizing the system’s dynamics in this regime will require more
sophisticated numerical methods for solutions containing discontinuities.

## Acknowledgements

The authors gratefully acknowledge funding from the US National Science
Foundation under Grant No. CBET-1705377.

## References

  * Basaran  _et al._ [2013] O. A. Basaran, H. Gao, and P. P. Bhat, Nonstandard inkjets, Annu. Rev. Fluid Mech. 45, 85 (2013).
  * Fernández de La Mora [2007] J. Fernández de La Mora, The fluid dynamics of Taylor cones, Annu. Rev. Fluid Mech. 39, 217 (2007).
  * Laser and Santiago [2004] D. J. Laser and J. G. Santiago, A review of micropumps, J. Micromech. Microeng. 14, R35 (2004).
  * Saville [1997] D. Saville, Electrohydrodynamics: the Taylor-Melcher leaky dielectric model, Annu. Rev. Fluid Mech. 29, 27 (1997).
  * Melcher [1969] J. R. Melcher, Electrohydrodynamics: a review of the role of interfacial shear stresses, Annu. Rev. Fluid Mech. 1, 111 (1969).
  * Schnitzer and Yariv [2015] O. Schnitzer and E. Yariv, The Taylor–Melcher leaky dielectric model as a macroscale electrokinetic description, J. Fluid Mech. 773, 1 (2015).
  * Mori and Young [2018] Y. Mori and Y. N. Young, From electrodiffusion theory to the electrohydrodynamics of leaky dielectrics through the weak electrolyte limit, J. Fluid Mech. 855, 67 (2018).
  * Taylor [1966] G. I. Taylor, Studies in electrohydrodynamics. I. The circulation produced in a drop by an electric field, Proc. R. Soc. A: Math. Phys. Eng. Sci. 291, 159 (1966).
  * Ajayi [1978] O. O. Ajayi, A note on Taylor’s electrohydrodynamic theory, Proc. R. Soc. A: Math. Phys. Eng. Sci. 364, 499 (1978).
  * Zabarankin [2013] M. Zabarankin, A liquid spheroidal drop in a viscous incompressible fluid under a steady electric field, SIAM J. Appl. Math. 73, 677 (2013).
  * Zhang  _et al._ [2013] J. Zhang, J. D. Zahn, and H. Lin, Transient solution for droplet deformation under electric fields, Phys. Rev. E 87, 043008 (2013).
  * Lanauze  _et al._ [2013] J. A. Lanauze, L. M. Walker, and A. S. Khair, The influence of inertia and charge relaxation on electrohydrodynamic drop deformation, Phys. Fluids 25, 112101 (2013).
  * Shkadov and Shutov [2002] V. Y. Shkadov and A. A. Shutov, Drop and bubble deformation in an electric field, Fluid Dyn. 37, 713 (2002).
  * Feng [2002] J. Q. Feng, A 2D electrohydrodynamic model for electrorotation of fluid drops, J. Colloid Interface Sci. 246, 112 (2002).
  * He  _et al._ [2013] H. He, P. F. Salipante, and P. M. Vlahovska, Electrorotation of a viscous droplet in a uniform direct current electric field, Phys. Fluids 25, 032106 (2013).
  * Das and Saintillan [2017a] D. Das and D. Saintillan, A nonlinear small-deformation theory for transient droplet electrohydrodynamics, J. Fluid Mech. 810, 225 (2017a).
  * Sherwood [1988] J. D. Sherwood, Breakup of fluid droplets in electric and magnetic fields, J. Fluid Mech. 188, 133 (1988).
  * Baygents  _et al._ [1998] J. C. Baygents, N. J. Rivette, and H. A. Stone, Electrohydrodynamic deformation and interaction of drop pairs, J. Fluid Mech. 368, 359 (1998).
  * Lac and Homsy [2007] E. Lac and G. M. Homsy, Axisymmetric deformation and stability of a viscous drop in a steady electric field, J. Fluid Mech. 590, 239 (2007).
  * Lanauze  _et al._ [2015] J. A. Lanauze, L. M. Walker, and A. S. Khair, Nonlinear electrohydrodynamics of slightly deformed oblate drops, J. Fluid Mech. 774, 245 (2015).
  * Das and Saintillan [2017b] D. Das and D. Saintillan, Electrohydrodynamics of viscous drops in strong electric fields: numerical simulations, J. Fluid Mech. 829, 127 (2017b).
  * Hu  _et al._ [2015] W. F. Hu, M. C. Lai, and Y. N. Young, A hybrid immersed boundary and immersed interface method for electrohydrodynamic simulations, J. Comput. Phys. 282, 47 (2015).
  * Bjørklund [2009] E. Bjørklund, The level-set method applied to droplet dynamics in the presence of an electric field, Comput. Fluids 38, 358 (2009).
  * Theillard  _et al._ [2019] M. Theillard, F. Gibou, and D. Saintillan, Sharp numerical simulations of incompressible two-phase flows, J. Comput. Phys. 391, 91 (2019).
  * Feng and Scott [1996] J. Q. Feng and T. C. Scott, A computational analysis of electrohydrodynamics of a leaky dielectric drop in an electric field, J. Fluid Mech. 311, 289 (1996).
  * Feng [1999] J. Q. Feng, Electrohydrodynamic behaviour of a drop subjected to a steady uniform electric field at finite electric reynolds number, Proc. R. Soc. A: Math. Phys. Eng. Sci. 455, 2245 (1999).
  * Supeene  _et al._ [2008] G. Supeene, C. R. Koch, and S. Bhattacharjee, Deformation of a droplet in an electric field: Nonlinear transient response in perfect and leaky dielectric media, J. Colloid Interface Sci. 318, 463 (2008).
  * Collins  _et al._ [2013] R. T. Collins, K. Sambath, M. T. Harris, and O. A. Basaran, Universal scaling laws for the disintegration of electrified drops, PNAS 110, 4905 (2013).
  * Wagoner  _et al._ [2021] B. W. Wagoner, P. M. Vlahovska, M. T. Harris, and O. A. Basaran, Electrohydrodynamics of lenticular drops and equatorial streaming, J. Fluid Mech. 925 (2021).
  * Veerapaneni [2016] S. Veerapaneni, Integral equation methods for vesicle electrohydrodynamics in three dimensions, J. Comput. Phys. 326, 278 (2016).
  * Sorgentone  _et al._ [2019] C. Sorgentone, A. K. Tornberg, and P. M. Vlahovska, A 3D boundary integral method for the electrohydrodynamics of surfactant-covered drops, J. Comput. Phys. 389, 111 (2019).
  * Sorgentone  _et al._ [2021] C. Sorgentone, J. I. Kach, A. S. Khair, L. M. Walker, and P. M. Vlahovska, Numerical and asymptotic analysis of the three-dimensional electrohydrodynamic interactions of drop pairs, J. Fluid Mech. 914 (2021).
  * Sorgentone and Vlahovska [2022] C. Sorgentone and P. M. Vlahovska, Tandem droplet locomotion in a uniform electric field, arXiv e-prints , arXiv (2022).
  * Rallison and Acrivos [1978] J. M. Rallison and A. Acrivos, A numerical study of the deformation and burst of a viscous drop in an extensional flow, J. Fluid Mech. 89, 191 (1978).
  * Pozrikidis  _et al._ [1992] C. Pozrikidis _et al._ ,  _Boundary integral and singularity methods for linearized viscous flow_ (Cambridge university press, 1992).
  * Zhao  _et al._ [2010] H. Zhao, A. H. G. Isfahani, L. N. Olson, and J. B. Freund, A spectral boundary integral method for flowing blood cells, J. Comput. Phys. 229, 3726 (2010).
  * Bryngelson and Freund [2018a] S. H. Bryngelson and J. B. Freund, Floquet stability analysis of capsules in viscous shear flow, J. Fluid Mech. 852, 663 (2018a).
  * Bryngelson and Freund [2019] S. H. Bryngelson and J. B. Freund, Non-modal Floquet stability of a capsule in large amplitude oscillatory extension, Eur. J. Mech. B Fluids. 77, 171 (2019).
  * Bryngelson  _et al._ [2019] S. H. Bryngelson, F. Guéniat, and J. B. Freund, Irregular dynamics of cellular blood flow in a model microvessel, Phys. Rev. E 100, 012203 (2019).
  * Bryngelson and Freund [2018b] S. H. Bryngelson and J. B. Freund, Global stability of flowing red blood cell trains, Phys. Rev. Fluids 3 (2018b).
  * Freund [2013] J. B. Freund, The flow of red blood cells through a narrow spleen-like slit, Phys. Fluids 25, 110807 (2013).
  * Firouznia and Saintillan [2021] M. Firouznia and D. Saintillan, Electrohydrodynamic instabilities in freely suspended viscous films under normal electric fields, Phys. Rev. Fluids 6, 103703 (2021).
  * Firouznia  _et al._ [2022] M. Firouznia, M. J. Miksis, P. M. Vlahovska, and D. Saintillan, Instability of a planar fluid interface under a tangential electric field in a stagnation point flow, J. Fluid Mech. 931 (2022).
  * Saad and Schultz [1986] Y. Saad and M. H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Comput. 7, 856 (1986).
  * Boyd [2001] J. P. Boyd,  _Chebyshev and Fourier spectral methods_ (Courier Corporation, 2001).
  * Adams and Swarztrauber [1999] J. C. Adams and P. Swarztrauber, Spherepack 3.0: A model development facility, MWR 127, 1872 (1999).
  * Swarztrauber and Spotz [2000] P. N. Swarztrauber and W. F. Spotz, Generalized discrete spherical harmonic transforms, J. Comput. Phys. 159, 213 (2000).
  * Rahimian  _et al._ [2015] A. Rahimian, S. K. Veerapaneni, D. Zorin, and G. Biros, Boundary integral method for the flow of vesicles with viscosity contrast in three dimensions, J. Comput. Phys. 298, 766 (2015).
  * Canuto  _et al._ [2012] C. Canuto, M. Y. Hussaini, A. Quarteroni, A. Thomas Jr, _et al._ ,  _Spectral methods in fluid dynamics_ (Springer Science & Business Media, 2012).
  * Bruno and Kunyansky [2001] O. P. Bruno and L. A. Kunyansky, A fast, high-order algorithm for the solution of surface scattering problems: basic implementation, tests, and applications, J. Comput. Phys. 169, 80 (2001).
  * Ying  _et al._ [2006] L. Ying, G. Biros, and D. Zorin, A high-order 3D boundary integral equation solver for elliptic PDEs in smooth domains, J. Comput. Phys. 219, 247 (2006).
  * Veerapaneni  _et al._ [2011] S. K. Veerapaneni, A. Rahimian, G. Biros, and D. Zorin, A fast algorithm for simulating vesicle flows in three dimensions, J. Comput. Phys. 230, 5610 (2011).
  * Sorgentone and Tornberg [2018] C. Sorgentone and A. K. Tornberg, A highly accurate boundary integral equation method for surfactant-laden drops in 3D, J. Comput. Phys. 360, 167 (2018).
  * Salipante and Vlahovska [2010] P. F. Salipante and P. M. Vlahovska, Electrohydrodynamics of drops in strong uniform dc electric fields, Phys. Fluids 22, 112110 (2010).
  * Ha and Yang [2000] J. W. Ha and S. M. Yang, Electrohydrodynamics and electrorotation of a drop with fluid less conductive than that of the ambient fluid, Phys. Fluids 12, 764 (2000).
  * Sato  _et al._ [2006] H. Sato, N. Kaji, T. Mochizuki, and Y. H. Mori, Behavior of oblately deformed droplets in an immiscible dielectric liquid under a steady and uniform electric field, Phys. Fluids 18, 127101 (2006).
  * Quincke [1896] G. Quincke, Ueber rotationen im constanten electrischen felde, Annalen der Physik 295, 417 (1896).
  * Jones [1984] T. B. Jones, Quincke rotation of spheres, IEEE Trans. Ind. Appl. , 845 (1984).
  * Brosseau and Vlahovska [2017] Q. Brosseau and P. M. Vlahovska, Streaming from the equator of a drop in an external electric field, Phys. Rev. Lett. 119, 034501 (2017).
  * Dommersnes  _et al._ [2013] P. Dommersnes, Z. Rozynek, A. Mikkelsen, R. Castberg, K. Kjerstad, K. Hersvik, and J. Otto Fossum, Active structuring of colloidal armour on liquid drops, Nat. Commun. 4, 1 (2013).
  * Ouriemi and Vlahovska [2014] M. Ouriemi and P. M. Vlahovska, Electrohydrodynamics of particle-covered drops, J. Fluid Mech. 751, 106 (2014).
  * Gelb [1997] A. Gelb, The resolution of the Gibbs phenomenon for spherical harmonics, Math. Comput. 66, 699 (1997).
  * Chung  _et al._ [2005] M. K. Chung, S. M. Robbins, K. M. Dalton, R. J. Davidson, A. L. Alexander, and A. C. Evans, Cortical thickness analysis in autism with heat kernel smoothing, NeuroImage 25, 1256 (2005).