  We present a highly accurate numerical method based on a boundary integral
formulation and the leaky dielectric model to study the dynamics of
surfactant-covered drops in the presence of an applied electric field. The
method can simulate interacting 3D drops (no axisymmetric simplification) in
close proximity, can consider different viscosities, is adaptive in time and
able to handle substantial drop deformation. For each drop global
representations of the variables based on spherical harmonics expansions are
used and the spectral accuracy is achieved by designing specific numerical
tools: a specialized quadrature method for the singular and nearly singular
integrals that appear in the formulation, a general preconditioner for the
implicit treatment of the surfactant diffusion and a reparametrization
procedure able to ensure a high-quality representation of the drops also under
deformation. Our numerical method is validated against theoretical, numerical
and experimental results available in the literature, as well as a new
second-order theory developed for a surfactant-laden drop placed in a
quadrupole electric field.
