Immobilisation and remobilisation of surfactant-covered drop interfaces in electric fields: Simulations

Introduction

Surfactants are common in industrial and natural multiphase systems, where they are used to stabilise emulsions and foams, control coating and wetting processes, and regulate the motion and coalescence of drops and bubbles. [1, 2] Even when present at low concentrations, surfactants can strongly affect the flow because they preferentially adsorb at fluid interfaces. Spatial variations in surfactant concentration produce variations in surface tension and hence Marangoni stresses, which modify the motion of the interface and the surrounding fluids. [3]

A well-known consequence of this coupling is the immobilisation of a fluid interface by Marangoni stresses, as described in the classical theory of Frumkin and Levich 4. Tangential flow redistributes adsorbed surfactant and can cause it to accumulate near regions where the surface flow converges. The resulting surface-tension gradients oppose this motion, reducing the surface velocity and, when sufficiently strong, effectively immobilising the interface. The degree of immobilisation depends on both surfactant redistribution and the sensitivity of surface tension to concentration. Combining theory and experiments on rising surfactant-laden bubbles, Palaparthi et al. 5 examined the stagnant-cap regime, in which surfactant accumulation immobilises the rear of the bubble while the front end of the surface remains mobile. Their work identifies how surfactant transport and exchange control the extent of the stagnant region. For an insoluble surfactant on a drop in an extensional flow, Eggleton et al. 6 found stagnant surfactant-rich regions adjoining mobile, depleted regions at low coverage. At high coverage, the steep dependence of surface tension on concentration allowed strong Marangoni stresses and complete stagnation at steady state despite relatively small concentration variations. For flow over superhydrophobic surfaces, Mayer and Crowdy 7 likewise distinguished stagnant-cap formation from immobilisation associated with the approach to maximum surfactant packing and the resulting steep nonlinear equation of state.

Surfactant transport also influences drop deformation and migration. Stone and Leal 8 examined the competition between surfactant convection and dilution caused by the increase in interfacial area during deformation. For soluble surfactants, Chen and Stebe 9 showed how the adsorption framework and exchange kinetics affect the retardation of thermocapillary migration. Eggleton and Stebe 10 studied a deforming drop in an extensional flow and showed that adsorption and desorption can alter both the Marangoni stresses and the distribution of surface tension, producing a non-monotonic dependence of deformation on the mass-transfer rate at high coverage.

Interfacial mobility can be restored when transport relaxes the concentration gradients generated by the flow. Stebe et al. 11 demonstrated remobilisation for surfactants with rapid adsorption–desorption kinetics and sufficiently effective bulk transport. In their experiments, micellar solutions of Triton X-100 and Brij-35 recovered mobile interfaces at high concentrations, whereas the slowly desorbing protein bovine serum albumin did not. Stebe and Maldarelli 12 extended this approach to mixtures containing a retarding surface-active component. By adding a rapidly exchanging surfactant in sufficient quantity, competitive adsorption allowed it to dominate the interface and restore mobility. Surface diffusion provides another means of relaxing concentration gradients. The possibility that different combinations of surfactant properties and transport processes can produce similar observable flows underlies the “hidden-variable” perspective of Manikantan and Squires [2]. Measurements of interfacial motion alone may therefore be insufficient to identify the mechanisms controlling that motion.

Electric fields drive interfacial flow directly, providing a setting in which to study how surfactants modify circulation and interfacial mobility. In the Taylor–Melcher model of a leaky-dielectric drop, the bulk liquids are assumed to be electroneutral, so free charge is confined to the liquid–liquid interface and is supplied by weak but finite Ohmic currents in the two phases. 13–15 The electric field acting on this interfacial charge produces a tangential Maxwell stress that drives Taylor circulation inside and outside the drop. Depending on the electrical properties of the two liquids, the surface flow is directed either from the equator towards the poles or from the poles towards the equator. This circulation defines the non-rotating Taylor regime; normal electric stresses additionally deform a freely suspended drop. In the complete leaky-dielectric model, interfacial charge conservation balances the local rate of change of surface-charge density and surface-charge convection against the jump in normal Ohmic current. Charge convection can produce sharp charge-density gradients in the Taylor regime. 16–19 A surfactant-covered leaky-dielectric drop subject to a uniform electric field therefore contains two coupled interfacial fields: surface charge and surfactant concentration. The same interfacial velocity redistributes both quantities, while electric and Marangoni stresses feed back on that velocity. This creates a strongly nonlinear coupling between the electrical, hydrodynamic and surfactant dynamics.

The electrohydrodynamics of surfactant-covered drops has been studied mainly in regimes where surface-charge convection is neglected or treated as a weak perturbation. Ha and Yang 20 analysed the deformation and stability of a drop covered by insoluble surfactant in an applied electric field. Teigen and Munkejord 21 subsequently used numerical simulations to show that surfactant can either increase or decrease drop deformation depending on the direction of the Taylor circulation. Nguangia et al. 22 developed small-deformation and spheroidal models for the equilibrium electrodeformation of a surfactant-covered drop. In their non-diffusing equilibrium limit, the interface becomes completely immobilised: Marangoni stress balances the tangential electric stress and suppresses the Taylor circulation. Nganguia et al. 23 studied nonlinear surfactant advection and diffusion in electrodeformation, while Nganguia et al. 24 included adsorption, desorption and bulk surfactant transport; both studies neglected surface-charge convection. Three-dimensional boundary-integral simulations have also treated general drop shapes and interacting drops in surfactant electrohydrodynamics. 25,26 Charge convection has been retained together with surfactant transport in perturbative calculations; for example, Poddar et al. 27 considered a sedimenting drop using a small electric Reynolds-number expansion. These studies establish the influence of surfactant transport on electrically driven deformation and circulation, primarily under conditions of negligible or weak charge convection. Interfacial rheology also modifies electrically driven drop dynamics. Nganguia et al. 28 used a spheroidal model to show how surface shear and dilatational viscosities alter prolate drop deformation and circulation, with their effects depending on the electrical properties and viscosity ratio.

At stronger electric fields, surface-charge convection can destabilise the Taylor state and produce spontaneous Quincke rotation. 29 In the Quincke regime, the drop has a nonzero mean rotational motion, with a tilted electric dipole and an electric torque that balances viscous resistance. Straining circulation can coexist with this mean rotation. When charge relaxes more slowly inside the drop than in the suspending liquid, the induced dipole can oppose the applied field. A small tilt then creates an electric torque, and charge convection can reinforce the tilt until sustained rotation develops. Quincke rotation also provides a versatile platform for synthetic active matter. Near a solid boundary, rotating solid spheres translate as Quincke rollers, whose collective dynamics include directed motion30 and coherent vortices. 31 Subsequent developments involving Quincke rollers include run-and-tumble motion, 32 chiral pattern formation, 33 colloidal oscillations, 34 dynamic roller clusters, 35 magnetic control, 36 and anisotropic rolling. 37 Geometrical asymmetry can also convert Quincke rotation into propulsion in a bulk fluid. 38

Quincke rotation is also observed in liquid drops and experiments have reported steady rotation, tilted shapes, and, for sufficiently low drop viscosities, hysteretic onset. 39–42 Fully nonlinear three-dimensional simulations of the complete leaky-dielectric model reproduced the transition from Taylor circulation to Quincke rotation and demonstrated the central role of charge convection. 18 The three-dimensional small-deformation model of He et al. 43 retained charge convection associated with the rotational part of the interfacial velocity but neglected convection by the straining flow. The later theory of Das and Saintillan 44 retained the leading-order contribution of the straining flow to surface-charge convection. More recently, McDougall et al. 45 extended this framework to an insoluble surfactant-covered drop by coupling charge convection to surfactant transport and Marangoni stresses. Their near-spherical model retains the first few interfacial harmonics, applies when the straining-flow contribution to charge convection remains weak, and predicts that surfactant can shift the critical field for Quincke rotation.

To our knowledge, a corresponding fully resolved three-dimensional calculation coupling nonlinear charge convection and surfactant transport in the low-viscosity-ratio Quincke regime has not yet been reported. Such a calculation would require tracking a deforming interface while resolving sharp surface-charge gradients and strongly non-uniform surfactant distributions. Systematic parameter studies using pseudo-arclength continuation would consequently be computationally demanding. We therefore adopt a two-dimensional model with a fixed circular interface, motivated by the strong surface-tension regime in which deformation is small, while retaining nonlinear charge convection and surfactant transport. The surface charge and surfactant concentration are resolved around the interface without assuming weak charge convection or restricting the fields to a few harmonics. The circular-drop model of Feng 46 and the large electric Reynolds-number analysis of Yariv and Frankel 47 established that this geometry retains the nonlinear mechanism responsible for electrorotation. Lanauze et al. 16 and Das and Saintillan 17 reported the development of sharp surface-charge gradients under strong charge convection. More recently, Peng et al. 19 analysed charge-density blowup in the circular-drop problem without charge diffusion, identifying the equatorial singularity and the associated balance of charge transport and annihilation. Peng and Schnitzer 48 demonstrated that the onset of Quincke rotation at low viscosity ratio can become subcritical, producing a bistable field interval, qualitatively consistent with the experimental hysteresis reported by Salipante and Vlahovska 41. Despite these advances for clean drops, to our knowledge the fully nonlinear coupling between surface-charge convection and surfactant transport through the Taylor–Quincke transition has not yet been investigated in a spatially resolved calculation even within the fixed-circle model. The circular geometry therefore provides a tractable setting in which to address this problem, resolve the coupled interfacial distributions and follow steady solution branches through the onset of rotation.

