1. Introduction
Electric fields are widely used to steer particles and droplets for applications in directed assembly (van Blaaderen et al. Reference van Blaaderen, Dijkstra, van Roij, Imhof, Kamp, Kwaadgras, Vissers and Liu2013; Harraq, Choudhury & Bharti Reference Harraq, Choudhury and Bharti2022), microfluidics (Link et al. Reference Link, Grasland-Mongrain, Duri, Sarrazin, Cheng, Cristobal, Marquez and Weitz2006; Hartmann, Schuer & Hardt Reference Hartmann, Schuer and Hardt2022), ink-jet printing (Basaran, Gao & Bhat Reference Basaran, Gao and Bhat2013), modulation of emulsion microstucture and rheology (Eow & Ghadiri Reference Eow and Ghadiri2002; Tao et al. Reference Tao, Tang, Tawhid-Al-Islam, Du and Kim2016) and electrosprays (Ganan-Calvo et al. Reference Ganan-Calvo, Lopez-Herrera, Herrada, Ramos and Montanero2018). An important issue in practical applications is the droplet interactions due to electric polarization and electrohydrodynamic flows. In the canonical case of an applied uniform electric field, the induced dipoles promote particle chaining along the applied field direction (Zukoski Reference Zukoski1993; Sheng & Wen Reference Sheng and Wen2012). In addition to the electrostatic interactions, particles may interact electrohydrodynamically due to induced-charge electrophoretic flows in the case of ideally polarizable particles (Squires & Bazant Reference Squires and Bazant2004) or electric-shear-driven flows around droplets (Melcher & Taylor Reference Melcher and Taylor1969). These flows can be either cooperative or antagonistic to the dipolar interactions (Baygents, Rivette & Stone Reference Baygents, Rivette and Stone1998; Saintillan Reference Saintillan2008; Park & Saintillan Reference Park and Saintillan2010; Sorgentone et al. Reference Sorgentone, Kach, Khair, Walker and Vlahovska2021) and prevent chaining (Ha & Yang Reference Ha and Yang2000). Recently, the three-dimensional interactions of a pair of identical droplets were investigated by means of numerical simulations using the boundary integral method, asymptotic theory for large separations and spherical droplets (Sorgentone, Tornberg & Vlahovska Reference Sorgentone, Tornberg and Vlahovska2019; Sorgentone et al. Reference Sorgentone, Kach, Khair, Walker and Vlahovska2021; Sorgentone & Vlahovska Reference Sorgentone and Vlahovska2021) and experiments (Kach, Walker & Khair Reference Kach, Walker and Khair2022). The systematic exploration of the effects of fluid properties and the droplet initial configuration revealed intricate relative motions that eventually lead to either droplet coalescence or indefinite repulsion; only if the droplets line-of-centres were initially perpendicular to the applied field direction and the electrohydrodynamic flow along the droplet surface were equator-to-pole, the drops’ motion is eventually arrested and the drops remain at an equilibrium separation.
Asymmetry in terms of droplet size or properties is expected to increase the complexity of the droplet interactions; however, the problem has been studied only to a limited extent for small droplet deformations (Kach et al. Reference Kach, Walker and Khair2022) or only configurations where droplets are aligned with the field (Zabarankin Reference Zabarankin2020). Here, we analyse the three-dimensional interactions of dissimilar drops using both theory and simulations. We focus on droplets with different electrical properties (conductivity and permittivity). We find non-trivial dynamics, such as droplets ‘dancing’, where droplets execute complex trajectories before coming into contact or separating, or ‘swimming’, where droplets reach stable separation and the pair translates in a direction either parallel or perpendicular to the applied field. The intricate behaviours are linked to non-reciprocal interactions generated by the electrohydrodynamic flows.
2. Problem formulation
Let us consider two neutrally buoyant and charge-free drops with radii 
$a_i$ and different viscosities 
$\eta _{{d},i}$, conductivities 
$\sigma _{{d},i}$ and permittivities 
${\varepsilon }_{{d},i}$, suspended in a fluid with viscosity 
$\eta _{{m}}$, conductivity 
$\sigma _{{m}}$ and permittivity 
${\varepsilon }_{{m}}$. The mismatch of drop and suspending fluid properties is characterized by the conductivity, permittivity and viscosity ratios:
\begin{equation} R_i=\frac{\sigma_{{d},i}}{\sigma_{{m}}},\quad S_i=\frac{{\varepsilon}_{{d},i}}{{\varepsilon}_{{m}}},\quad \lambda_i=\frac{\eta_{{d},i}}{\eta_{m}}, \quad i=1,2. \end{equation}
The difference in drop size introduces one more parameter, 
$\nu =a_2/a_1$. The distance between the drops’ centroids is 
$d$ and the angle between the drops’ line-of-centres with the applied field direction is 
$\varTheta$. The unit separation vector between the drops is defined by the difference between the position vectors of the drops’ centres of mass 
$\hat {\boldsymbol {d}}=({\boldsymbol {x}}_2^c-{\boldsymbol {x}}^c_1)/d$. The unit vector normal to the drops’ line-of-centres and orthogonal to 
$\hat {\boldsymbol {d}}$ is 
$\hat {\boldsymbol {t}}$. The problem geometry is sketched in figure 1.
Figure 1. Two initially spherical drops with different radii 
$a_i$, fluid viscosities 
$\eta _{{d}, i}$, permittivities 
${\varepsilon }_{{d},i}$ and conductivities 
$\sigma _{{d},i}$ (
$i=1,2$) suspended in a fluid with viscosity 
$\eta _{m}$, permittivity 
${\varepsilon }_{m}$ and conductivity 
$\sigma _{m}$ and subjected to a uniform DC electric field 
${\boldsymbol {E}}^\infty =E_0{\boldsymbol {\hat z}}$. The angle between the line-of-centres vector and the field direction is 
$\varTheta =\arccos ({\boldsymbol {\hat z}}{\boldsymbol{\cdot}} \hat{\boldsymbol {d}})$.
We adopt the leaky dielectric model, which is widely used to describe the electrohydrodynamics of weakly conducting, viscous fluids (Melcher & Taylor Reference Melcher and Taylor1969; Saville Reference Saville1997; Vlahovska Reference Vlahovska2019). Fluid motion is described by Stokes equations:
\begin{equation} \eta \nabla^2{{\boldsymbol u}}-\boldsymbol{\nabla} p=0,\quad \boldsymbol{\nabla}\boldsymbol{\cdot}{{\boldsymbol u}}=0, \end{equation}
where 
${{\boldsymbol u}}$ and 
$p$ are the fluid velocity and pressure. The electric field 
${\boldsymbol {E}}$ obeys
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{E}}=0,\quad \boldsymbol{\nabla}\times{\boldsymbol{E}}=0. \end{equation}
Far away from the drops, 
${\boldsymbol {E}}\rightarrow {\boldsymbol {E}}^\infty =E_0{\boldsymbol {\hat z}}$ and 
${{\boldsymbol u}}\rightarrow 0$.
At the drop interfaces, normal electric current is continuous (assuming negligible charge convection, as originally proposed by Taylor Reference Taylor1966), 
$E_n^{m}={\it R}_i E_n^{d}$, where 
$E_n={\boldsymbol {E}}\boldsymbol {\cdot}\boldsymbol {n}$, and 
$\boldsymbol {n}$ is the outward pointing normal vector to the drop interface. The surface charge density adjusts to satisfy the current balance, leading to a discontinuity of the displacement field 
${\varepsilon }_{m}(E_n^{m}-{\it S}_i E_n^{d})=q$. The conductivity and permittivity ratios, 
${\it R}_i$ and 
${\it S}_i$, are defined by (2.1a–c).
The electric field acting on the induced surface charge 
$q$ gives rise to electric shear stress at the interface. The tangential stress balance yields
\begin{equation} \left({\boldsymbol{I}}-\boldsymbol{n}\boldsymbol{n}\right)\boldsymbol{\cdot} ( {\boldsymbol{T}}^{m}- {\boldsymbol{T}}^{d}) \boldsymbol{\cdot}\boldsymbol{n}+q{\boldsymbol{E}}_t=0,\quad {\boldsymbol{x}}\in \mathcal{D}, \end{equation}
where 
$T_{ij}=-p\delta _{ij}+\eta (\partial _j u_i+\partial _i u_j)$ is the hydrodynamic stress and 
$\delta _{ij}$ is the Kronecker delta function. The electric tractions are calculated from the Maxwell stress tensor 
$T^{el}_{ij}={\varepsilon } (E_iE_j-E_kE_k\delta _{ij}/2)$. Here 
${\boldsymbol {E}}_t={\boldsymbol {E}}-E_n\boldsymbol {n}$ is the tangential component of the electric field, which is continuous across the interface, and 
${\boldsymbol {I}}$ is the idemfactor. The normal stress balance is
\begin{equation} \boldsymbol{n} \boldsymbol{\cdot}({\boldsymbol{T}}^{m}- {\boldsymbol{T}}^{d}) \boldsymbol{\cdot}\boldsymbol{n}+\tfrac{1}{2}((E_n^{{m}})^2-S (E_n^{{d}})^2-(1-S)E_t^2)=\gamma \boldsymbol{\nabla}_s\boldsymbol{\cdot} \boldsymbol{n},\quad {\boldsymbol{x}}\in \mathcal{D}, \end{equation}
where 
$\gamma$ is the interfacial tension.
Henceforth, all variables are non-dimensionalized using the radius of the undeformed drop 
$a_1$, the undisturbed field strength 
$E_0$, a characteristic applied stress 
$\tau _c={\varepsilon }_{m} E_0^2$, and the properties of the suspending fluid. Accordingly, the time scale is 
$t_c=\eta _{m}/\tau _c$ and the velocity scale is 
$u_c=a_1\tau _c/\eta _{m}$. The ratio of the magnitude of the electric stresses and surface tension defines the electric capillary number 
${\it Ca}_i={{\varepsilon }_{m} E_0^2 a_i}/{\gamma }$.
3. Methodology
We utilize a boundary integral method for the numerical simulations and analytical results derived from an asymptotic theory for large separations and small deformations. These approaches were presented and validated for identical drops in Sorgentone et al. (Reference Sorgentone, Kach, Khair, Walker and Vlahovska2021). Here we summarize the extension of the theory and simulations to treat dissimilar drops.
3.1. Integral representation for the electric field and fluid velocity
The boundary integral formulation taking into account the fact that the two drops may have different permittivities and conductivities is
\begin{align} & {\boldsymbol{E}}^\infty({\boldsymbol{x}})+\sum_{i=1}^2 \int_{{\mathcal{D}}_i} \frac{\hat{{\boldsymbol{x}}}}{4{\rm \pi} r^3} {({\boldsymbol{E}}^{m}({\boldsymbol{y}})-{\boldsymbol{E}}^{{d}, i}({\boldsymbol{y}})) \boldsymbol{\cdot}\boldsymbol{n}({\boldsymbol{y}})}\,{\rm d} A({\boldsymbol{y}}) \nonumber\\ &\quad = \begin{cases} {\boldsymbol{E}}^{{d}, i}({\boldsymbol{x}}) & {\boldsymbol{x}} $ inside $ {\mathcal{D}}_i, \\ \dfrac{1}{2} ({\boldsymbol{E}}^{{d}, i}({\boldsymbol{x}})+{\boldsymbol{E}}^{m}({\boldsymbol{x}})) & {\boldsymbol{x}} \in{\mathcal{D}}_i, \\ {\boldsymbol{E}}^{m}({\boldsymbol{x}}) & {\boldsymbol{x}} \mbox{ outside } {\mathcal{D}}_i, \end{cases} \end{align}
where 
$i=1,2$ denotes drop 
$i$, 
$\hat {\boldsymbol {x}}={\boldsymbol {x}}-{\boldsymbol {y}}$ and 
$r=|\hat {\boldsymbol {x}}|$ and 
${\rm d} A$ is the surface area element. Accordingly, the normal and tangential components of the electric field at the drop interface, 
${\boldsymbol {x}} \in {\mathcal {D}}_i$, are calculated as
\begin{gather} \frac{(R_i+1)}{2R_i}E^{m}_n({\boldsymbol{x}})={\boldsymbol{E}}^\infty({\boldsymbol{x}}) \boldsymbol{\cdot} \boldsymbol{n}({\boldsymbol{x}})+\sum_{j=1}^2\frac{R_j-1}{R_j}\boldsymbol{n}({\boldsymbol{x}})\boldsymbol{\cdot} \int_{{\mathcal{D}}_j}\frac{\hat{\boldsymbol{x}} }{4{\rm \pi} r^3} E^{m}_n({\boldsymbol{y}})\,{\rm d} A({\boldsymbol{y}}), \end{gather}
\begin{gather}{\boldsymbol{E}}_t({\boldsymbol{x}})={\boldsymbol{E}}^\infty({\boldsymbol{x}})-\frac{1+R_i}{2R_i}E^{m}_n({\boldsymbol{x}})\boldsymbol{n}({\boldsymbol{x}})+\sum_{j=1}^2 \int_{{\mathcal{D}}_j} \frac{\hat{{\boldsymbol{x}}}}{4{\rm \pi} r^3}{\left(\frac{R_j-1}{R_j}\right)} E^{m}_n({\boldsymbol{y}})\,{\rm d} A({\boldsymbol{y}}). \end{gather}For the flow field, we have developed the method for fluids of arbitrary viscosity, but for the sake of brevity here we list the equations in the case of equiviscous drops and suspending fluids. The velocity is given by
\begin{equation} 2{{\boldsymbol u}}({\boldsymbol{x}})={-}\sum_{j=1}^2 \left( \frac{1}{4{\rm \pi}}\int_{{\mathcal{D}}_j} \left(\frac{{\boldsymbol{f}}({\boldsymbol{y}})}{{Ca}_j}-{\boldsymbol{f}}^E({\boldsymbol{y}})\right)\boldsymbol{\cdot} \left(\frac{{\boldsymbol{I}}}{r}+\frac{\hat{\boldsymbol{x}}\hat{\boldsymbol{x}}}{r^3} \right){\rm d} A({\boldsymbol{y}})\right), \end{equation}
where 
${\boldsymbol {f}}$ and 
${\boldsymbol {f}}^E$ are the interfacial stresses due to surface tension and electric field:
\begin{equation} \left.\begin{gathered} {\boldsymbol{f}}=\boldsymbol{n}\boldsymbol{\nabla}_s\boldsymbol{\cdot} \boldsymbol{n},\quad {\boldsymbol{f}}^E=({\boldsymbol{E}}^{m}\boldsymbol{\cdot} \boldsymbol{n}){\boldsymbol{E}}^{m}-\tfrac{1}{2} ({\boldsymbol{E}}^{m}\boldsymbol{\cdot}{\boldsymbol{E}}^{m})\boldsymbol{n} \\ -S_i(({\boldsymbol{E}}^{{d},i}\boldsymbol{\cdot} \boldsymbol{n}){\boldsymbol{E}}^{{d},i}-\tfrac{1}{2} ({\boldsymbol{E}}^{{d},i}\boldsymbol{\cdot}{\boldsymbol{E}}^{{d},i})\boldsymbol{n}). \end{gathered}\right\} \end{equation}Drop velocity and centroid are computed from the volume averages
\begin{equation} {\boldsymbol{U}}_j=\frac{1}{V}\int_{V_j}{{\boldsymbol u}} \,{\rm d} V=\frac{1}{V}\int_{{\mathcal{D}}_j}\boldsymbol{n}\boldsymbol{\cdot} \left({{\boldsymbol u}}{\boldsymbol{x}}\right){\rm d} A,\quad {\boldsymbol{x}}^c_j=\frac{1}{V}\int_{V_j}{\boldsymbol{x}} \,{\rm d} V=\frac{1}{2V} \int_{{\mathcal{D}}_j}\boldsymbol{n}\left({\boldsymbol{x}}\boldsymbol{\cdot}{\boldsymbol{x}}\right){\rm d} A. \end{equation} where 
${\rm d}V$ is the volume element. To solve the system of equations (3.2) and (3.4) we use a Galerkin formulation based on a spherical harmonics representation presented in Sorgentone et al. (Reference Sorgentone, Tornberg and Vlahovska2019). All variables (position vector, velocities, electric field) are expanded in spherical harmonics which provides an accurate representation even for relatively low expansion order. In order to deal with the singular and nearly singular integrals that appear in the formulation we evoke specialized quadrature methods able to control the quadrature errors (af Klinteberg, Sorgentone & Tornberg Reference af Klinteberg, Sorgentone and Tornberg2022), and a reparametrization procedure able to ensure a high-quality representation of the drops also under deformation is used to ensure the spectral accuracy of the method (Sorgentone & Tornberg Reference Sorgentone and Tornberg2018).
3.2. Asymptotic theory for large separations
Here, we modify the asymptotic theory developed in Sorgentone et al. (Reference Sorgentone, Kach, Khair, Walker and Vlahovska2021); Kach et al. (Reference Kach, Walker and Khair2022) to dissimilar drops. We first evaluate the electrostatic interaction of two widely separated spherical drops. In this case, the drops can be approximated by point dipoles. The disturbance field 
${\boldsymbol {E}}_1$ of the drop dipole 
${\boldsymbol {P}}_1$ induces a dielectrophoretic (DEP) force on the dipole 
${\boldsymbol {P}}_2$ located at 
${\boldsymbol {x}}^c_2=d\hat {\boldsymbol {d}}$, given by 
${\boldsymbol {F}}_2(d)=({\boldsymbol {P}}_2\boldsymbol {\cdot} \boldsymbol {\nabla } {\boldsymbol {E}}_1)|_{r=d}$. Likewise, dipole 
${\boldsymbol {P}}_2$ induces a force on dipole 1 that is of equal magnitude and opposite sign 
${\boldsymbol {F}}_1=-{\boldsymbol {F}}_2$ (i.e. the DEP interaction is reciprocal). The drop velocity under the action of this force can be estimated from Stokes law, 
${\boldsymbol {U}}_i={\boldsymbol {F}}_i /\zeta _i$ where 
$\zeta$ is the friction coefficient 
$\zeta _i=6{\rm \pi} (3\lambda _i+2)/(3(\lambda _i+1))$. Thus,
\begin{equation} \left.\begin{gathered} {\boldsymbol{U}}_i^{dep}= \frac{2\beta_D}{d^4}\left(\frac{3(1+\lambda_i)}{2+3\lambda_i}\right) [(1-3\cos^2\varTheta)\hat{\boldsymbol{d}}-\sin(2\varTheta)\hat{\boldsymbol{t}}],\\ \beta_D= \left(\frac{R_1-1}{R_1+2}\right)\left(\frac{R_2-1}{R_2+2}\right). \end{gathered}\right\} \end{equation}
If 
$(\textit {R}_1-1)(\textit {R}_2-1)>0$, as in the case of identical droplets, droplets attract if 
$\varTheta <\varTheta _c=\arccos (\frac {1}{\sqrt {3}})\approx 54.7^{\circ }$, e.g. when the drops are lined up with the field, and repel if the line-of-centres of the two drops is perpendicular to the applied field. The droplets line-of-centres rotates to align with the applied field. However, this situation reverses if 
$(\textit {R}_1-1)(\textit {R}_2-1)<0$: the droplets repel if their line-of-centres is parallel to the applied field direction, and attract if their line-of-centres is perpendicular to the field. The DEP interaction in this case rotates the droplet line-of-centres away from the applied field direction.
The electrohydrodynamic (EHD) flow due to droplet 1 moves droplet 2 and, vice versa, the flow due to droplet 2 moves droplet 1. The velocities of the droplets are
\begin{equation} {{\boldsymbol{U}}}^{ehd}_{2} = \beta_{T,1} {\boldsymbol{U}}^{ehd}\left(d,\lambda_2\right),\quad {{\boldsymbol{U}}}^{ehd}_{1} ={-}\beta_{T,2} {\boldsymbol{U}}^{ehd}\left(d,\lambda_1\right) \end{equation}where at leading order in separation
\begin{equation} {\boldsymbol{U}}^{ehd}\left(d,\lambda\right)= \left(\frac{1}{d^2}-\frac{2}{d^4} \left(\frac{1+3\lambda}{2+3\lambda}\right)\right)({-}1+3 \cos^2\varTheta)\hat{\boldsymbol{d}}-\frac{2}{d^4} \left(\frac{1+3\lambda}{2+3\lambda}\right)\sin(2\varTheta)\hat{\boldsymbol{t}}, \end{equation}and the stresslet magnitude is
\begin{equation} \beta_{T,i}=\frac{9}{10}\frac{R_i-S_i}{\left(1+\lambda_i\right)\left(R_i+2\right)^2},\quad i=1,2. \end{equation}
For equiviscous droplets, the relative velocity 
${\boldsymbol {U}}_2-{\boldsymbol {U}}_1$ shows that the EHD interaction changes sign (attractive to repulsive or vice versa depending on 
$\beta _{T,2}+\beta _{T,1}$) at the same critical angle 
$\varTheta _c$ as the DEP case. However, the EHD interaction also changes sign at separation 
$d^2_c=2({1+3\lambda })/({2+3\lambda })$. Here 
$d_c$ ranges from 1 for bubbles (
$\lambda =0$) to 
$\sqrt {2}$ for very viscous drops (
$\lambda \rightarrow \infty$), both corresponding to centre-to-centre distance smaller than the minimal separation of 2 for spherical drops. Accordingly, in reality the sign of the EHD interactions does not vary with drop-drop separation. For droplets aligned with the field, both 
$\beta _T$ negative results in EHD attraction, since the surface flow around each drop is pole-to-equator and the fluid is being drawn away from the space between the droplets. Both 
$\beta _T$ positive results in repulsion because the surface flow around the droplets is equator-to-pole and the fluid is being drawn into the space between the droplets, effectively pushing them away. Dissimilar droplets can either attract or repel depending on the relative strength of their stresslets. These scenarios reverse for droplet with line-of-centres perpendicular to the applied field direction.
4. Results
Here, we apply the asymptotic theory to dissimilar drops and show that the EHD interactions turn the pair into a self-propelled entity. Using the analytical theory and numerical simulations, we find that, depending on the initial configuration, droplets can form a translating pair with steady separation between the droplets, can ‘chase’ each other or come into contact.
4.1. Droplet cooperative propulsion
An isolated, charge-neutral drop in a uniform electric field does not move. The proximity of a boundary (Yariv Reference Yariv2006) or another drop breaks the symmetry and can cause droplet motion. However, if the drops are identical there is no net migration, i.e. their centre of mass remains stationary. Difference in droplet size and properties gives rise to cooperative droplet propulsion and motion of the pair's centre of mass. The ‘swimming’ velocity of the pair at leading order in separation and same-size droplets is
\begin{equation} {\boldsymbol{U}}_{{cm}}=\tfrac{1}{2}\left({\boldsymbol{U}}_2+{\boldsymbol{U}}_1\right)=f(d)({-}1+3 \cos^2\varTheta)\hat{\boldsymbol{d}}+g(d) \sin(2\varTheta)\hat{\boldsymbol{t}}. \end{equation}A simple difference in droplets’ viscosities, if all other properties and drop radii are the same, gives rise to centre-of-mass motion with
\begin{equation} f(d)=g(d)=\frac{3(\lambda_2-\lambda_1)}{d^4(2+3\lambda_1)(2+3\lambda_2)}\left(\beta_T+\beta_D\right). \end{equation}
The pair's migration direction and speed are controlled by the relative importance of the induced dipole and the EHD stresslet. The EHD flow weakens with increasing conductivity, and for 
$\textit {R}\rightarrow \infty$, 
$\beta _T+\beta _D\rightarrow 1$. The DEP interaction vanishes at 
$\textit {R}=1$, and in this case the pair translation is solely driven by the interaction of the droplets’ stresslet flows.
Here, we focus on droplets with same size and viscosity but different conductivities and permittivities. In this case, the DEP contributions to the centre-of-mass velocity cancel and the pair migration speed is set by the droplet stresslets:
\begin{equation} {\boldsymbol{U}}_{cm}=\tfrac{1}{2}\left({\boldsymbol{U}}_2+{\boldsymbol{U}}_1\right)=\tfrac{1}{2}(\beta_{T,1}-\beta_{T,2}){\boldsymbol{U}}^{ehd}. \end{equation}
Hence, the non-reciprocal EHD interaction is the source of the droplet tandem locomotion; the pair migration speed vanishes if the droplet stresslets are the same. The direction of motion is determined by the stresslets’ difference. For example, droplets with 
$\textit {R}_1=0.1$, 
$\textit {R}_2=100$ and same permittivity ratio (
$\textit {S}_1=\textit {S}_2=1$) that are initially aligned with the field translate antiparallel to the field (see supplementary movie 1 available at https://doi.org/10.1017/jfm.2022.875); swapping the droplets reverses the direction of propulsion. In this case, droplets settle into a stable separation. In general, however, the drop pair dynamics is complex because the centre-of-mass motion is superimposed on changes in separation and rotation of the line-of-centres relative to applied field direction.
4.2. Droplet trajectories
Here we examine the conditions to form a stable locomoting tandem. According to the theory, the droplet separation and line-of-centre orientation evolve as
\begin{gather} \dot d={\boldsymbol{U}}\boldsymbol{\cdot}\hat{\boldsymbol{d}}=\left[(\beta_{T,1}+\beta_{T,2}) \frac{1}{d^2}-\frac{2}{d^4}\varPhi(\lambda, \beta_{T,1}, \beta_{T,2}, \beta_{D})\right]({-}1+3\cos^2\varTheta), \end{gather}
\begin{gather}\dot \varTheta=\frac{1}{d}{\boldsymbol{U}}\boldsymbol{\cdot}\hat{\boldsymbol{t}}={-}\frac{2}{d^5} \varPhi(\lambda, \beta_{T,1}, \beta_{T,2}, \beta_{D})\sin\left(2\varTheta\right), \end{gather}
where 
${\boldsymbol {U}}={\boldsymbol {U}}_2-{\boldsymbol {U}}_1$ is the relative velocity and
\begin{equation} \varPhi(\lambda, \beta_{T,1}, \beta_{T,2}, \beta_{D})= \left(\frac{1+3\lambda}{2+3\lambda}\right)\left[(\beta_{T,1}+\beta_{T,2})+2\beta_D \left(\frac{3(1+\lambda)}{1+3\lambda}\right)\right]. \end{equation}
Examination of this dynamical system shows that there are two equilibrium points: 
$\varTheta _*=0$ and 
$d_*=d_{eq}$, and 
$\varTheta _*={\rm \pi} /2$ and 
$d_*=d_{eq}$, where (for viscosity ratio 1)
\begin{equation} d^2_{eq}=\frac{8}{5}\!+\!\frac{32 (R_1-1) (R_1\!+\!2) (R_2-1) (R_2\!+\!2)}{3 (R_1^2 (R_2-S_2)\!+\!R_1 (R_2 (R_2+8)-4 S_2\!+\!4)+R_2 (4\!-\!(R_2\!+\!4) S_1)\!-\!4 (S_1\!+\!S_2))}. \end{equation}
The equilibrium points are saddles, as seen from the phase plane plotted in figure 2(a). If the droplet line-of-centres is initially aligned with the applied field direction, the droplets attain a steady separation for values of the droplet conductivities corresponding to figure 2(c) and the left branch of figure 2(b). If the droplet line-of-centres is initially perpendicular to the applied field direction, the steady separation is given by the right branch of figure 2(b). In these scenarios, the DEP is repulsive and stronger than the EHD at short separations and the EHD is attractive at large separations. Accordingly, the drops attract or repel until they reach the equilibrium separation 
$d_{eq}$.
Figure 2. (a) Phase plane of droplet trajectories for 
$\textit {R}_1=0.1$, 
$\textit {R}_2=100$ and 
$\textit {S}_1=\textit {S}_2=1$, corresponding to 
$\varPhi <0$. The dashed lines correspond to 
$d_{eq}$ and 
$\varTheta _c$. (b,c) Equilibrium separation 
$d_{eq}$ for drops with different conductivity but same permittivity 
$\textit {S}_1=\textit {S}_2=1$.
Any misalignment of the drops line-of-centres drives the droplets away from the equilibrium configurations towards contact or infinite separation. The trajectories 
$d(\varTheta )$ are given by
\begin{equation} d^2(\varTheta, d_0, \varTheta_0)=\frac{f(\varTheta, d_0, \varTheta_0)}{1+2 b f(\varTheta, d_0, \varTheta_0)}, \end{equation}where
\begin{equation} f(\varTheta, d_0, \varTheta_0)=\frac{d_0^2}{1-2 b d_0^2} \left(\frac{\cos\varTheta\sin^2 \varTheta}{\cos\varTheta_0\sin^2 \varTheta_0}\right),\quad b=\frac{5(\beta_{T,1}+\beta_{T,2})}{16(\beta_{T,1}+\beta_{T,2}+3\beta_D)}. \end{equation}
The trajectories for different initial configurations – angle and separation – are plotted in figure 2(a). In the case 
$\varPhi <0$, if the initial separation 
$d_0>d_{eq}$, the droplets either initially attract but then separate indefinitely if 
$\varTheta _0<\varTheta _c$ or monotonically separate if 
$\varTheta _0>\varTheta _c$. This scenario is reversed if 
$d_0< d_{eq}$, where droplets ultimately come into contact. However, if the drops are misaligned but separated exactly by 
$d_{eq}$, i.e. 
$d=d_0$, their separation remains constant while their line-of-centres rotates continuously towards the equilibrium points, either 
$\varTheta _*=0$ , if 
$\varPhi >0$, or 
$\varTheta _*={\rm \pi} /2$, in the opposite case, see (4.5).
Figures 3 and 5 illustrate the droplet dynamics for the cases 
$\varPhi <0$ and 
$\varPhi >0$, respectively. While for 
$\varTheta _0=0$ and 
$\varPhi <0$ the drops’ line-of-centres remains aligned with the applied field and the drops form a pair with steady separation, if 
$\varTheta _0\ne 0$, the misaligned droplets migrate towards a configuration where the line-of-centres is nearly perpendicular to the field, see figures 3 and 4(a,e). For this system, the induced dipoles are in opposite directions, i.e. 
$\beta _D<0$. The stresslets also have opposite sign (
$\beta _{T,1}<0$ and 
$\beta _{T,2}>0$); however, the EHD flow for the 
$\textit {R}=0.1$ droplet is stronger (
$|\beta _{T,1}|>\beta _{T,2}$). If 
$\varTheta _0<\varTheta _c$, the DEP is repulsive, while the pole-to-equator surface EHD flow of drop 1 results in attraction between the drops. If the initial distance between the drops is greater than the equilibrium separation, 
$d_0>d_{eq}$, the interaction is initially dominated by the EHD and droplets attract. However, as they get closer the DEP repulsion intensifies. Moreover, while the relative distance between the drops is changing, the line-of-centres is also rotating and the angle with the applied field direction enters the 
$\varTheta >\varTheta _c$ range, where the EHD is repulsive. As a result, the droplets indefinitely separate, while drop 1 is ‘chasing’ drop 2, with decreasing propulsion speed. The trajectory in this case is shown in figure 3(a). The non-monotonic evolution of the separation, relative radial and propulsion velocities is illustrated in figure 4(a–d). The long-time behaviour is dominated by the stresslets’ flow, which results in separation increasing as 
$\sim t^{1/3}$, see (4.4) and inset in figure 4(b). Accordingly, the centre-of-mass and relative radial speeds approach zero as the power law 
$\sim t^{-2/3}$. For example, the radial relative velocity decreases as 
${\boldsymbol {U}}\boldsymbol {\cdot}\hat {\boldsymbol {d}}= (\beta _{T,1}+\beta _{T,2})^{1/3}(3t)^{-2/3}$ for 
$t\gg 1$. If the initial distance between the drops is less than the equilibrium separation, 
$d_0< d_{eq}$, the interactions are reversed: the droplets initially repel and then attract (if 
$\varTheta _0<\varTheta _c$), or monotonically attract (if 
$\varTheta _0>\varTheta _c$), and eventually come in contact. The trajectory in this case is shown in figure 3(b). The droplets’ relative velocity increases rapidly as they approach each other, see figure 4(g). The centre-of-mass speed varies along the trajectories and it is minimal when the radial velocity is close to zero, see figure 4(d,h). Comparison of the numerical and theoretical results shows that the asymptotic theory qualitatively captures the drop dynamics. The agreement between simulations and theory is better for droplets that are initially farther apart. Thus, given the high computational costs of the simulations, the theory can be used to estimate droplet interactions. Droplet deformation increases the 
$d_{eq}$, above which drops evolve towards separating state. In the considered example, we found by numerical simulations that 
$d_0=6$ also leads to contact since the deformation causes the drops to get too close and unable to escape the DEP attraction which ultimately leads to contact.
Figure 3. Droplet trajectories for 
$\textit {R}_1=0.1$ (red), 
$\textit {R}_2=100$ (blue), 
$\textit {S}_1=\textit {S}_2=1$, 
$\lambda _1=\lambda _2=1$ and initial angle 
$\varTheta _0=10^{\circ }$. For this system 
$d_{eq}=4.94$. Lines are computed from the asymptotic theory. Symbols correspond to the numerical simulations. Initial separation is 
$d_0=7$ (a), see movie 2, and 
$d_0=4$ (b), see movie 3. Simulation time is 4161 (a) and 967 (b). Theoretical calculation ends at time 10 000 (a) and 2827 (b).
Figure 4. Droplet-pair dynamics for 
$\textit {R}_1=0.1$ (a–d), 
$\textit {R}_2=100$ (e–h), 
$\textit {S}_1=\textit {S}_2=1$, 
$\lambda _1=\lambda _2=1$ and initial angle 
$\varTheta _0=10^{\circ }<\varTheta _c$. For this system 
$d_{eq}=4.94$. Initial separation is 
$d_0=7$ (a–d) and 
$d_0=4$ (e–h). Angle (a,e), separation (b,f), relative radial velocity (c,g) and absolute value of the pair migration speed (d,h) as a function of time scaled by the EHD scale 
$\eta_{m}/({\varepsilon }_{m} E_0^2)$. Lines are computed from the asymptotic theory. Symbols correspond to the numerical simulations. Insets in (b,d) show the long-time asymptotic behaviour (red).
Figure 5. Droplet-pair dynamics for 
$\textit {R}_1=2$ (red), 
$\textit {R}_2=100$ (blue), 
$\textit {S}_1=\textit {S}_2=1$, 
$\lambda _1=\lambda _2=1$, and initial angle 
$\varTheta _0=80^{\circ }>\varTheta _c$. Initial separation (a) 
$d_0=4$, leading to contact, and (b) 
$d_0=7$, leading to indefinite separation. Symbols correspond to the numerical simulations (end time 2200). Trajectories computed from the asymptotic theory are over period of time 2227 (a) and 500 000 (b).
In the case of droplets pairing in the transverse direction, 
$\varPhi >0$, the droplets exhibit the opposite orientational behaviour to the 
$\varPhi <0$ case, and move to align with the field. If the initial separation is smaller than the equilibrium one, drops come in contact, and otherwise separate indefinitely, while both drops move in opposite direction (‘run away’ from each other), see figure 5. In the latter case, the interaction is extremely weak. The trajectory time is 500 000, which is prohibitively expensive to simulate numerically.
5. Conclusions
We analyse the interactions of dissimilar droplets in a uniform electric field by means of an asymptotic theory, assuming spherical droplets (
$\textit {Ca}\ll 1$) and large separations, and numerical simulations, using a three-dimensional boundary integral method. The simulations for 
$\textit {Ca}=0.1$ qualitatively agree with the theory, and thus the theory can be used for a fast estimate of the drop trajectories. Our study focuses on the effect of the mismatch in the electric properties, considering drops with different conductivity and permittivity but same size and viscosity. In this case, the non-reciprocal EHD interactions give rise to a net motion of the drop pair. The centre-of-mass motion is accompanied by changes in drop separation and angle between their line-of-centres to the applied field direction, which gives rise to intricate trajectories. Depending on the droplet stresslets and dipoles in the function 
$\varPhi$, defined by (4.6), drops tend to orient their line-of-centres either parallel, if 
$\varPhi >0$, or perpendicular to the applied field direction, if 
$\varPhi <0$. Initial separation determines if drops will coalesce or indefinitely separate. For drops with 
$\varPhi <0$, if 
$d_0>d_{eq}$ and 
$\varTheta _0<\varTheta _c$, drops initially attract and then separate indefinitely, while chasing each other. If 
$d_0< d_{eq}$ and 
$\varTheta _0<\varTheta _c$, droplets repel and then attract until contact; the interaction is purely attractive and separation decreases monotonically if 
$\varTheta _0>\varTheta _c$. In the particular case of drops aligned with the field and 
$\varPhi <0$, the drops reach steady separation and ‘swim’ along the applied field direction, if 
$\textit {R}_1<\textit {R}_2$; direction of motion is reversed if 
$\textit {R}_1>\textit {R}_2$. If 
$\varPhi >0$ and the drops’ line-of-centres is perpendicular to the applied field direction, the droplets form a tandem translating transversely to the field. If instead the drops line-of-centres is initially misaligned with the applied field direction, with 
$d_0>d_{eq}$ and 
$\varTheta _0>\varTheta _c$, drops initially attract, then repel indefinitely while moving in opposite directions to each other. If 
$d_0< d_{eq}$ and 
$\varTheta _0>\varTheta _c$, droplets first repel and then attract until contact. In both cases, the separation changes monotonically if 
$\varTheta _0<\varTheta _c$.
Our work represents the first numerical study of the three-dimensional dynamics of electrically dissimilar drops and opens new directions of exploration of how to manipulate droplets and direct assembly of particles with electric fields. For example, the preferential migration to or away from the applied field direction could be exploited to separate droplets with different conductivities and/or permittivities. The self-propelled droplet pair, which is powered but not directed by the applied electric field, can inspire new designs of active units and, more generally, non-equilibrium structures driven by non-reciprocal interactions (Ivlev et al. Reference Ivlev, Bartnick, Heinen, Du, Nosenko and Löwen2015; Meredith et al. Reference Meredith, Moerman, Groenewold, Chiu, Kegel, van Blaaderen and Zarzar2020; Bowick et al. Reference Bowick, Fakhri, Marchetti and Ramaswamy2022).
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2022.875.
Funding
P.V. has been supported in part by NSF award CBET-2126498.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Migration of a pair of droplets with different size
The most common source of droplet dissimilarity is a difference in size. The centre-of-mass velocity in this case is (note that since the droplets are assumed to be neutrally buoyant, their densities are the same)
\begin{align} {\boldsymbol{U}}_{cm} &= \frac{a_1^3{\boldsymbol{U}}_1+a_2^3{\boldsymbol{U}}_2}{a_1^3+a_2^3} \nonumber\\ &=(\beta_{T} f^{ehd}(d)+\beta_D f^{{dep}}(d)) ({-}1+3 \cos^2\varTheta)\hat{\boldsymbol{d}} \nonumber\\ &\quad +(\beta_{T} g^{ehd}(d)+\beta_D g^{{dep}}(d))\sin(2\varTheta)\hat{\boldsymbol{t}}. \end{align}Let us consider the simplest scenario of droplets with same material properties (viscosity, conductivity and permittivity). In this case, even though the DEP forces are reciprocal, the droplet velocities are different due to different drag. Accordingly,
\begin{equation} f^{{dep}}(d)={-}g^{{dep}}(d)= \frac{6\nu^3(1+\lambda)(\nu^2-1)}{d^4(2+3\lambda)(1+\nu^3)}, \end{equation}
where 
$\nu =a_2/a_1$ is the ratio of droplet radii. The EHD contribution is
\begin{equation} \left.\begin{gathered} f^{{{ehd}}}(d)= \frac{(\nu^3-1)}{d^2(1+\nu^3)}+g^{{ehd}}(d),\\ g^{{ehd}}(d)={-}\frac{1}{d^4(1+\nu^3)}\left((\nu^3-1)+\frac{3\lambda(\nu^7-1)}{2+3\lambda}\right). \end{gathered}\right\} \end{equation}
Open access