We investigate how surfactant modifies these dynamics when charge convection is appreciable and Quincke rotation is possible. We examine the suppression of Taylor circulation by Marangoni stresses, the recovery of interfacial motion and straining circulation when rotation begins at weak surfactant diffusion, and the effects of surfactant on the critical field and the nature of the rotational bifurcation. Motivated by the hidden-variable perspective of Manikantan and Squires[2], we then compare surface diffusion with exchange through a well-mixed bulk reservoir to determine whether the two transport mechanisms can produce similar flow responses and whether their concentration distributions distinguish them. The model couples surface-charge conservation, including nonlinear charge convection, to surfactant transport and a Langmuir equation of state. Steady solutions and their linear stability are used to characterise immobilisation, remobilisation and the transition between the Taylor and Quincke regimes. These calculations extend the reduced descriptions discussed above to spatially resolved nonlinear interfacial transport within the fixed-circle model, providing a test of whether similar flow responses can conceal distinct surfactant-transport mechanisms.

Dimensional problem formulation

All quantities in this section are dimensional. Dimensionless variables are introduced in Sec. 3, where the same symbols are reused for brevity.

Geometry and electric problem

We consider an infinitely long cylindrical drop of radius $a$ suspended in an immiscible, unbounded liquid. The problem is two-dimensional in the plane normal to the cylinder axis. The drop occupies $V^{-}=\{r<a\}$ and the suspending liquid occupies $V^{+}=\{r>a\}$. Superscripts $-$ and $+$ denote fields inside and outside the drop, respectively. Subscripts $r$, $\theta$ and $n$ denote vector or tensor components and do not denote differentiation; derivatives are written explicitly using $\partial/\partial r$, $\partial/\partial\theta$ and $\partial/\partial t$. The two liquids are Newtonian leaky dielectrics with permittivities $\epsilon^{\pm}$, electrical conductivities $\sigma^{\pm}$ and dynamic viscosities $\mu^{\pm}$. A uniform direct-current electric field $E_{0}\mathbf{e}_{z}$ is imposed at infinity. Polar coordinates $(r,\theta)$ are centred on the drop, with $\theta=0$ along the positive $z$-direction, so that $z=r\cos\theta$ and $x=r\sin\theta$. The outward unit normal is $\mathbf{n}=\mathbf{e}_{r}$ and the unit tangent in the direction of increasing $\theta$ is $\mathbf{e}_{\theta}$. For any field with limiting values on the two sides of the interface we use the jump convention $\llbracket f\rrbracket=f^{+}-f^{-}$. We denote the angular average of an

Vector (PDF)

caption: Schematic of the two-dimensional problem. A surfactant-covered circular leaky-dielectric drop of radius $a$ is placed in a uniform electric field $E_0\mathbf{e}_z$. Representative surface charges and the internal Taylor circulation are shown for electrical parameters with $R_S>1$, for which the interfacial flow is directed from the poles towards the equator. The interface is fixed at $r=a$.

interfacial quantity by

The bulk liquids are assumed charge-free, so that the electric fields are irrotational and the electric potentials $\phi^\pm$ satisfy

The imposed field gives the far-field condition

The electric potential and tangential electric field are continuous at $r=a$,

and the free surface-charge density $q(\theta,t)$ is determined by Gauss’s law,

Here $E_n^\pm=\mathbf{E}^\pm\cdot\mathbf{n}$. The charge evolves through Ohmic currents from the two bulk phases, convection by the interfacial flow and a surface-diffusion term,

Here $\nabla_s=(\mathbf{I}-\mathbf{n}\mathbf{n})\cdot\nabla$ and $D_q^{\mathrm{dim}}$ is the dimensional surface-charge diffusivity. Surface-charge diffusion has been retained in earlier electrohydrodynamic calculations, including those of Collins et al. 49. In the circular-drop problem, weak charge diffusion also regularises the sharp charge-density structures associated with the singular diffusionless limit. 48 We retain this term to resolve the charge distribution at finite diffusivity under strong convection. On the fixed circle, Eq. (6) becomes

Creeping flow and interfacial stresses

Fluid inertia is neglected because the Reynolds numbers of the small drops in the viscous liquids of interest are typically very small. The velocity and pressure fields therefore obey the Stokes equations

with hydrodynamic stress

The exterior fluid is quiescent at infinity and the interior velocity is regular at the origin. Velocity is continuous across the interface, and the prescribed geometry imposes no penetration:

The electric stress is described by the Maxwell tensor

For a freely deforming interface, normal traction balance determines the evolving shape. Here the shape is prescribed and the normal velocity is set to zero instead. We represent the normal traction required to maintain this constraint by $\Lambda(\theta,t)$, so that the total traction balance is

where $\kappa=1/a$. Projecting the same traction balance onto the tangent $\mathbf{e}_\theta$, and using continuity of the tangential electric field together with Eq. (5), gives

Equation (13) displays the three tangential stresses central to the present problem: the hydrodynamic stress $\mathbf{e}_\theta\cdot\llbracket\mathbf{T}_{H}\rrbracket\cdot\mathbf{e}_r$, the electric shear stress $qE_\theta$, and the Marangoni stress $a^{-1}\partial\gamma/\partial\theta$.

Surfactant transport and equation of state

The interface is coated by surfactant of concentration $\Gamma(\theta,t)$. We consider both an insoluble layer and a soluble-surfactant exchange model in which the interface exchanges surfactant with a prescribed, spatially uniform bulk reservoir of concentration $c_b$. For finite-site adsorption with rate constant $k_a$ and desorption rate $k_d$, the net flux into the interface is

where

Thus the soluble model represents kinetic exchange with a well-mixed reservoir; no bulk surfactant concentration field is needed in this model. The insoluble limit is recovered by setting $k_{\mathrm{ex}}=0$ while retaining the chosen mean surface concentration.

Surface advection, surface diffusion and exchange with the well-mixed bulk reservoir give

where $D_\Gamma^{\mathrm{dim}}$ is the dimensional surface diffusivity of the surfactant. On the circle, this reduces to

Its angular mean satisfies

We initialise the mean concentration at $\Gamma_{\mathrm{eq}}$, so the mean remains fixed even though adsorption and desorption may occur locally. For a finite-coverage monolayer we use a Langmuir equation of state,

where $R_{\mathrm{g}}$ is the gas constant, $T$ the absolute temperature, $\Gamma_{\infty}$ the maximum surface concentration and $\gamma_0$ the clean-interface surface tension. We denote by $\gamma_{\mathrm{eq}}$ the surface tension evaluated at $\Gamma=\Gamma_{\mathrm{eq}}$.

Nondimensional formulation

Scales and dimensionless groups

Lengths are scaled with $a$ and time with the exterior charge-relaxation time

The electroviscous velocity scale is

This scale follows from balancing electric and viscous stresses. Although it depends on the applied field, we retain it for consistency with previous electrohydrodynamic studies. For comparisons between different applied fields, we also introduce the fixed electroviscous velocity

Here $E_0$ is the applied field and $E_{\mathrm{c,s}}$ is the critical field for Quincke rotation of the reference rigid cylinder, specified in Eq. (32). The tilded velocity is nondimensionalised by $U_{\mathrm{c,s}}$, whereas $\boldsymbol{u}$ is nondimensionalised by $U_{\mathrm{E}}$. The fixed scale is used in Figure 6. Figures S1–S7 of the Supplementary Material reproduce Figure 2, Figure 3, Figure 5, Figure 8, Figure 11, Figure 12 and Figure 13, respectively, with velocities, rotation rates and strain intensities expressed using the corresponding fixed reference scales. The electrohydrodynamic flow time associated with $U_{\mathrm{E}}$ is

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Electric potential, electric field and surface charge are scaled with $E_0a$, $E_0$ and $\varepsilon^{+}E_0$, respectively. Pressure and both hydrodynamic and electric stresses are scaled with $\varepsilon^{+}E_0^2$. Surface tension is scaled with $\gamma_{\mathrm{eq}}$ and surfactant concentration with $\Gamma_{\mathrm{eq}}$. The same symbols are used below for the corresponding dimensionless variables, and the interface is at $r=1$. The material ratios are

The inner charge-relaxation time is $\tau^{-}=\varepsilon^{-}/\sigma^{-}$. Hence

so $RS$ directly compares charge relaxation inside and outside the drop. The Quincke rotational instability requires $RS>1$, or equivalently $\tau^{-}>\tau^{+}$.

The electric Reynolds number compares the charge-relaxation time of the outer liquid with the electroviscous advection time,

The dimensionless surface diffusivities and exchange rate are

We also define the electric capillary number and surfactant elasticity number,

and their ratio, the Marangoni number,

Here Ma denotes the Marangoni number, and must not be confused with the Mason number used for the viscous-to-electric stress ratio in earlier electrohydrodynamic studies. 17,18,44,50 The fixed-circle approximation is motivated by $\mathrm{Ca}_{\mathrm{E}}\ll1$, provided that the local surface tension remains positive and does not become asymptotically small relative to $\gamma_{\mathrm{eq}}$. Circularity is imposed exactly in the present model, and the deformation correction is not computed. A useful physical interpretation follows by defining the Marangoni–viscous time

so that $\mathrm{Ma}=\tau_{\mathrm{EHD}}/\tau_{\mathrm{Ma}}$. At fixed material properties, drop size and reference surface concentration, $\operatorname{Re}_{\mathrm{E}}\propto E_0^2$ whereas $\mathrm{Ma}\propto E_0^{-2}$. It is therefore convenient to introduce the field-independent surfactant-strength parameter

Thus a sweep in applied field for a single drop–fluid system is made at fixed $B$, with $\mathrm{Ma}=B/\operatorname{Re}_{\mathrm{E}}$; keeping fixed field while changing $\operatorname{Re}_{\mathrm{E}}$ would instead change the material system. The rigid-cylinder critical field, including surface-charge diffusion, is

where $D_q^{\mathrm{dim}}$ denotes the dimensional surface-charge diffusivity. This threshold exists when $\varepsilon^{-}\sigma^{+}>\varepsilon^{+}\sigma^{-}$, equivalently $RS>1$. The dimensionless applied field is then

where $\operatorname{Re}_{\mathrm{E,c,s}}$ is the corresponding rigid-cylinder critical electric Reynolds number. For the baseline electrical system, $R=S=2$ and $D_q=10^{-3}$, $\operatorname{Re}_{\mathrm{E,c,s}}=1.5060$. The same rigid-cylinder critical field $E_{\mathrm{c,s}}$ is used as the reference throughout.

We denote by $E_{\mathrm{c,d}}$ the field at which the non-rotating Taylor state first becomes unstable to Quincke rotation. For a subcritical transition, stable rotating states also exist below $E_{\mathrm{c,d}}$, down to the saddle-node shown in Figure 3. The corresponding critical Reynolds numbers are $\operatorname{Re}_{\mathrm{E,c,s}}$ and $\operatorname{Re}_{\mathrm{E,c,d}}$. The threshold ratio $E_{\mathrm{c,d}}/E_{\mathrm{c,s}}$ depends on the electrical material ratios, the viscosity ratio, the surfactant coverage and elasticity, and the charge and surfactant transport parameters.

Finally, the equilibrium surface coverage is

For an insoluble monolayer, $\Gamma_{\mathrm{eq}}$ is the prescribed mean surface concentration. Increasing the mean concentration of a given surfactant at fixed $\Gamma_{\infty}$ changes both $B$ and $\chi$; the fixed-$\chi$ sweeps therefore do not represent variation of surfactant concentration alone. The model parameters $B$ and $\chi$ can be specified independently, and choosing $B=1$ does not determine the coverage. The quantity $B$ is an electric field-independent measure of surfactant strength. The dimensional Gibbs (or Marangoni) elasticity, evaluated at the mean concentration, is

Thus $B/(1-\chi)$, rather than $B$ alone, measures the surface elasticity at the mean concentration in units set by viscous stresses and the charge-relaxation time. [2] [51] The dimensionless Langmuir equation of state is

The equation of state requires $0\le\chi\Gamma<1$. For dilute surface coverage, with $\chi\ll1$ and $\Gamma=O(1)$, this reduces to the linear equation of state $\gamma\simeq1-\mathrm{El}(\Gamma-1)$. Table S1 in the Supplementary Material gives an illustrative dimensional realisation of selected parameters for stated fluid properties, temperature and drop radius; the accompanying text gives the conversion relations for the transport parameters.

dimensionless parameter	value	dimensional value
$\mathrm{Re}_{E}$	1	$E_{0}\simeq119\ \mathrm{V\,mm}^{-1}$
$B$	1	$\Gamma_{\mathrm{eq}}\simeq0.24\ \mathrm{\mu mol\,m}^{-2}$
$D_{\Gamma}$	$10^{-3}$	$D_{\Gamma}^{\mathrm{dim}}\simeq7.8\times10^{-10}\ \mathrm{m}^{2}\,\mathrm{s}^{-1}$
$K_{s}$	0.004	$k_{\mathrm{ex}}\simeq3.8\times10^{-3}\ \mathrm{s}^{-1}$
$\chi$	0.2	$\Gamma_{\infty}\simeq1.2\ \mathrm{\mu mol\,m}^{-2}$

Dimensionless equations and boundary conditions

The electric potentials satisfy

with

and, at $r=1$,

The dimensionless surface-charge balance is

For an initially neutral drop, integration of Eq. (41) together with Gauss’s law preserves zero net charge,

The Stokes equations become

The exterior flow tends to zero at infinity, and the interfacial conditions are

The dimensionless hydrodynamic and Maxwell stresses are

The complete interfacial traction balance is

Projecting the dimensionless traction balance onto the tangent and using the Langmuir equation of state gives

Because the interface radius is prescribed, the normal component of Eq. (50) determines only the constraint traction $\Lambda$. The equations for tangential surface motion contain $\mathrm{El}$ and $\mathrm{Ca}_{E}$ only through their ratio $\mathrm{Ma}$.

The surfactant transport equation is

with

for the equilibrium-mean states considered here. Setting $K_s=0$ recovers an insoluble surfactant with conserved mean concentration; setting $\mathrm{Ma}=0$ removes the Marangoni feedback and recovers the clean-interface electrohydrodynamic problem. Equations (37)–(53) define the nonlinear circular-drop model studied below. Because the interface is circular, the bulk electric and Stokes problems can be eliminated analytically, leaving a closed interfacial problem for $q(\theta,t)$ and $\Gamma(\theta,t)$. The numerical methods and validation checks are described in the Supplementary Material.

Results and discussion

We first examine the steady branches and the recovery of circulation at rotational onset, then quantify the flow and its temporal evolution, and finally examine the parameter dependence of the onset and compare surface diffusion with bulk exchange.

Parameter families and representative interfacial states

We first consider an insoluble surfactant ($K_s=0$) and then return to soluble exchange in § 4.6. Unless otherwise stated, the electrical ratios are $R=S=2$, the surface-charge diffusivity is $D_q=10^{-3}$ and the equilibrium coverage is $\chi=0.2$. The baseline surfactant-covered family has $M=1$ and $B=1$, while the baseline surface diffusivity is varied over $D_\Gamma=1$, $0.1$, $0.01$ and $0.001$. The rotation-branch comparison additionally includes $D_\Gamma=10$. Clean-interface calculations use the same electrical and viscous material properties with the Marangoni coupling removed ($\mathrm{Ma}=0$). Additional sweeps vary $M$, $B$ and $D_\Gamma$ independently in order to separate viscous, elastic and surfactant-transport effects. The reported thresholds and branches refer to the finite charge diffusivity $D_q=0.001$.

It is important to distinguish changing the electric field in one physical drop–fluid system from varying dimensionless parameters independently.

From (26), (29) and (31), varying only $E_0$ for a given drop–fluid pair increases $Re_E$ while decreasing $Ma$ so that $B=Re_E Ma$ remains constant. All field-dependent branches shown against $E_0/E_{c,s}$ therefore follow this fixed-material path. Surface diffusion $D_{\Gamma}$ and exchange $K_s$ are also defined using the outer charge-relaxation time and remain fixed when only the field is varied. The same rigid-cylinder value $E_{c,s}$ is used as the reference field throughout.

The examples below are supplemented by sweeps of viscosity ratio, surfactant strength and surface diffusivity in Figure 10–12. These calculations establish the parameter dependence within the chosen electrical system, $R=S=2$, at $\chi=0.2$ and $D_q=10^{-3}$. The comparison of surface diffusion with exchange is examined separately in Sec. 4.6.

Vector (PDF)

Immobilisation of Taylor circulation and remobilisation on rotation

The effect of surfactant is already clear in the interior flow fields of Figure 2. The top row shows the clean drop at four applied electric fields. Throughout the displayed range the clean $M=1$ drop remains in the non-rotating Taylor regime and retains the familiar Taylor circulation. At the same electric fields, an insoluble surfactant with $B=1$ and $D_{\Gamma}=0.1$ strongly suppresses the Taylor circulation before the onset of Quincke rotation. The electrically driven surface motion advects surfactant, and the resulting Marangoni stress opposes the electric shear stress. Once the surfactant-covered drop begins to rotate, a strong mean azimuthal flow appears and the straining circulation also increases. Quincke rotation therefore remobilises the nearly immobilised interface and restores a substantial straining circulation. We quantify the net rotation by the mean angular velocity

Because the drop is fixed and centred, translation is excluded; the zeroth Fourier mode of $u_{\theta}$ therefore represents rigid rotation.

The steady branches in Fig. 3 show that surfactant changes the critical field for Quincke rotation, the supercritical or subcritical character of the rotational bifurcation, and the shape of the rotating branch. Here, bistability denotes the coexistence, at the same applied field and parameter values, of two stable states with different rotation magnitudes: either a non-rotating Taylor state and a rotating state, or two rotating states with different amplitudes. For a clean drop, reducing $M$ raises the critical field above the rigid-cylinder value and, at sufficiently low viscosity ratio, produces a subcritical rotating branch whose unstable segment extends towards lower electric fields and joins a stable rotating segment at a saddle-node bifurcation, creating the bistable interval found in the clean circular-drop problem. [48] This coexistence is qualitatively consistent with the experimental hysteresis reported by Salipante and Vlahovska [41]. The right panel fixes $M=1$ and $B=1$ and varies $D_{\Gamma}$. At large $D_{\Gamma}$, the rotating branch remains close to the clean-drop branch. Decreasing $D_{\Gamma}$ progressively lowers the critical field for Quincke rotation. At $D_{\Gamma}=0.01$ and $0.001$ rotation begins below the rigid-cylinder reference field $E_{c,s}$. In addition to suppressing Taylor circulation, surfactant transport changes the rotational instability: for this family, reducing $D_{\Gamma}$ lowers the critical field and changes the onset from subcritical in the clean and sufficiently high-diffusion cases ($D_{\Gamma}=10$) to supercritical at lower diffusion. At $D_{\Gamma}=1$, the onset is supercritical, but the rotating branch develops an S-shaped segment at finite rotation amplitude, bounded by two saddle-node (fold) bifurcations. This produces a narrow field interval containing two stable rotating states of different amplitudes, separated by an unstable branch. A magnified view is shown in Fig. S8 of the Supplementary Material.

Figure 4 shows the interfacial fields underlying these changes for the baseline surfactant-covered drop. In the weak Taylor regime the charge and surfactant distributions are smooth and retain the reflection symmetries of the applied field. Closer to onset, charge convection produces a much sharper charge profile and the surfactant redistribution becomes stronger. In the Quincke regime the charge distribution is tilted relative to the applied field and the surfactant pattern is shifted around the interface, consistent with the broken reflection symmetry and finite mean rotation. These profiles provide a direct view of the two affected interfacial fields before introducing global measures of their effect on the flow.

Vector (PDF)

Movies S1 and S2 show the development of Taylor circulation and Quincke rotation for a clean drop. Movies S3 and S4 show the corresponding evolution with a weakly diffusing insoluble surfactant, including immobilisation and subsequent remobilisation in Movie S4. The movie parameters and initial conditions are given in the Supplementary Material.

Measures of interfacial and interior flow

We quantify interfacial motion using the RMS tangential speed and the RMS speed remaining after subtraction of mean rotation,

A signed average cannot measure the Taylor circulation because oppositely directed surface motions cancel. We decompose $u_{\theta}=\Omega+u_{\mathrm{str}}$ with $\langle u_{\mathrm{str}}\rangle=0$. The mean tangential surface velocity generates rigid-body rotation inside the drop, $u_{\mathrm{rot}}=\Omega r\mathbf{e}_{\theta}$. We call the remaining non-rigid flow the straining circulation; it can contain local vorticity as well as strain. The interfacial measure $U_{\mathrm{str}}$ quantifies its surface-speed amplitude and includes all resolved nonzero velocity harmonics. The decomposition gives $U_{\mathrm{rms}}^{2}=\Omega^{2}+U_{\mathrm{str}}^{2}$.

Figure 5 separates the total surface speed from its straining component. In the Taylor regime $\Omega=0$, so the two measures coincide. Low diffusion can suppress this interfacial flow by more than an order of magnitude relative to the clean interface. At the onset of Quincke rotation $U_{\mathrm{rms}}$ rises sharply because of the newly acquired mean rotation, but the change is not purely rigid-body motion: $U_{\mathrm{str}}$ also recovers. This is clearest for $D_{\Gamma}=0.01$ and $0.001$: the Taylor circulation is nearly arrested before onset, but a finite straining component reappears after rotation begins. The rotational instability therefore restores both mean azimuthal motion and a substantial circulation suppressed by Marangoni stresses.

These speeds are nondimensionalised by the field-dependent electro-viscous velocity $U_E=\epsilon^{+}E_0^2a/\mu^{+}$. A decrease in $U_{\mathrm{rms}}$ or $U_{\mathrm{str}}$ therefore means a decrease relative to this scale and need not imply slower dimensional speeds. For example, on the clean $M=1$ Taylor branch, $U_{\mathrm{rms}}$ decreases from $0.1164$ at $E_0/E_{c,s}=0.4$ to $0.0969$ at $1.0$, whereas the same speeds on the fixed scale $U_{c,s}$ increase from $0.0186$ to $0.0969$. Charge convection redistributes the electric traction and reduces the large-scale circulation relative to $U_E$. Its spatial distribution matters as well as its magnitude: finer angular variations produce weaker velocity responses. Indeed, the RMS electric traction in units of $\epsilon^{+}E_0^2$ increases slightly from $0.941$ to $0.967$ between these states, while the quadrupolar contribution to the surface speed decreases from $0.1160$ to $0.0938. For a surfactant-covered interface, Marangoni traction also contributes to the stress balance. Corresponding speed plots using the fixed scale $U_{c,s}$ are provided in the Supplementary Material.

Figure 6 compares the interfacial concentration and interior flow before and after the onset of rotation at $D_{\Gamma}=10^{-3}$. At $E_0/E_{c,s}=0.85$, just below $E_{c,d}/E_{c,s}=0.8663$, surfactant accumulates towards the equator and the Taylor circulation is nearly arrested. At $E_0/E_{c,s}=1.30$, the surfactant pattern shifts around the interface, while both mean rotation and straining circulation are present. To compare speeds across these fields, the figure uses the fixed scale $U_{c,s}$. The straining RMS speeds on this scale are $6.44\times10^{-4}$ and $8.50\times10^{-2}$ in the Taylor and Quincke states, respectively. The increase also remains substantial when each speed is nondimensionalised by its own $U_E$, giving $8.92\times10^{-4}$ and $5.03\times10^{-2}$. Mean rotation accounts for about $89\%$ of the squared total surface speed in the Quincke state. Subtracting it reveals the recovered straining circulation more clearly.

The recovered circulation coexists with a strongly nonuniform surfactant distribution. The concentration ranges are $0.411$–$1.555$ and $0.588$–$1.710$ in the Taylor and Quincke states, respectively. In the rotating state, electric and Marangoni tractions no longer nearly cancel, allowing substantial circulation despite the persistent concentration gradients. The spatial correlation coefficient of the mean-subtracted electric and Marangoni tractions changes from approximately $-1.000$ in the Taylor state to $-0.720$ in the Quincke state, indicating weaker anticorrelation of their spatial variations. The same Quincke concentration profile is shown in the total flow and its two velocity components in Fig. 6(b–d).

Figures 4 and 6 show that the rotating state has displaced charge and surfactant patterns together with an inclined straining circulation. The charge distribution, the surfactant distribution and the circulation need not have the same orientation. Their different evolution laws help explain why the electric and Marangoni stresses cease to cancel closely. Surface charge is supplied and relaxed by Ohmic currents, whereas insoluble surfactant is conserved. The electric traction is $\tau_E=qE_\theta$, with $E_\theta$ containing both the applied and induced fields, whereas the Marangoni traction is $\tau_M=-Ma\partial_\theta\Gamma/(1-\chi\Gamma)$. A displacement of the charge relative to the applied field therefore need not change the electric traction in the same way as a displacement of surfactant changes the Marangoni traction. The orientation of the resulting circulation is quantified below.

Figure 7 quantifies the orientation using the extensional eigenvector of the rate-of-strain tensor at the drop centre. Writing the second harmonic of the straining surface velocity as

Vector (PDF)

the extensional-axis angle is defined, for $A_2=(a_2^2+b_2^2)^{1/2}>0$, by

The angle is measured from $+z$ towards $+x$ and is defined modulo $180^\circ$. When $A_2=0$, the central rate of strain vanishes and this orientation is undefined. Only the second surface-velocity harmonic contributes to the rate-of-strain tensor at the drop centre. Higher harmonics remain in the complete circulation and contribute to strain away from the centre. The Taylor solutions shown here have $\alpha=0$. For the clean drop and the surfactant cases with $D_\Gamma=1$ and $0.1$, the central extensional axis has positive tilt along the positive-rotation branches over the displayed fields. At lower diffusivities, $D_\Gamma=0.01$ and $0.001$, the tilt initially becomes negative after the respective onsets at $E_{c,d}/E_{c,s}=0.8876$ and $0.8663$. The sampled minima are approximately $-22.9^\circ$ and $-37.4^\circ$, respectively. The tilt subsequently crosses zero near $E_0/E_{c,s}=1.14$ and becomes positive. The mean rotation remains counterclockwise throughout: the reversal concerns the orientation of strain, rather than the direction of rotation.

We also measure the area-averaged strain intensity in the interior disk. With the rate-of-strain tensor

we define

Rigid rotation has zero strain, so $\mathcal{S}$ measures the straining part of the interior flow. In the $(z,x)$ coordinates, the rotation tensor associated with the mean rotational component is

Its squared norm defines the mean-rotation intensity,

We then define

Thus $\Theta=1$ in the non-rotating Taylor regime and $\Theta=-1$ for pure rigid-body rotation; intermediate values compare the relative contributions of interior straining and mean rotation. The same recovery is visible in these measures of the interior flow. Figure 8 shows the area-averaged strain intensity $\mathcal{S}$ and $\Theta$. For weakly diffusing surfactant the Taylor regime has very small strain, consistent with the nearly immobilised interface. Upon the onset of Quincke rotation, $\mathcal{S}$ rises from this suppressed value before decaying at larger fields in the electric velocity scaling. At the same time, $\Theta$ decreases from its Taylor value $+1$ and becomes negative. The rotating state is therefore dominated by mean rotation but still contains a finite straining circulation, and the relative importance of the two changes along the branch.

Temporal evolution of the tangential stress balance

We follow the temporal evolution of the tangential stress balance in Figure 9 using the baseline electrical parameters $R=S=2$, with $M=1$, $B=1$, $\operatorname{Re}_E=1$, $\chi=0.2$, $D_q=10^{-3}$ and $K_s=0$. The field is applied at $t=0$ to an initially uncharged interface with uniform surfactant concentration.

Vector (PDF)

Since $RS>1$, the Taylor circulation is directed from the poles towards the equator.

The left panel gives the RMS Marangoni-to-electric traction ratio and the right panel the hydrodynamic-to-electric ratio. For a clean interface these ratios are zero and one, respectively. The ratios quantify how surfactant redistribution transfers the tangential stress balance from viscous to Marangoni stresses as the concentration gradients develop. These ratios of norms need not sum exactly to one because the local tractions have spatial structure and signs.

This behaviour is the finite-diffusion counterpart of the equilibrium immobilisation identified by Nganguia et al. 22. For a non-diffusing insoluble surfactant, they showed that surfactant redistribution can bring the interface to complete stagnation, so that the viscous tangential stress vanishes and the tangential electric stress is balanced solely by the Marangoni stress. The present calculations recover the approach to that stress balance as $D_{\Gamma}$ decreases. For any finite $D_{\Gamma}>0$, however, exact stagnation cannot persist in a steady state with a non-uniform surfactant distribution: in the absence of advection, surface diffusion would relax $\Gamma$ towards uniformity and eliminate the Marangoni stress. A small residual Taylor circulation therefore remains and continuously maintains the surfactant gradient.

Surfactant advection generates concentration gradients whose Marangoni stresses balance most of the tangential electric traction, leaving weak viscous stresses and a small circulation. Increasing surfactant diffusion or introducing sufficiently rapid exchange with the bulk weakens these gradients and restores interfacial motion. Quincke rotation provides a further route to remobilisation and recovery of straining circulation, as demonstrated in Figure 5 and Figure 6.

How surfactant transport controls the onset of Quincke rotation

The dependence of the critical field for Quincke rotation on surfactant strength is summarised in Figure 10. The trends depend strongly on both $M$ and $D_{\Gamma}$. At $M=1$, large $D_{\Gamma}$ smooths the surfactant distribution and the critical field remains well above the rigid-cylinder value. Increasing $B$ nevertheless lowers the critical field. At smaller $D_{\Gamma}$ the reduction is much stronger, and rotation can begin below $E_{c,s}$. The lowest-diffusion curves are non-monotonic in $B$: the critical field reaches a minimum at intermediate surfactant strength and then rises again. A similar, though progressively weaker, non-monotonicity persists for $M=10$ and 50. The effect of viscosity ratio is substantial; note that the $y$-axis range changes between the three panels of Figure 10. The non-monotonic dependence shows that surfactant cannot be represented simply as progressively rigidifying the interface. Changing $B$ modifies both the Taylor-regime base state and the coupled charge–surfactant dynamics at the onset of Quincke rotation.

For $M=1$, the sampled minimum is $E_{c,d}/E_{c,s}=0.7472$ at $B=0.4$ for $D_{\Gamma}=10^{-3}$, and 0.8476 at $B=0.5$ for $D_{\Gamma}=10^{-2}$. Throughout these sweeps $\chi=0.2$ is held fixed. Consequently, this non-monotonicity cannot be explained by simultaneous variation of $\Gamma_{\mathrm{eq}}$ in the numerator and denominator of the Langmuir stress law. Even at fixed $\chi$, changing $B$ changes the spatial base state and the linearised coupling that sets $E_{c,d}$. A small value of $B$ is also not necessarily a weak-feedback limit at small $D_{\Gamma}$: the leading immobilisation parameter scales as $B/D_{\Gamma}$, as derived below.

The two-parameter organisation is made explicit in Figure 11, which maps the rotation magnitude in the $(B,D_{\Gamma})$ plane at six fixed electric fields for $M=1$. The cyan points mark the onset of Quincke rotation. At $E_0/E_{c,s}=0.8$ the rotating region is confined to a very small pocket at low diffusion and intermediate $B$; the inset magnifies this pocket and shows that it is centred around $B\simeq0.4$ and $D_{\Gamma}$ of a few $10^{-3}$, consistent with the minimum in the low-diffusion onset curves of Figure 10. By $E_0/E_{c,s}=1$ the pocket has expanded markedly, and for fields of 1.2–1.8 the rotating state occupies an increasingly broad portion of the parameter plane. The onset boundary is strongly curved: neither increasing surfactant strength nor decreasing diffusion changes the onset monotonically over the whole parameter range. The boundary reflects the competition between surfactant advection and diffusion, coupled to the flow through Marangoni stresses.

Figure 12 takes cuts through the same parameter space along the $D_{\Gamma}$ direction at fixed $B$. The upper panel plots the critical field for Quincke rotation against $D_{\Gamma}$ for four surfactant strengths. At large $D_{\Gamma}$ all curves approach the clean-drop value because diffusion prevents appreciable surfactant gradients. Decreasing $D_{\Gamma}$ then separates the families: weak surfactant ($B=0.1$) still delays the onset relative to the rigid-cylinder value, whereas $B=0.3$–$3$ can drive it below $E_{c,s}$. The two middle panels show the finite rotation present at fixed electric fields $E_0/E_{c,s}=1.05$ and $1.1$. Remarkably, $D_{\Gamma}$ itself acts as a bifurcation parameter: at fixed electric field and surfactant strength $B$, reducing surfactant diffusion can trigger Quincke rotation. The lower panels show the corresponding area-averaged strain intensity. A state can therefore rotate strongly while its strain intensity remains small: rotation can begin before the Taylor circulation is fully restored.

Soluble exchange and the hidden-variable problem

All results above use an insoluble monolayer, with surface diffusion as its only explicit surfactant-relaxation process. The interpretation of surface-speed measurements nevertheless extends beyond this particular transport model. In the hidden-variable perspective of Manikantan and Squires [2], the same macroscopic response can arise from different combinations of elasticity, transport and surface rheology. The translating-disk analysis of Elfring et al. 52, for example, shows how Marangoni stresses can complicate the interpretation of interfacial drag as a surface-viscosity measurement. Li and Manikantan 53 demonstrated this ambiguity in the stability of liquid jets: when exchange is rapid compared with instability growth and surface diffusion is negligible, Marangoni stresses coupled to finite adsorption and desorption contribute to the dispersion relation in the same form as an intrinsic surface viscosity. They also showed that faster exchange and stronger surface diffusion both relax concentration gradients and weaken the Marangoni resistance to interfacial deformation. The present model has no intrinsic surface viscosity, and so tests diffusion versus exchange.

Dip coating provides useful examples of both the ambiguity and ways to constrain it. Scheid et al. 54 attributed residual film thickening under conditions of rapid surfactant supply to surface rheology. A thickness change is therefore not automatically evidence for a persistent Marangoni gradient. Delacotte et al. 55 varied withdrawal speed, concentration and surfactant chemistry, and connected changes in thickening to confinement of the surfactant supply and, for some systems, adsorption barriers. Neither finite-film reservoir depletion nor intrinsic surface viscosity is included in our local-exchange model. Conversely, Champagney et al. 51 compared free-film drawing and dip coating using a model in which surfactant exchange is negligible over the flow timescale, and obtained a consistent surface elasticity that also agreed with independent Langmuir-trough measurements. Their comparison illustrates how different forcing geometries and independent measurements of surface elasticity in a Langmuir trough can constrain an interpretation that film thickness alone leaves ambiguous.

Here the imposed electric field changes the coupled charge and surfactant distributions and can drive a symmetry-breaking rotational instability. This offers a complementary, spatially resolved setting for testing whether similar surface speeds and rotation thresholds conceal distinct surfactant transport mechanisms. To examine whether exchange produces changes similar to those caused by surface diffusion, we set $D_{\Gamma}=0$ and allow exchange with a well-mixed reservoir. Figure 13 shows the resulting branches for $M=1$ and $B=1$. When exchange is rapid ($K_s=40$ and $4$), concentration variations are continually relaxed and the rotating branch remains close to the clean-drop response, including a subcritical unstable segment connected to the stable rotating branch at a saddle-node bifurcation. As exchange is slowed, the critical field for Quincke rotation shifts to lower values and the onset becomes supercritical, as also found for weak surface diffusion in Figure 3. Movies S5 and S6 illustrate the development of Taylor circulation and Quincke rotation with rapid exchange ($K_s=4$) and no surface diffusion; their parameters are listed in the Supplementary Material.

The correspondence between diffusion and exchange is nevertheless incomplete. In particular, the subcritical branch at $K_s=4$ differs from the $D_{\Gamma}=1$ case, where rotation emerges supercritically and a much narrower bistable interval develops at finite rotation amplitude. The paired values $K_s=40$, $4$, $0.4$, $0.04$ and $0.004$ and $D_{\Gamma}=10$, $1$, $0.1$, $0.01$ and $0.001$, respectively, match the relaxation rate of the quadrupolar surfactant mode, as derived below. This matching gives the same leading-order response in the weak-field Taylor regime, but does not imply identical nonlinear bifurcation diagrams. Charge and surfactant transport are nonlinearly coupled through the interfacial flow, and diffusion and exchange relax the higher spatial harmonics at different rates. Increasing the viscosity ratio weakens the straining circulation and its feedback on the interfacial distributions, making the leading-order quadrupolar description more appropriate. We therefore use $M = 5$ for the subsequent comparisons in Figure 14 and Figure 15, where the two transport mechanisms produce closely similar responses.

Vector (PDF)

Vector (PDF)Vector (PDF)

To establish the leading-order correspondence, we first linearise the surfactant equation about the uniform state. The diffusion and local-exchange operators damp a Fourier perturbation $\delta\Gamma e^{\mathrm{i}n\theta}$ at the rate due to diffusion or exchange alone,

The comparison below keeps only one of these two transport mechanisms active. The Taylor circulation is dominated by the quadrupolar surfactant mode, so $\Lambda_2 = 4D_{\Gamma}$ for the diffusion-only family and $\Lambda_2 = K_s$ for the exchange-only family.

The collapse expected from $\Lambda_2$ can be made explicit in the weak-field Taylor regime. Retaining only the dominant quadrupolar mode, write

For small concentration variations about $\Gamma= 1$, we linearise the surfactant transport equation and the Langmuir equation of state while retaining finite $\chi$. The steady transport equation then gives

The Marangoni traction opposes the surface motion, with amplitude $-4BU_2/[(1-\chi)\Lambda_2]$, while the hydrodynamic traction jump has amplitude $-4(1+M)U_2$. At leading order in weak electric forcing, the electric traction is unchanged by surfactant. Comparing with the clean Taylor solution at the same field therefore gives

We define $\mathcal{N}_{\mathrm{imm}}$ as a dimensionless immobilisation parameter for this weak-field quadrupolar approximation. It measures the strength

of Marangoni suppression relative to the clean Taylor circulation: $\mathcal{N}_{\mathrm{imm}} \ll1$ gives weak suppression, whereas $\mathcal{N}_{\mathrm{imm}} \gg1$ gives strong suppression. (65) predicts the same leading-order surface-speed reduction for diffusion and exchange when their quadrupolar relaxation rates $\Lambda_2$ are equal. The derivation is given in the Supplementary Material. By matched relaxation we mean choosing paired diffusion-only and exchange-only cases with the same $\Lambda_2$, equivalently $K_s = 4D_{\Gamma}$ between the paired cases. These calculations use $M = 5$, $R = S = 2$, $D_q = 10^{-3}$, $B = 1$ and $\chi= 0.2$. The Taylor-regime speeds, concentration gradients and stresses are compared at the same pre-onset field $Re_{\mathrm{E}} = 1.0881$, corresponding to $E_0/E_{c,s} \simeq0.85$; the critical fields are determined separately by varying the applied field for each parameter set. The larger viscosity ratio is used here to compare the two surfactant-relaxation mechanisms away from the sharpest charge-density gradients at low $M$.

The two surfactant-transport mechanisms give almost the same straining surface speed and critical field for Quincke rotation in Figure 14. The straining surface speed, normalised by its clean value at the same field, nearly collapses for diffusion and exchange over more than three decades in $\Lambda_2$. The weak-field calculation in Eq. (65) shows that this agreement is expected for the dominant quadrupolar mode: at leading order the surface-speed ratio depends on the surfactant-transport parameters only through $\Lambda_2$. The corresponding critical fields also lie close together and approach the clean-drop value as $\Lambda_2$ increases. The nonlinear calculations show that the two mechanisms produce similar surface speeds and rotation thresholds even at finite field, while their surfactant distributions differ.

The associated surfactant gradients and Marangoni stresses are shown in Figure 15. Both $\lVert\partial_{\theta}\Gamma\rVert_{\mathrm{rms}}$ and $\lVert\tau_M\rVert_{\mathrm{rms}}$ decrease with $\Lambda_2$, following broadly similar trends for diffusion and exchange. At small $\Lambda_2$, the Taylor circulation produces large concentration gradients and a strong Marangoni stress. Either faster surface diffusion or faster exchange with the reservoir reduces those gradients, weakens the Marangoni stress and remobilises the interface. At the level of the dominant mode, it matters little whether the decay is produced by surface diffusion or by local exchange with the bulk reservoir.

The surfactant distributions are not the same, however. Surface diffusion damps mode $n$ at rate $n^2D_{\Gamma}$, whereas local exchange with the well-mixed reservoir damps every mode at $K_s$. At equal $\Lambda_2$, the diffusion and exchange calculations give nearly identical straining surface speeds and critical fields, but different amplitudes of the higher Fourier modes of surfactant concentration. At $\Lambda_2 = 1$, for example, the computed Taylor surface-speed ratios are 0.8326 and 0.8332 for diffusion and exchange, respectively, while $|\Gamma_4|/|\Gamma_2|$ is 0.0743 and 0.2326, a difference of more than a factor of three (Figure 16). Surface speed and rotation threshold alone may therefore be insufficient to distinguish diffusion from exchange. The angular variation of surfactant concentration provides an additional measure that distinguishes the two mechanisms in these calculations.

Vector (PDF)

The distinction arises because matching the relaxation of the $n = 2$ mode does not match that of the higher modes. For example, setting $K_s = 4D_\Gamma$ makes the two rates equal for $n = 2$, but diffusion damps $n = 4$ four times faster than the paired exchange process. The ratio $|\Gamma_4|/|\Gamma_2|$ therefore distinguishes the calculated concentration distributions even when their surface speeds are very similar. Resolved concentration measurements could provide this additional information. The prescribed reservoir assumes that bulk transport maintains a spatially uniform concentration adjacent to the interface. Resolving finite bulk transport would relax this assumption and could modify the comparison between exchange and surface diffusion.

Conclusions

In this paper we studied a surfactant-covered leaky-dielectric drop in a uniform dc electric field using a fixed circular interface. The model retains the full nonlinear surface-charge convection together with surfactant advection, surface diffusion, Marangoni stresses and, separately, local exchange with a well-mixed bulk reservoir. This setting resolves the charge and surfactant fields around the entire interface without a small-$\mathrm{Re}_E$ expansion or a truncation to the first few harmonics.

The problem couples two advected interfacial fields. Surface charge sets the electric shear stress and is convected by the flow, while the same flow redistributes surfactant and generates a Marangoni stress. In the non-rotating Taylor regime, our results recover the Marangoni immobilisation mechanism previously identified for a non-diffusing surfactant: as surface diffusion decreases, the Marangoni stress carries most of the tangential electric-stress balance and the hydrodynamic stress and Taylor circulation become small. With finite surface diffusion the interface does not become exactly stagnant, because diffusion continually smooths the surfactant gradient, but it can become nearly immobilised. Retaining charge convection at arbitrary strength resolves the sharp charge-density gradients that develop in the Taylor regime and determine the base state from which Quincke rotation begins.

When charge convection destabilises the nearly immobilised Taylor state, Quincke rotation restores both interfacial motion and a substantial straining circulation. At weak surfactant diffusion, this recovery occurs despite persistent concentration gradients: electric and Marangoni tractions no longer nearly cancel. The circulation also changes orientation. For the clean drop and larger surfactant diffusivities, its central extensional axis has positive tilt on the positive-rotation branches. At low diffusivity, the tilt changes from negative to positive as the field increases, while the direction of rotation remains unchanged.

Surfactant also changes the Quincke rotational instability itself. It can lower the critical field below the rigid-cylinder value, produce a non-monotonic dependence of the critical field on surfactant strength, and change the onset from subcritical to supercritical. At fixed electric field and surfactant strength, the surface diffusivity can itself act as a bifurcation parameter, with decreasing $D_\Gamma$ driving the drop from the Taylor regime into Quincke rotation. The Taylor-regime base state and the rotational instability are therefore coupled through the redistributed surfactant and charge fields.

Surface diffusion and exchange with a well-mixed bulk reservoir act through different physical processes. In the $M = 5$ comparison, the straining surface speed and critical field for Quincke rotation nearly coincide when the dominant quadrupolar surfactant mode has the same relaxation rate in the two systems. A weak-field calculation predicts this agreement in surface speed. The nonlinear calculations retain higher spatial modes, whose amplitudes differ because diffusion and exchange damp them at different rates. Surface speed and rotation threshold alone can therefore conceal the transport mechanism, whereas the spatial variation of surfactant concentration supplies additional information.

The circular geometry suppresses shape changes but allows charge convection at arbitrary strength and resolves the sharp interfacial charge and surfactant gradients that make fully resolved three-dimensional continuation expensive, particularly at low viscosity ratios. The results are therefore most directly applicable when surface tension is sufficiently strong that deformation remains small. Extending the present diffusion–exchange comparison and rotation-induced remobilisation mechanism to a fully resolved deforming three-dimensional interface, or to a finite bulk surfactant reservoir, would require substantially greater computational effort.

Another important extension is to determine how interactions between drops modify the immobilisation and rotation-induced remobilisation identified here. Surfactant can weaken electrohydrodynamic interactions and increase the relative importance of dielectrophoresis when charge convection is negligible. 26 The competition between these interactions also governs drop trajectories and collective organisation in emulsions. 56,57 More recently, a reduced model retaining charge convection has captured interacting clean drops undergoing Quincke rotation in the regime of large separation, high viscosity ratio and small deformation.58 These studies raise the question of how surfactant modifies interactions when charge convection by the straining flow becomes appreciable, including whether rotation restores interactions weakened by Marangoni stresses. Addressing this question with spatially resolved nonlinear charge and surfactant transport would be computationally demanding, requiring sharp interfacial gradients to be resolved on multiple moving and deforming drops together with their coupled electric and hydrodynamic interactions.

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary Material is available. Raw numerical data are available from the corresponding author upon reasonable request.

Acknowledgements

This work was supported by the Leverhulme Trust (Research Project Grant RPG-2026-185).

We thank Harishankar Manikantan for helpful discussions and valuable insights into the physical interpretation of the problem and Alexander Wray for constructive comments on the manuscript.

Numerical method

Notation and nondimensionalisation in this section are as in the main text. For a circular interface, the bulk Laplace and Stokes problems can be solved analytically to express the interfacial electric field and velocity in terms of the surface charge and surfactant concentration. The resulting boundary-integral formulation uses linear operators that are diagonal in Fourier space and can be expressed using the periodic Hilbert transform. We follow the interfacial formulation used for the clean circular-drop problem by Peng and Schnitzer [1], express it in terms of $R$, $S$, $M$ and $\mathrm{Re}_{E}$, and include surfactant transport and Marangoni stresses.

Boundary operators and interfacial evolution equations

Let $\mathcal{H}$ denote the circular Hilbert transform,

with the convention $\mathcal{H}[\cos(n\theta)]=\sin(n\theta)$. Eliminating the harmonic electric potentials gives the tangential electric field at the interface as

Similarly, eliminating the two Stokes fields gives a relation between the derivative of the surface velocity and the total non-hydrodynamic tangential traction,

where

The mean surface velocity is fixed by the electric–hydrodynamic torque balance,

Equations (S2)–(S5) therefore reconstruct the electric shear and interfacial velocity directly from $q(\theta,t)$ and $\Gamma(\theta,t)$. Substituting the interfacial normal electric fields into the charge-conservation equation gives

The forcing term changes sign with $RS-1=\tau^{-}/\tau^{+}-1$, which is also the material distinction relevant to the direction of the Taylor circulation. The surfactant equation is

The steady problem is therefore obtained by setting $\partial_t q = \partial_t \Gamma= 0$ in Eqs. (71) and (72), with $E_\theta$, $\tau_M$ and $u_\theta$ determined from Eqs. (67)–(70). The clean-drop formulation is recovered by suppressing the Marangoni term and omitting Eq. (72).

Periodic discretisation and reconstruction

We discretise $0 \le\theta< 2\pi$ on an even periodic grid

On a uniform grid the circular Hilbert transform is represented by its exact trigonometric collocation matrix,

while first and second derivatives in the transport equations are evaluated with periodic centred differences,

Indices are periodic. We combine centred differences for conservative interfacial transport with Fourier evaluation of the circular electric and hydrodynamic boundary operators on a uniform periodic grid. The Hilbert operator is spectrally exact on the retained trigonometric modes, while the transport discretisation determines the overall second-order spatial accuracy for smooth solutions. The same spatial discretisation is used for steady continuation, transient integration and linear stability. Spatial accuracy is assessed by successive grid refinement of the interfacial fields and their gradients.

In the implementation, the Hilbert transform is evaluated by its Fourier multiplier $-i\,\operatorname{sgn}(n)$, with the mean and Nyquist coefficients set to zero. The Nyquist mode is the highest grid wavenumber, $n = N/2$. The Marangoni traction is evaluated by applying the discrete derivative directly to the Langmuir expression for $(\gamma-1)/\mathrm{Ca}_E$,

This expression has zero discrete mean on the periodic grid and converges to Eq. (69). At finite spacing it differs from applying a pointwise chain rule to $D_1\Gamma$. Its analytic variation, $D_1[-\mathrm{Ma}\,\delta\Gamma/(1-\chi\Gamma)]$, is used in the discrete Jacobian. The velocity is reconstructed directly in Fourier space. For each Fourier mode with $0 < |n| < N/2$, the mobility multiplier is $[2(1+M)|n|]^{-1}$ acting on $\tau_E+\tau_M$; the mean is supplied by Eq. (70). This is equivalent to spectral integration of Eq. (68). The Nyquist velocity coefficient is set to zero. The steady residuals are

The local exchange term therefore adds no new nonlocal operator: it enters directly in the surfactant residual and Jacobian. The steady Taylor- and Quincke-regime solutions considered here obey the half-turn symmetries

These symmetries are imposed in the nonlinear solve. For an $N$-point grid only $N/2$ charge values are independent. For $K_s=0$, the prescribed mean $\langle\Gamma\rangle=1$ removes the conserved surfactant-mass mode; for $K_s>0$ the same mean is fixed by the reservoir equilibrium.

The branch and stability calculations use uniform grids with $N=128$, 256 and 512, with refined clean and exchange-only branches using $N=1024$ and 2048. Selected states on the clean and $K_s=40$ branches are additionally corrected at $N=4096$ where required, using the same acceptance tolerances. Selected interfacial profiles were also recomputed on finer grids. Both low-diffusion states in Figure 6 of the main text use $N=512$, whereas the clean Taylor states at $E_0/E_{c,s}=1.30$ and 1.80 in Figure 2 use $N=4096$. Narrow charge layers near convergent stagnation points make convergence particularly demanding at low viscosity ratio and weak charge diffusion.

Newton iteration and pseudo-arclength continuation

The non-rotating Taylor branch is initiated from the weak-field solution, for which $q \propto\cos\theta$ and $\Gamma= 1$, and continued in $\operatorname{Re}_{E}$ using the converged solution at the previous parameter value as the initial guess. At fixed $\operatorname{Re}_{E}$ the nonlinear system $\boldsymbol{F}(\boldsymbol{x},\operatorname{Re}_{E}) = \boldsymbol{0}$ is solved by Newton iteration with an analytically differentiated discrete Jacobian. A backtracking line search is applied whenever the full Newton step does not reduce the residual. Continuation using $\operatorname{Re}_{E}$ as the parameter cannot pass a fold, where the steady solution ceases to be a single-valued function of $\operatorname{Re}_{E}$. Rotating branches are therefore followed by pseudo-arclength continuation in the state and logarithmic field parameter $\rho= \ln\operatorname{Re}_{E}$. With $n_x$ independent real state coordinates, the metric is

Two converged rotating states define a unit secant tangent $(\boldsymbol{t}_x,t_\rho)$ in this metric. A predictor $(\boldsymbol{x}_p,\rho_p)$ is corrected by solving

The rotating branch is initiated using the symmetry-breaking eigenmode at the Taylor-state instability. Two nearby converged solutions with small nonzero rotation provide the initial continuation direction. The continuation step is adjusted automatically to maintain convergence. Accepted solutions satisfy the discrete steady equations and the concentration constraints $\Gamma\ge0$ and $\chi\Gamma< 1$. Pseudo-arclength continuation allows both stable and unstable segments to be followed through folds.

Parameter sweeps use the previously converged solution as the initial guess at the next parameter value. The rotation maps are constructed from steady solution branches, with stability determined by the eigenvalue analysis described below.

Transient calculations and tangential stresses

Transient calculations use the same spatial discretisation as the steady calculations and the dimensionless time $t = t_{\mathrm{dim}}/\tau^{+}$ used in the main text. The histories in Figure 9 describe relaxation in the non-rotating Taylor regime following application of a constant electric field at $\operatorname{Re}_{E} = 1$.

Time integration uses a second-order exponential Runge–Kutta scheme (ETD2RK){sup}2`. The linear operator $L$ includes charge relaxation, diffusion, exchange and the linearised Marangoni contribution about uniform concentration, and is diagonal in Fourier space. Its exponential is evaluated directly, avoiding the stability restrictions associated with explicit treatment of these linear terms. The remaining terms are retained in $\mathcal{N}$, so the splitting does not alter the nonlinear equations.

Writing the Fourier-space system as $\mathrm{d}\hat{\boldsymbol{y}}/\mathrm{d}t = L\hat{\boldsymbol{y}}+\mathcal{N}(\boldsymbol{y})$, where $\mathcal{N}(\boldsymbol{y})$ denotes the Fourier coefficients of the remaining terms, a step of length $\Delta t$ is

where $\varphi_1(z)=(e^z-1)/z$ and $\varphi_2(z)=(e^z-1-z)/z^2$, evaluated by their regular series near zero. Electric fields, tractions and velocity are reconstructed at each stage. One full step is compared with two half steps, using their difference divided by three as the local error estimate. The two-half-step state is accepted when the error and concentration-admissibility tests pass. The time step also satisfies an advective restriction.

The three tangential tractions plotted in the main text are evaluated directly on the interface,

The Marangoni traction is discretised using Eq. (76). The hydrodynamic traction is evaluated from the Fourier Stokes operator, with multiplier $-2(1+M)|n|$ for nonzero modes and $-2$ for the mean. The final equality in Eq. (85) supplies a check of the tangential balance on resolved modes.

Linear stability and numerical checks

The stability routines express growth rates using the electroviscous time $\hat{t}=\operatorname{Re}_{E}t$. A computed growth rate $\hat{\lambda}$ therefore corresponds to $\lambda=\operatorname{Re}_{E}\hat{\lambda}$ in the main-text time variable. Since $\operatorname{Re}_{E}>0$, this conversion leaves the stability classification unchanged.

For linear stability we linearise the semi-discrete time-dependent system represented by Eqs. (71) and (72) about each steady state. For $K_{s}=0$ the conserved surfactant mean defines a neutral mass mode, which is removed by restricting perturbations to zero mean; for $K_{s}>0$ the corresponding mean perturbation decays at rate $K_{s}$ in the time variable $t$, or $K_{s}/\operatorname{Re}_{E}$ in $\hat{t}$, and is not neutral. The steady solve and stability calculation both use the half-turn symmetry in Eq. (79). The perturbations satisfy $\delta q(\theta+\pi)=-\delta q(\theta)$ and $\delta\Gamma(\theta+\pi)=\delta\Gamma(\theta)$, with both sine and cosine components included. The reported stability classifications therefore apply only to this symmetry class.

A state is classified as linearly stable to perturbations preserving the imposed half-turn symmetries when the largest real part of the computed eigenvalues, in the time variable $\hat{t}$, is below $-10^{-7}$, unstable when it is above $10^{-7}$, and marginal otherwise. The onset of Quincke rotation on the Taylor branch is located from the zero crossing of the leading real eigenvalue. Along a rotating branch, changes between solid and dotted curves in the figures use the same criterion. Here bistability denotes coexistence of stable states with different rotation magnitudes over the same field interval. These may be a Taylor state and a rotating state, or two rotating states with different amplitudes. Figure S8 illustrates the latter case through a zoom of the rotating branch.

Vector (PDF)

Vector (PDF)

Several numerical checks are used throughout. The analytic steady Jacobian is compared with centred finite-difference directional derivatives, and the stability spectrum is checked against a direct finite-difference Jacobian for representative states. The clean-interface limit reproduces the circular-drop branches and stability changes of Peng and Schnitzer [1]; the weak-field solution is recovered as $\operatorname{Re}_{E}\to0$; and the surfactant solver approaches the clean response when Marangoni coupling is removed or when diffusion or exchange becomes sufficiently rapid.

Plots with a fixed reference velocity

Figures S1–S7 reproduce Figure 2, Figure 3, Figure 5, Figure 8, Figure 11, Figure 12 and Figure 13 of the main text using the fixed velocity scale

Thus velocities and rotation rates are multiplied by $e^{2}$, and the strain intensity by $e^{4}$:

The underlying solutions, stability classifications, critical fields and flow-type parameter $\Theta$ are unchanged. Figure 6 of the main text already uses $U_{c,s}$. The straining angle in Figure 7 is invariant under this change of velocity scale and is not repeated.

Illustrative dimensional realisation

To give dimensional context for the parameters used in this study, we use the castor-oil suspending-phase and silicone-oil material properties reported in the Quincke experiments of Salipante and Vlahovska [3] and take an illustrative drop radius $a=0.9\,\mathrm{mm}$. Specifically, we use $\varepsilon^{+}/\varepsilon_{0}=5.3$, $\sigma^{+}=4.5\times10^{-11}\,\mathrm{S\,m^{-1}}$, $\mu^{+}=0.69\,\mathrm{Pa\,s}$ and $T=298\,\mathrm{K}$. These values give $\tau^{+}\simeq1.04\,\mathrm{s}$ and $\operatorname{Re}_{E}=1$ at $E_{0}\simeq1.19\times10^{5}\,\mathrm{V\,m^{-1}}=119\,\mathrm{V\,mm^{-1}}$. Table S1 gives one dimensional realisation of selected parameters for these specified properties; the mapping is not unique when

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

the other dimensional quantities can vary. For example, $B = 1$ corresponds to an equilibrium surface excess of about $0.24\,\mu\mathrm{mol}\,\mathrm{m}^{-2}$, while $D_{\Gamma} = 10^{-3}$ corresponds to a surface diffusivity of order $10^{-9}\,\mathrm{m}^{2}\,\mathrm{s}^{-1}$.

The conversions for the transport and concentration parameters are

They apply after the remaining dimensional properties have been specified. The illustrative outer-fluid properties do not fix the independently chosen electrical ratios $R$ and $S$ of the model.

Weak-field derivation of the surface speed in Taylor regime

We derive the leading-order suppression of Taylor circulation due to surfactant. Surface-charge convection is weak in this limit, so the leading electric traction is unchanged by surfactant. Writing $\Gamma= 1+g$, we assume $|g| \ll1$ and $|\chi g| \ll1-\chi$, while retaining finite $\chi$. The steady surfactant equation is

Write the leading quadrupolar fields as

Linearisation about $\Gamma= 1$ gives

The main-text comparisons use either diffusion alone or exchange alone. The linearised Langmuir traction is

The hydrodynamic coefficient follows from the Stokes solution inside and outside a circular interface. Using $u_{r}=r^{-1}\partial_{\theta}\psi$ and $u_{\theta}=-\partial_{r}\psi$, the streamfunctions that are regular inside, decay in velocity outside and satisfy $u_{r}(1)=0$ are

Thus

At $r=1$, $\partial_{\theta}u_{r}=0$ and the shear stress is proportional to $\partial_{r}u_{\theta}-u_{\theta}/r$. In units of $\epsilon^{+}E_{0}^{2}$, this gives

Let $\tau_{E}=T_{2}\sin2\theta$ be the leading electric traction. Its amplitude is the same in the covered and clean systems at this order because charge convection enters at higher order. The balance $\tau_{H}+\tau_{E}+\tau_{M}=0$ gives

For the clean interface, $T_{2}=4(1+M)U_{2,\mathrm{clean} }$. Eliminating $T_{2}$ therefore yields

For the quadrupolar field, $U_{\mathrm{str}}=|U_{2}|/\sqrt{2}$, so the same ratio applies to the RMS straining surface speed. Matching $\Lambda_{2}$ therefore gives the same leading-order surface-speed ratio for diffusion and exchange. At finite electric field, changes in charge convection alter the electric traction, and higher harmonics of velocity and concentration contribute. Equation (97) is compared with the full numerical results in Figure 14(a) of the main text.

Supplementary movies

Movies S1–S6 show the evolution of surface charge, surfactant concentration where present, and the interior flow following application of a constant electric field. The common parameters are $R = S = 2$, $M = 1$ and $D_q = 10^{-3}$. The surfactant-covered cases have $B = 1$ and $\chi= 0.2$; the clean cases have no surfactant and $B = 0$. The applied field, transport parameters and number of uniform angular grid points are listed in Table S2. The field reference is the same rigid-cylinder value as in the main text, with $Re_{E,c,s} = 1.5060$ and $Re_E = 1.5060(E_0/E_{c,s})^2$.

Movie	Interface	Regime	$E_0/E_{c,s}$	$D_\Gamma$	$K_s$	$N$
S1	Clean	Taylor	1.50	–	–	2048
S2	Clean	Quincke	2.00	–	–	2048
S3	Insoluble surfactant	Taylor	0.85	$10^{-3}$	0	1024
S4	Insoluble surfactant	Quincke	1.00	$10^{-3}$	0	1024
S5	Soluble surfactant	Taylor	0.85	0	4	2048
S6	Soluble surfactant	Quincke	2.00	0	4	4096

All movies begin with an uncharged interface, $q(\theta,0) = 0$. For the surfactant-covered cases, $\Gamma(\theta,0) = 1$. Movies S1, S3 and S5 follow relaxation to a non-rotating Taylor state. In Movies S2 and S6, a small perturbation $\delta q = -10^{-3}\sin\theta$ is applied at $t = 6$ to initiate the rotational instability. Movie S4 first approaches the Taylor state while reflection symmetry is imposed. This reflection constraint is then removed and a perturbation $\delta q = -10^{-4}\sin\theta$ initiates rotation. The applied electric field and all material parameters remain unchanged during each movie, including at the perturbation. The perturbations preserve zero net charge and do not change the surfactant concentration.

The displayed time is $t = t_{\mathrm{dim}}/\tau^+$, labelled $T$ in the movies. Surface speeds are scaled by $U_E$, and angular velocities by $U_E/a$. The interior velocity is reconstructed from the interfacial solution at each time. Interior tracers follow this velocity. The charge symbols represent weighted surface-charge patches whose weights change through Ohmic currents and diffusive transport. For the insoluble surfactant, surface markers represent the combined advective and diffusive surfactant flux. In the exchange movies, their weights also change through the local source $K_s(1-\Gamma)$.

Movies S3 and S5 compare weak surface diffusion with rapid exchange at the same field. Movies S2 and S6 compare clean and rapidly exchanging interfaces at the same field. These movie pairs illustrate the different transport regimes; the matched-relaxation comparisons are those in Figure 14–Figure 16 of the main text.