# Critical charges for droplet collisions

A. Dubey1, G. P. Bewley2, K. Gustavsson1, B. Mehlig1 1Department of Physics,
Gothenburg University, 41296 Gothenburg, Sweden
2Sibley School of Mechanical and Aerospace Engineering, Cornell University,
USA

###### Abstract

The collision efficiency of uncharged micron-sized water droplets in air is
determined by the breakdown of hydrodynamics at droplet separations of the
order of the mean-free path, by van-der-Waals forces, or a combination of the
two. In contrast, electrostatic forces determine the collision efficiency of
charged droplets if the charge is large enough. To find the charge for which
the transition to charge-dominated collisions occurs, we computed the
collision efficiency of charged, hydrodynamically-interacting droplets
settling in quiescent air, including the breakdown of hydrodynamics at small
interfacial distances. For oppositely charged droplets, the transition occurs
when a saddle point of the relative droplet-dynamics exits the region where
the hydrodynamics breaks down. For droplets with radii $16\,\mu$m and
$20\,\mu$m, this occurs at $\sim 10^{3}$ elementary charges $e$. For smaller
charges, the collision efficiency depends upon the Kn number (defined as the
ratio of the mean-free-path of air to the mean droplet radius), whereas for
larger charges it does not. For droplets charged with the same polarity, the
critical charge is $\sim 10^{4}\,e$ for the above radii.

##  I Introduction

How do micron-sized cloud droplets grow? This question relates to the
fundamental mechanisms that determine droplet-size distributions in
atmospheric clouds Pruppacher and Klett (1997); Shaw (2003); Devenish  _et
al._ (2012). In clouds with droplets of different sizes, droplet collisions
occur due to differential settling. This mechanism can cause rapid droplet
growth Berry and Reinhardt (1974); Grabowski and Wang (2013). An open question
is how this process is initiated, how droplet-size differences develop in the
first place. Saffman & Turner explained that micron-sized droplets of similar
sizes can collide if turbulent strains bring them together Saffman and Turner
(1956), but this process is very slow on average. Since the sequence of
droplet collisions in a cloud is random, the collision times are essentially
Poisson distributed, and fluctuations may nevertheless result in some droplets
that grow very rapidly Kostinski and Shaw (2005).

An uncertainty in these considerations is the collision efficiency, the ratio
of the collision cross section of hydrodynamically interacting droplets to the
geometric collision cross section $\propto(a_{1}+a_{2})^{2}$, for two droplets
with radii $a_{1}$ and $a_{2}$. The point is that hydrodynamic interactions
cause the paths of approaching droplets to curve around each other Kim and
Karrila (1991), reducing the cross section. This hydrodynamic repulsion
diverges at very small distances, apparently preventing the approaching
droplets from colliding at all. The divergence is unphysical though, because
the hydrodynamic approximation breaks down at droplet separations of the order
of the mean free path of air Sundararajakumar and Koch (1996). As a result of
this non-continuum effect, the collision efficiency depends on the Knudsen
number Kn, equal to the mean free path divided by the droplet radius.

Hydrodynamic interactions and the breakdown of hydrodynamics result in a
significant reduction of the Saffman-Turner collision rate Dhanasekaran  _et
al._ (2021). Its dependence on Kn and on the non-dimensional relative settling
speed can be understood in terms of an intricate sequence of bifurcations of
the collision dynamics Dubey  _et al._ (2022). Many other factors are thought
to play a role for the collision dynamics of settling droplets. Refs.
Dhanasekaran  _et al._ (2021); Dubey  _et al._ (2022) neglected particle
inertia. This works for very small droplets, but particle inertia allows
larger droplets to detach from the flow, leading to larger collision rates
Sundaram and Collins (1997); Falkovich  _et al._ (2002); Gustavsson and Mehlig
(2016). Klett & Davis Klett and Davis (1973) showed that the collision rate of
similar-sized droplets depends sensitively on the particle Reynolds number, a
measure of the effect of convective fluid inertia upon the disturbance flow
caused by the settling droplets. The reason is that the Stokes problem for the
collision dynamics is degenerate: equal-sized droplets neither approach nor
diverge. Convective fluid inertia breaks this degeneracy Candelier (2007).
Van-der-Waals forces are short-range attractive forces due to polarisation of
the water molecules comprising the droplets Rother  _et al._ (2022) that can
increase the collision rate at small Knudsen numbers.

Small water droplets tend to carry excess electrical charges, sometimes
positive, sometimes negative. It is clear that electrostatic forces either
reduce or increase the collision efficiency, depending on whether the two
droplets have excess charges of the same polarity or not. But how much charge
is required to make a significant difference to the collision dynamics? This
question is motivated in part by the fact that water droplets in atmospheric
clouds are charged Harrison  _et al._ (2015). The question has a long history,
but it has not been resolved. Davis Davis (1965) found for instance that more
than 800 elementary charges $e$ are required for the Coulomb force to change
the collision efficiency for radii of the order of $20\,\mu$m, the typical
size of cloud droplets, yet Tinsley and Zhou (2014) suggest that it may take
at least $10^{4}\,e$. These differences are not surprising, given that the
estimates are based on highly idealised models that either do not account for
how the hydrodynamic forces change as the droplets come close Tinsley and Zhou
(2014); Paluch (1970); Schlamp  _et al._ (1976, 1979), or do not consider how
the hydrodynamic approximation breaks down below the mean free path Davis
(1965); Tinsley and Zhou (2014); Paluch (1970); Schlamp  _et al._ (1976,
1979), causing the predictions to fail Abbott (1974).

Magnusson et al. Magnusson  _et al._ (2022) used models for hydrodynamic and
electrical interactions valid at large droplet separations to analyse
collisions of strongly charged $20\mu$m-droplets, with opposite polarities,
settling in quiescent air. In this regime, the collision efficiency is
determined by a stable manifold of a saddle point of the relative droplet
dynamics. At the saddle point, the larger droplet travels below the smaller
one – so that Coulomb attraction, hydrodynamic interactions, and the
difference of the gravity forces cancel, resulting in a steady state. When the
charges are strong enough, the droplets are far apart at this equilibrium
point. In this case, non-continuum effects do not matter, and it is sufficient
to consider the large-distance asymptotes of the hydrodynamic interactions.
These approximations break down at smaller charges. Therefore we cannot use
the model analysed by Magnusson et al. Magnusson  _et al._ (2022) to estimate
the smallest amount of charge needed to significantly change the collision
dynamics of micron-sized droplets.

In order to determine the collision efficiencies of weakly charged droplets
settling under gravity, we performed numerical simulations of the collision
dynamics, taking into account particle inertia, hydrodynamic interactions not
only at large but also at small droplet separations, as well as their
regularisation by non-continuum effects due to the breakdown of hydrodynamics
at interfacial distances smaller than the mean free path. We modeled the
droplets as perfect conductors Lekner (2012), which allowed us to consider the
effect of induced charges that affect the electrostatic force at small
separations. This is important, because induced charges can cause droplets to
attract each other even though they are charged with the same polarity Khain
_et al._ (2004); Fletcher (2013). We found that the collision efficiency
exhibits two qualitatively different regimes: a small-charge regime where non-
continuum effects determine collisions and a large-charge regime where
electrical forces determine collisions. We determined how much charge is
needed to affect this qualitative change in the collision dynamics. The
transition region is quite broad as the charge varies, but the qualitative
change occurs when a saddle point of the relative droplet-dynamics exits the
spatial region where non-continuum effects dominate. The critical charge
depends on whether the droplets are charged with the same or with opposite
polarities. For two droplets with radii $16\,\mu$m and $20\,\mu$m charged with
the opposite polarities, the transition occurs at $\sim 10^{3}$ elementary
charges $e$. When the droplets have excess charges with the same polarity,
$\sim 10^{4}\,e$ are required. In addition, a bifurcation occurs at about
$1.5\times 10^{4}$ elementary charges, leading to two saddle points above the
collision sphere. Droplets approaching near the symmetry axis – the
$R_{3}$-axis in Fig. 1(a) – spend a long time in the vicinity of one of the
saddle points, before they diverge along its unstable manifold.

##  II Model

###  II.1 Equations of motion

\begin{overpic}[scale={.4}]{Fig1.pdf} \end{overpic} Figure 1:  Panel (a) shows
a schematic of two droplets with radii $a_{1}<a_{2}$ settle in a quiescent
fluid. The separation vector between their centres-of-mass is $R$. Also shown
is the collision sphere around the smaller droplet (dashed). The droplets
collide if the centre-of-mass of the larger one hits the collision sphere.
Gravity points in the negative $R_{3}$-direction. Panel (b) illustrates two
regimes of the collision dynamics, distinguished by the location of the fixed
point S$1$ (see text). For weak charges, S$1$ is inside the region where
hydrodynamics breaks down (hashed area), at non-dimensional distance
$s^{\ast}\equiv R^{\ast}-2<{\rm Kn}$. For strong charges, the fixed point lies
outside, $s^{\ast}>{\rm Kn}$.

The motions of two spherical droplets $i=1,2$ are determined by Newton’s
equations,

$\displaystyle\dot{\mbox{\boldmath$x$}}^{(i)}$
$\displaystyle=\mbox{\boldmath$v$}^{(i)}\,,$ (1a)
$\displaystyle\dot{\mbox{\boldmath$v$}}^{(i)}$
$\displaystyle=[\mbox{\boldmath$F$}_{\rm g}^{(i)}+\mbox{\boldmath$F$}_{\rm
h}^{(i)}+\mbox{\boldmath$F$}_{\rm e}^{(i)}]/m_{i}\,,$ (1b)
$\displaystyle\dot{\mbox{\boldmath$\omega$}}^{(i)}$
$\displaystyle=\mbox{\boldmath$T$}_{\rm h}^{(i)}/I_{i}\,.$ (1c)

Here $\mbox{\boldmath$x$}^{(i)}$ and $\mbox{\boldmath$v$}^{(i)}$ are the
position and velocity of droplet $i$, and ${\mbox{\boldmath$\omega$}}^{(i)}$
is its angular velocity. Droplet $i$ has mass
$m_{i}=\tfrac{4\pi}{3}\varrho_{\rm w}a_{i}^{3}$ and moment of inertia
$I_{i}={\tfrac{2}{5}m_{i}a_{i}^{2}}$, where $a_{i}$ is the radius of droplet
$i$, and $\varrho_{\rm w}$ is the mass density of water. On the r.h.s. of Eq.
(1b), $\mbox{\boldmath$F$}_{\rm g}^{(i)}=m_{i}\mbox{\boldmath$g$}$ is the
gravity force with gravitational acceleration $g$, and
$\mbox{\boldmath$F$}_{\rm e}^{(i)}$ is the electrostatic force on droplet $i$,
$\mbox{\boldmath$F$}_{\rm h}^{(i)}$ is the hydrodynamic force, and
$\mbox{\boldmath$T$}_{\rm h}^{(i)}$ is the hydrodynamic torque. The droplets
settle through quiescent air. We denote their centre-of-mass separation by
$\mbox{\boldmath$R$}=\mbox{\boldmath$x$}^{(2)}-\mbox{\boldmath$x$}^{(1)}$
(Fig. 1), and their relative velocity by
$\mbox{\boldmath$V$}=\mbox{\boldmath$v$}^{(2)}-\mbox{\boldmath$v$}^{(1)}$.

We non-dimensionalise the problem as follows. As a length scale we use
$\overline{a}=(a_{1}+a_{2})/2$. Velocities are non-dimensionalised with the
differential settling speed of non-interacting droplets,
$V_{0}=g{|}\gamma_{2}^{-1}-\gamma_{1}^{-1}{|}$. Here
$\gamma_{i}=(9/2)(\varrho_{\rm f}/\varrho_{\rm w})\nu/a_{i}^{2}$ is the Stokes
constant of droplet $i$, with fluid-mass density $\varrho_{\rm f}$ and
kinematic viscosity $\nu$. Time is non-dimensionalised with the characteristic
time scale $\tau_{\rm c}=\overline{a}/V_{0}$. The non-dimensional parameters
are the Stokes number ${\rm St}=(\widetilde{\gamma}\tau_{\rm c})^{-1}$
(particle inertia), and the Coulomb number ${\rm Cu}=2k_{\rm
e}|q_{1}q_{2}|/[\widetilde{m}(a_{1}+a_{2})V_{0}^{2}]$ (electrical charge).
Here $\widetilde{\gamma}=(9/2)(\varrho_{\rm f}/\varrho_{\rm
w})\nu/\overline{a}^{2}$ and
$\widetilde{m}=(4/3)\pi\overline{a}^{3}\varrho_{\rm w}$ are Stokes constant
and droplet mass based on the mean droplet radius $\overline{a}$. The Coulomb
number is the ratio of the Coulomb energy upon contact and
$\widetilde{m}V_{0}^{2}/2$. Related parameters (ratio of Coulomb to kinetic
energy) were used by Hidy and Brock (1965) to quantify the effect of electric
charge on Brownian coagulation, and in Refs. Lu  _et al._ (2010a, b); Lu and
Shaw (2015); Magnusson  _et al._ (2022) to parameterise the importance of
charge on the dynamics of micron-sized water droplets in turbulent air. The
remaining three non-dimensional parameters are the Knudsen number ${\rm
Kn}=\ell/\overline{a}$, the ratio of the mean free path of air $\ell$ to the
mean droplet radius, the radius ratio $\lambda=a_{1}/a_{2}$, and the charge
ratio $\lambda_{q}=q_{1}/q_{2}$. The charge ratio can take positive or
negative values, depending on the polarities of the excess charge on the
droplets.

###  II.2 Electrostatic force

We assume that the droplets are good conductors, so that we may use Lekner’s
results for the electrostatic force between two charged conducting spheres
Lekner (2012). We outline the main steps of the derivation in Ref. Lekner
(2012) here, using dimensional units. The electrical potential energy
$\mathscr{W}$ of two conductors with charges $q_{i}$ and potentials $V_{i}$ is
given by

$\displaystyle\mathscr{W}=\frac{k_{\rm e}}{2}\,(q_{1}V_{1}+q_{2}V_{2})\,,$ (2)

where $k_{\rm e}$ is the Coulomb constant. The charges $q_{i}$ are linearly
related to the electrical potentials $V_{j}$,

$\displaystyle\begin{bmatrix}q_{1}\\\
q_{2}\end{bmatrix}=\begin{bmatrix}C_{11}&C_{12}\\\
C_{21}&C_{22}\end{bmatrix}\begin{bmatrix}V_{1}\\\ V_{2}\end{bmatrix}\,,$ (3)

with capacitance coefficients $C_{ij}$ Jackson (1975). The capacitance
coefficients depend on the radii of the spheres, $a_{1}$ and $a_{2}$, and upon
the centre-of-mass distance between them as follows:

$\displaystyle C_{11}$
$\displaystyle=a_{1}a_{2}\sinh(U)\sum_{n=0}^{\infty}[a_{1}\sinh(nU)+a_{2}\sinh((n+1)U)]^{-1}\,,$
(4a) $\displaystyle C_{22}$
$\displaystyle=a_{1}a_{2}\sinh(U)\sum_{n=0}^{\infty}[a_{2}\sinh(nU)+a_{1}\sinh((n+1)U)]^{-1}\,,$
(4b) $\displaystyle C_{12}$
$\displaystyle=C_{21}=-a_{1}a_{2}\frac{\sinh(U)}{R}\sum_{n=1}^{\infty}[\sinh(nU)]^{-1}\,,$
(4c)

with $\cosh(U)=(R^{2}-a_{1}^{2}-a_{2}^{2})/(2a_{1}a_{2})$. Eliminating the
potentials in Eq. (2) by means of (3), one obtains the electrostatic potential
energy $\mathscr{W}$, and by differentiation of $\mathscr{W}$ w.r.t. $R$ one
obtains the electrostatic force. Truncating the sums in Eqs. (4) allows us to
compute the electrical force when droplets are far apart. We denote this
asymptote by $F_{{\rm e},\rm far}$. When $R\gg a_{1}+a_{2}$, the sums in Eq.
(4) converge rapidly. When the centre-of-mass distance is close to
$a_{1}+a_{2}$, however, this is not the case. Lekner derived an approximate
expression for the force valid at small distances $s=R-a_{1}-a_{2}$, Eq. (3.4)
in Ref. Lekner (2012):

$F_{{\rm e},{\rm
near}}=-k_{e}\,\frac{a_{1}+a_{2}}{a_{1}a_{2}s}\frac{\left[\left(\psi(\tfrac{1}{1+\lambda})+\gamma\right)q_{2}-\left(\psi(\tfrac{\lambda}{1+\lambda})+\gamma\right)q_{1}\right]^{2}}{\left[\left(\psi(\tfrac{1}{1+\lambda})+\gamma\right)\left\\{2\gamma+\log[\tfrac{2a_{1}a_{2}}{(a_{1}+a_{2})s}]\right\\}+\left(\psi(\tfrac{\lambda}{1+\lambda})+\gamma\right)\left\\{-2\psi(\tfrac{1}{1+\lambda})+\log[\tfrac{2a_{1}a_{2}}{(a_{1}+a_{2})s}]\right\\}\right]^{2}}.$
(5)

valid as $s\to 0$. Here $\psi(x)$ is the digamma function, and
$\gamma=0.5772\ldots$ is Euler’s constant. In order to obtain an approximation
valid uniformly in $R$, the two asymptotes are matched at the interfacial
distance $s_{\rm m}$ where the absolute difference between the two asymptotes
reaches a minimum, using exponential smoothing of the form

$F_{{\rm e},\rm uv}=F_{{\rm e},\rm near}\,{\rm e}^{-s/{s_{\rm m}}}+F_{{\rm
e},\rm far}\,\left(1-{\rm e}^{-s/{s_{\rm m}}}\right)\,.$ (6)

###  II.3 Hydrodynamic forces

Since the droplets settle in still air, there is no undisturbed fluid velocity
or vorticity. The hydrodynamic forces $\mbox{\boldmath$F$}_{\rm h}^{(i)}$ and
torques $\mbox{\boldmath$T$}_{\rm h}^{(i)}$ are calculated in the steady
Stokes approximation as linear functions of the translational and angular
velocities Kim and Karrila (1991),

$\displaystyle\begin{bmatrix}\mbox{\boldmath$F$}_{\rm h}^{(1)}\\\
\mbox{\boldmath$F$}_{\rm h}^{(2)}\\\ \mbox{\boldmath$T$}_{\rm h}^{(1)}\\\
\mbox{\boldmath$T$}_{\rm
h}^{(2)}\end{bmatrix}=\begin{bmatrix}\mathbb{A}_{11}&\mathbb{A}_{12}&\widetilde{\mathbb{B}}_{11}&\widetilde{\mathbb{B}}_{12}\\\
\mathbb{A}_{21}&\mathbb{A}_{22}&\widetilde{\mathbb{B}}_{21}&\widetilde{\mathbb{B}}_{22}\\\
\mathbb{B}_{11}&\mathbb{B}_{12}&\mathbb{C}_{11}&\mathbb{C}_{12}\\\
\mathbb{B}_{21}&\mathbb{B}_{22}&\mathbb{C}_{21}&\mathbb{C}_{22}\end{bmatrix}\begin{bmatrix}\mbox{\boldmath$v$}^{(1)}\\\
\mbox{\boldmath$v$}^{(2)}\\\ \mbox{\boldmath$\omega$}^{(1)}\\\
\mbox{\boldmath$\omega$}^{(2)}\end{bmatrix}.$ (7)

The elements of the resistance tensors $\mathbb{A},\mathbb{B},$ and
$\mathbb{C}$ are not independent. They are related by symmetry. Some of the
symmetry relations follow from the reciprocal theorem, while others are due to
the spherical shape of the droplets Jeffrey and Onishi (1984). A consequence
of the symmetries is that the elements. of each resistance tensor can be
written in terms of at most two scalar functions: the radial resistance
function $X_{{\alpha\beta}}$ and the tangential resistance function
$Y_{{\alpha\beta}}$. In the notation of Refs. Jeffrey and Onishi (1984); Kim
and Karrila (1991) the resistance tensors read

$\displaystyle A^{\alpha\beta}_{ij}$
$\displaystyle=-6\pi\mu\frac{a_{\alpha}+a_{\beta}}{2}[X^{A}_{\alpha\beta}\,\widehat{R}_{i}\widehat{R}_{j}+Y^{A}_{\alpha\beta}(\delta_{ij}-\widehat{R}_{i}\widehat{R}_{j})],$
(8a) $\displaystyle B^{\alpha\beta}_{ij}$
$\displaystyle=-4\pi\mu\left(\frac{a_{\alpha}+a_{\beta}}{2}\right)^{2}[Y^{B}_{\alpha\beta}\epsilon_{ijk}\widehat{R}_{k}],$
(8b) $\displaystyle\widetilde{B}^{\alpha\beta}_{ij}$
$\displaystyle=B^{\beta\alpha}_{ji},$ (8c) $\displaystyle
C^{\alpha\beta}_{ij}$
$\displaystyle=-8\pi\mu\left(\frac{a_{\alpha}+a_{\beta}}{2}\right)^{3}[X^{C}_{\alpha\beta}\,\widehat{R}_{i}\widehat{R}_{j}+Y^{C}_{\alpha\beta}(\delta_{ij}-\widehat{R}_{i}\widehat{R}_{j})].$
(8d)

The non-dimensional resistance functions $X_{{\alpha\beta}}$ and
$Y_{{\alpha\beta}}$ depend upon the non-dimensional centre-of-mass distance
$R$ between the droplets, non-dimensionalised by the mean droplet radius
$\overline{a}=(a_{1}+a_{2})/2$, and upon the radius ratio
$\lambda=a_{1}/a_{2}$. Moreover, $\widehat{R}_{j}\equiv R_{j}/R$. In total
there are 20 resistance functions in Eqs. (8):
$X^{A}_{\alpha\beta},Y^{A}_{\alpha\beta},Y^{B}_{\alpha\beta},X^{C}_{\alpha\beta}$,
and $Y^{C}_{\alpha\beta}$, all with indices $\alpha,\beta=1,2$. Ten out of
these are given in the main text of Ref. Jeffrey and Onishi (1984). The
remaining ten functions are obtained using the relations Jeffrey and Onishi
(1984)

$\displaystyle X^{A}_{\alpha\beta}(R,\lambda)$
$\displaystyle=X^{A}_{\beta\alpha}(R,\lambda)=X^{A}_{(3-\alpha)(3-\beta)}(R,1/\lambda)\,,$
(9a) $\displaystyle Y^{A}_{\alpha\beta}(R,\lambda)$
$\displaystyle=Y^{A}_{\beta\alpha}(R,\lambda)=Y^{A}_{(3-\alpha)(3-\beta)}(R,1/\lambda)\,,$
(9b) $\displaystyle Y^{B}_{\alpha\beta}(R,\lambda)$
$\displaystyle=-Y^{B}_{(3-\alpha)(3-\beta)}(R,1/\lambda)\,,$ (9c)
$\displaystyle X^{C}_{\alpha\beta}(R,\lambda)$
$\displaystyle=X^{C}_{\beta\alpha}(R,\lambda)=X^{C}_{(3-\alpha)(3-\beta)}(R,1/\lambda)\,,$
(9d) $\displaystyle Y^{C}_{\alpha\beta}(R,\lambda)$
$\displaystyle=Y^{C}_{\beta\alpha}(R,\lambda)=Y^{C}_{(3-\alpha)(3-\beta)}(R,1/\lambda)\,.$
(9e)

Jeffrey and Onishi (1984) provided series expansions in $1/R$ for the
resistance functions using the twin-multipole method, as well as small-$s$
asymptotes. Following Ref. Dhanasekaran  _et al._ (2021), the two asymptotes
are matched together using exponential smoothing as described in Section II.2
for the electrostatic force.

Townsend Townsend (2018) corrected a number of errors in the small-$s$
expressions derived by Jeffrey & Onishi Jeffrey and Onishi (1984). We used
Townsend’s corrections for the small-$s$ asymptotics of the functions
$X_{12}^{A},Y_{12}^{B},Y_{11}^{C}$ and $Y_{12}^{C}$. In particular, the
corrections mentioned by Townsend (2018) were used for the functions
$A_{12}^{X}$, $B_{12}^{Y}$, $C_{11}^{Y}$, $C_{12}^{Y}$ and the function
$g_{5}$ (see Townsend (2018) for definitions of these functions, Eqs. (15),
(23), (30), (31), and (32) in that paper. In order to verify that these
corrections provide consistent approximations for the hydrodynamic forces and
torques, we made the following checks. First, for all ten functions, we summed
the large-$R$ series expansions derived by Jeffrey and Onishi (1984),
including $150$, $200$, $250$,and $300$ terms, and confirmed that the results
were consistent with the small-$s$ asymptotes. Second, we compared with the
results for the radial and tangential mobility functions of Ref. Dhanasekaran
_et al._ (2021). We inverted Eq. (7), and set hydrodynamic forces and torques
to zero, to obtain the mobility functions for ${\rm St}\\!=\\!0$-droplets. We
obtained the functions $L$ and $M$ shown in Figs. 5 and 9 of Ref. Dhanasekaran
_et al._ (2021). This constrains the resistance functions
$X^{A}_{\alpha\beta}$, and a linear combination of $Y^{A}_{\alpha\beta}$ and
$Y^{B}_{\alpha\beta}$.

When the interfacial separations between two droplets are comparable to the
mean free path $\ell$ of air, the hydrodynamic approximation breaks down.
Instead, the droplets move in a non-continuum flow described by the Boltzmann
equation. The resulting non-continuum corrections to the radial hydrodynamic
forces between two spheres were first computed by Sundararajakumar and Koch
(1996). They used solutions to the linearised Boltzmann equation to calculate
the non-continuum radial resistance functions for interfacial separations
smaller than the mean free path. The corresponding tangential corrections were
evaluated by Li Sing How  _et al._ (2021). In our numerical computations we
used the uniformly valid resistance functions (denoted by the superscript
${\rm uv}$) quoted in Ref. Li Sing How (2021), their Eqs. (4.12) to (4.15):

$\displaystyle X_{\alpha\beta}^{A,{\rm uv}}$
$\displaystyle=X_{\alpha\beta}^{A}+(-1)^{\alpha+\beta}a_{\alpha\beta}\Big{(}2\frac{f^{||}_{\rm
fit}}{\rm Kn}-\frac{1}{\xi}\Big{)}\,,$ (10a) $\displaystyle
Y_{\alpha\beta}^{A,{\rm uv}}$
$\displaystyle=Y_{\alpha\beta}^{A}+(-1)^{\alpha+\beta}a_{\alpha\beta}\Big{[}\frac{1}{3\sqrt{\pi}}\widetilde{W}+\Big{(}\frac{1-\lambda}{1+\lambda}\Big{)}^{2}\widetilde{Q}\Big{]}\,,$
(10b) $\displaystyle Y_{\alpha\beta}^{B,{\rm uv}}$
$\displaystyle=Y_{\alpha\beta}^{B}+2a_{\alpha\beta}\Big{(}\frac{a_{\alpha}}{\overline{a}}\Big{)}^{|\alpha-\beta|}\Big{[}\frac{(-1)^{\beta}}{4\sqrt{\pi}}\widetilde{W}-(-1)^{(\alpha+\beta)}\frac{3}{4}\Big{(}\frac{1-\lambda}{1+\lambda}\Big{)}^{2}\widetilde{Q}\Big{]}\,,$
(10c) $\displaystyle Y_{\alpha\beta}^{C,{\rm uv}}$
$\displaystyle=Y_{\alpha\beta}^{C}+a_{\alpha\beta}\Big{(}\frac{2a_{H}}{\overline{a}}\Big{)}^{|\alpha-\beta|}\Big{[}\frac{1}{4\sqrt{\pi}}\widetilde{W}+(-1)^{(\alpha+\beta)}\frac{3}{4}\widetilde{Q}\Big{]}\,.$
(10d)

The above equations describe corrections to 16 out of the 20 resistance
functions mentioned above.The radial rotational resistance function
$X^{C}_{\alpha\beta}$ is not modified Li Sing How (2021). The functions
$f^{||}_{\rm fit},\widetilde{W},\widetilde{Q}$ are given in Ref. Li Sing How
(2021). In the above equations, $\overline{a}$ is the mean droplet radius
defined above, $a_{H}=a_{1}a_{2}/(a_{1}+a_{2})$ is the harmonic mean, and
$a_{\alpha\beta}=2a_{H}/(a_{\alpha}+a_{\beta})$ Li Sing How (2021). The non-
continuum resistance functions mentioned above are computed using the
linearised Boltzmann equation, valid when the Mach number remains small. An
additional simplification used is that the collision term in the Boltzmann
equation is replaced by a collision operator which relaxes the molecular
velocity to a Maxwellian distribution Bhatnagar  _et al._ (1954). Finally, the
results are valid for asymptotically small Kn numbers.

The resistance functions coupling translations and rotations in Eq. (7) are
subleading in $1/R$ to the translational resistance functions Kim and Karrila
(1991). When the droplets are far apart, the angular velocities may be
therefore neglected, as in Ref. Magnusson  _et al._ (2022). However, for
droplets traveling at small interfacial separations $s$ the angular dynamics
must be considered, because hydrodynamic torques tend to cause one droplet to
roll over the other one.

###  II.4 Numerical integration

To make sure that droplets do not overlap in the numerical simulations when
they travel at very small interfacial distances for finite times, we followed
Ref. Li Sing How (2021) and employed a $4$-th order Runge-Kutta algorithm with
an adaptive time step. Our prescription for the adaptive time step differs
slightly from the one in Ref. Li Sing How (2021). We used $\Delta
t=\varepsilon\min\left\\{{s}/{V_{R}},{s}/{V_{T}}\right\\}$. Here, $V_{R}$ and
$V_{T}$ are the radial and tangential speeds, and $\varepsilon\ll 1$ is a
numerical constant. The numerical results shown in Fig. 2 were obtained for
$\varepsilon=0.08$. We verified that reducing the numerical constant to
$\varepsilon=0.001$ did not change the results shown in Fig. 2(d).

We recorded a collision and stopped the numerical integration as soon as the
interfacial distance became smaller than $10^{-14}$ m. We verified that using
a larger cutoff, $10^{-13}$, did not change the results shown in Fig. 2(d).

##  III Results

\begin{overpic}[scale={0.55}]{Fig2.pdf} \end{overpic} Figure 2:  Collision
efficiencies and phase portraits of charged droplets settling in still air.
(a) Collision efficiency $\mathcal{E}$ for droplets with equal amounts of
charge of opposite parity, as a function of charge $q=|q_{1,2}|$ over
elementary charge $e$. Droplet radii: $a_{1}=16\,\mu$m and $a_{2}=20\,\mu$m
(blue), and $a_{1}=8\,\mu$m and $a_{2}=10\,\mu$m (red). Kn $=10^{-3}$ (solid
lines), $10^{-2}$ (dashed), and $5\times 10^{-2}$ (dash-dotted). The vertical
solid lines indicate the locations where Eqs. (14) and (15) predict a
transition. Arrows denote the charges in panels (b,c). Panel (b) shows the
relative droplet dynamics in the $R_{1}$-$R_{3}$-plane (where $R_{i}$ is non-
dimensionalised by $\bar{a}$), for $a_{1}=16\,\mu$m, $a_{2}=20\,\mu$m,
$q_{1,2}=\pm 908\,e$, and Kn $=10^{-3}$. Colliding trajectories (blue), non-
colliding trajectories (red), The solid black line is the separatrix between
colliding and non-colliding trajectories, and $\bullet$ represents the saddle
point S$1$. (c) Same as (b), but for $q_{1,2}=\pm 15469\,e$. Manifolds of the
saddle point S$1$ ($\bullet$) are shown as black solid lines. (d) Same as (a),
but for equal parity. Panel (e) shows the bifurcations (their locations are
denoted by vertical dashed lines) that give rise to saddle points S$2$ and
S$3$. Panel (f) shows the relative dynamics in the $R_{1}$-$R_{3}$-plane for
$a_{1}=16\,\mu$m, $a_{2}=20\,\mu$m, $q_{1,2}=6011\,e$ along with the saddle
point S$1$, and panel (g) corresponds to $q_{1,2}=21199\,e$, and shows the
saddle point S$2$.

In the following we discuss how the collision efficiency depends on the excess
droplet charge. The collision efficiency $\mathscr{E}$ is defined in terms of
the critical impact parameter $b_{\rm c}$ beyond which the droplets cease to
collide, $\mathscr{E}=b_{\rm c}^{2}/(a_{1}+a_{2})^{2}$. For droplets charged
with equal and with opposite polarities, we describe a transition from a
regime at small charges that is dominated by non-continuum effects, to a
charge-dominated regime at large excess charges. In both cases, the collision
efficiency is determined by a saddle point below the collision sphere, and the
transition between the regimes occurs when this saddle point leaves the region
where non-continuum effects dominate [Fig. 1(b)]. Our results for the
collision efficiency are shown in Fig. 2 for droplets that carry the same
amount of charge with different polarities (left), or with the same polarity
(right). We discuss these two cases separately in Sections III.1 and III.2.

The droplets were initialised sufficiently far apart that they settled
independently at their steady Stokes settling velocities. The initial angular
velocities were set to zero. Since the settling velocities and the centre-of-
mass separation $R$ were initially in the $R_{1}$-$R_{3}$-plane, and since the
torques and forces depend upon the cross product of the translational and
angular velocities, respectively, with the separation vector Kim and Karrila
(1991), the spatial dynamics of the droplets must remain confined to the
$R_{1}-R_{3}$ plane. The only non-zero component of the angular velocities was
$\omega_{2}^{(i)}$. Thus, due to the choice of initial conditions, one
component of the droplet velocities and two components of angular velocities
remain zero, $v_{2}^{(i)}=0,\,\omega_{1}^{(i)}=0,$ and $\omega_{3}^{(i)}=0$.

###  III.1 Opposite charges

Figure 2(a) shows that electrostatic attraction increases the collision
efficiencies $\mathscr{E}$ for droplet pairs of radii $8\,\mu$m and $10\,\mu$m
(red) more strongly than for droplets with radii $16\,\mu$m and $20\,\mu$m
(blue). The droplets were charged with the same amount of excess charge, but
with different polarities. For weak charges, $\mathscr{E}$ is approximately
independent of charge, but the collision efficiency depends on the Knudsen
number, Kn. In this limit, the collision dynamics resembles that of neutral
droplets settling in still air, as seen in Fig. 2(a) in Ref. Magnusson  _et
al._ (2022) (but note that there are subtle differences that we discuss
below). For strong charges, by contrast, the collision efficiency is
independent of Kn, as predicted by Magnusson et al. Magnusson  _et al._
(2022). The difference between these two regimes can be understood by
following the collision dynamics in the $R_{1}$-$R_{3}$-plane [Fig. 1(a)],
corresponding to a frame that translates with the smaller droplet.

Consider first the limit of large charges, panel (c). The collision efficiency
is determined by the stable manifold of a saddle point S$1$ located below the
collision sphere. The unstable manifold of the saddle is one-dimensional, with
a component along the $R_{3}-$axis. The stable manifold therefore has co-
dimension one. As a consequence, the stable manifold forms a separatrix for
the phase-space dynamics Wiggins (2013). The corresponding solid black line in
panel (c) represents the projection of a one-dimensional curve within the
stable manifold, the curve defined by the chosen initial condition.
Trajectories from above that approach the collision sphere above this line
must collide (blue). All other trajectories from above do not collide (red).
The saddle point S$1$ has four complex eigenvalues. The corresponding
eigenvectors span the $R_{1}-R_{2}-V_{1}-V_{2}$ space. This means that
trajectories approaching the fixed point in this plane spiral around the fixed
point in phase space. The trajectories plotted in black in panel (c) do
slightly overshoot the saddle point S$1$ , but this is not visible in the
Figure. Since the saddle point and its stable manifold are far from the
collision sphere, the local breakdown of the hydrodynamic approximation does
not matter. As a consequence, the collision efficiency is independent of Kn.
This is the limit analysed by Magnusson et al. Magnusson  _et al._ (2022), see
their Fig. 2(b).

The limit of small charges, shown in panel (b), is more subtle than the large-
charge limit. As mentioned above, the collision dynamics looks like that for
neutral droplets, where the separatrices delineating collisions from no
collisions (solid black lines) are grazing trajectories. However, our numerics
indicates that separatrices shown in panel (b) do not graze the collision
sphere. They lie very close to the collision sphere, but they appear to
connect to a saddle point S$1$ below the collision sphere, located on the
$R_{3}$-axis. The numerical integration becomes difficult when the trajectory
travels along the collision sphere at interfacial separations $s<{\rm Kn}$.
Therefore we could not resolve the saddle point in the numerical integration
(the trajectories corresponding to the black solid lines do not hit the saddle
point, but they pass very close to it). As outlined above, a transition
between the two different regimes shown in Fig. 1(b) occurs when the saddle
point S$1$ exits the region $s<{\rm Kn}$. In order to compute the charge at
which the transition occurs, we first start by determining the location of the
saddle points of the relative dynamics. In order to find the saddle points of
the relative droplet dynamics, consider the dynamical system defined by
$\dot{\mbox{\boldmath$R$}}=\mbox{\boldmath$v$}^{(2)}-\mbox{\boldmath$v$}^{(1)}$,
Eqs. (1b) and (1c). We look for zeros of this dynamical system,
$\dot{\mbox{\boldmath$R$}}=0,\dot{\mbox{\boldmath$v$}}^{(i)}=0$ and
$\dot{\mbox{\boldmath$\omega$}}^{(i)}=0$, $i=1,2$. At the saddle point, the
droplets fall at a constant interfacial separation with a common steady-state
settling velocity
${\mbox{\boldmath$v$}}^{(1)}={\mbox{\boldmath$v$}}^{(2)}=\mbox{\boldmath$v$}_{\rm
s}$. Due to the symmetry of the problem, $\mbox{\boldmath$v$}_{\rm s}$ must
point along the direction of gravity. Since there is no external torque acting
on the droplets, and because of the dissipative nature of the dynamics, the
angular velocities at the fixed point must vanish,
${\mbox{\boldmath$\omega$}}^{(i)}=0$. Thus, we find the location of this
equilibrium point by requiring that Magnusson  _et al._ (2022)

$\mbox{\boldmath$F$}^{(i)}_{\rm g}+\mbox{\boldmath$F$}^{(i)}_{\rm
h}+\mbox{\boldmath$F$}^{(i)}_{\rm e}=0\quad\mbox{for}\quad i=1,2\,.$ (11)

Given the forces described in Section II, Eq. (11) simplifies to

$\displaystyle-
m_{1}g-6\pi\mu\left(a_{1}X^{A}_{11}+\overline{a}X^{A}_{12}\right)v_{\rm
s}+{F}_{\rm e}$ $\displaystyle=0\,,$ (12a) $\displaystyle-
m_{2}g-6\pi\mu\left(a_{2}X^{A}_{22}+\overline{a}X^{A}_{21}\right)v_{\rm
s}-{F}_{\rm e}$ $\displaystyle=0\,.$ (12b)

Here $\mu=\varrho_{\rm f}\nu$ is the dynamic viscosity of air, and $v_{\rm
s}\equiv|\mbox{\boldmath$v$}_{\rm s}|$ is the common settling speed of the two
droplets in the steady state, and $F_{\rm e}=F_{\rm e}^{(1)}$. For given
droplet radii and charges, solutions to these equations give the saddle-point
location $[0,R_{3}^{\ast}]$, and centre-of-mass settling speed $V_{\rm s}$.
For the parameters corresponding to panel (b), we find $s^{\ast}\approx
6\times 10^{-4}$ . So the equilibrium point is very close to, but not on the
collision sphere. As mentioned above, we could not find the stable and
unstable manifolds by numerical integration, because numerical integration is
hard when the trajectories travel at interfacial separations $s<{\rm Kn}$
parallel to the collision sphere. However, the dynamics in the vicinity of the
equilibrium point [panel (b)] indicates that it is a saddle: the fixed point
has eigenvalues with positive and negative real parts. As the charge tends to
zero, the saddle approaches the collision sphere, $R_{3}^{\ast}\to-2$.
Conversely, as the charge increases, the saddle moves down the $R_{3}$-axis.
We conclude that the saddle point identified in Ref. Magnusson  _et al._
(2022) appears already at small Coulomb numbers, ${\rm Cu}\ll 1$. Since the
numerical integration is more difficult at smaller (but non-zero) charges, we
did not investigate this limit by numerical integration.

How much charge is needed to cross over from the small charge limit where non-
continuum effects determine collisions to the charge-dominated limit? There is
no bifurcation, the same saddle determines the collision dynamics in both
regimes. But the above considerations indicate that the transition occurs when
the saddle point leaves the region $s<{\rm Kn}$ in the $R_{1}$-$R_{3}$-plane,
the hashed region in Fig. 1(b). If the saddle point is close to the collision
sphere, its location and its manifolds depend on Kn. By contrast, when the
interfacial separation at the saddle point is much larger than the mean free
path $\ell$, the saddle point and its associated separatrix are independent of
Kn, and the collision efficiency is determined by electrostatic forces.

To find the critical charge, we set the interfacial separation equal to the
mean free path $\ell$ in Eqs. (12), or $s={\rm Kn}$ in non-dimensional units.
Substituting $s={\rm Kn}$ into Eqs. (12) and solving for $V_{\rm s}$ and
$F_{\rm e}$ we find the electrostatic force required to maintain a time-
independent relative settling velocity $v_{\rm s}$,

$\displaystyle v_{\rm s}$ $\displaystyle=\frac{4a_{1}^{2}g\varrho_{\rm
w}(1+\lambda^{3})}{9\mu\lambda^{2}}\frac{1}{2\lambda X^{A,{\rm
uv}}_{11}+(1+\lambda)(X^{A,{\rm uv}}_{12}+X^{A,{\rm uv}}_{21})+2X^{A,{\rm
uv}}_{22}}\,,$ (13a) $\displaystyle F_{\rm e}$
$\displaystyle=\frac{4a_{1}^{3}g\pi\varrho_{\rm
w}}{3\lambda^{3}}\frac{-2\lambda X^{A,{\rm uv}}_{11}-(1+\lambda)X^{A,{\rm
uv}}_{12}+\lambda^{3}[(1+\lambda)X^{A,{\rm uv}}_{21}+2X^{A,{\rm
uv}}_{22}]}{2\lambda X^{A,{\rm uv}}_{11}+(1+\lambda)(X^{A,{\rm
uv}}_{12}+X^{A,{\rm uv}}_{21}{)}+2X^{A,{\rm uv}}_{22}}\,.$ (13b)

To determine the corresponding amount of charge, we approximate $F_{\rm e}$
and the resistance functions by their small-$s$ asymptotes (Section II). We
find

$\displaystyle|q_{1}q_{2}|=\frac{1}{h(\lambda,{\lambda_{q}})}\frac{a_{1}^{5}g\,{\rm
Kn}\rho_{p}}{k_{e}\lambda}\left[f^{(1)}(\lambda)\log\left(\frac{4\lambda}{(1+\lambda)^{2}{\rm
Kn}}\right)^{2}+f^{(2)}(\lambda)\log\left(\frac{4\lambda}{(1+\lambda)^{2}{\rm
Kn}}\right)+f^{(3)}(\lambda)\right].$ (14)

In this expression, the function $h(\lambda,\lambda_{q})$ is defined as

$\displaystyle h(\lambda,{\lambda_{q}})$
$\displaystyle=\frac{1}{{|\lambda_{q}}|}\left[\left(\psi(\tfrac{1}{1+\lambda})+\gamma\right)-\left(\psi(\tfrac{\lambda}{1+\lambda})+\gamma\right){\lambda_{q}}\right]^{2}.$
(15a) The expressions for the functions $f^{(k)}$ are quite lengthy. Here we
only quote their Padé approximants around $\lambda=1$, of order $[3,3]$ (the
transition values shown in Fig. 2 were obtained using the exact functions),
$\displaystyle f^{(1)}$
$\displaystyle=\frac{8.98831(\lambda-1)-12.553(\lambda-1)^{2}+4.18158(\lambda-1)^{3}}{1-3.39659(\lambda-1)+3.95407(\lambda-1)^{2}-1.59769(\lambda-1)^{3}}\,,$
(15b) $\displaystyle f^{(2)}$
$\displaystyle=\frac{-45.6737(\lambda-1)+2.16728(\lambda-1)^{2}-0.229416(\lambda-1)^{3}}{1+1.95255(\lambda-1)+0.830927(\lambda-1)^{2}-0.112841(\lambda-1)^{3}}\,,$
(15c) $\displaystyle f^{(3)}$
$\displaystyle=\frac{-58.0221(\lambda-1)-5.79319(\lambda-1)^{2}-2.67774(\lambda-1)^{3}}{1+2.09984(\lambda-1)+1.39178(\lambda-1)^{2}+0.352657(\lambda-1)^{3}}\,.$
(15d)

Fig. 2 shows the predicted charge $q=|q_{1,2}|$ (assuming $\lambda_{q}=\pm 1$)
computed using Eqs. (14) and (15) as vertical solid lines. We see that the
theory predicts the transition from the small-charge limit where non-continuum
effects determine collisions to the charge-dominated regime for both droplet
sizes considered. This confirms the theoretical prediction that the transition
occurs when the saddle point leaves the Kn-dominated region in the
$R_{1}$-$R_{3}$-plane.

###  III.2 Same charges

Panel (d) illustrates that the collision efficiencies $\mathscr{E}$ for
droplet pairs with equal amounts of excess positive charge tend to zero as the
charge increases. Coulomb repulsion prevents collisions entirely at
sufficiently large charges. Furthermore, the collision efficiency vanishes
more quickly for the droplet pair with radii $8\,\mu$m and $10\,\mu$m than for
the droplet pair with radii $16\,\mu$m and $20\,\mu$m, as expected, because
the Coulomb number is larger for smaller droplets with the same charge. As in
panel (a), we identify two regimes. For small charges, the collision
efficiency depends on Kn but the Kn-dependence becomes much weaker as the
charge magnitude increases and the collision efficiency becomes charge
dominated. The transition mechanism is the same as for oppositely charged
droplets: despite the fact that the droplets carry charges with the same
polarity, the electrostatic force is attractive at small separations Lekner
(2012). This gives rise to a saddle-point S$1$ below the collision sphere, as
in panel (b), where the attractive force balances differential settling. The
corresponding predictions of Eq. (14) are shown as vertical solid lines in
panel (d). One difference to the collision dynamics for charged droplets with
equal polarity [panel (a)] is that the crossover between the small-charge
regime where non-continuum effects determine collisions and the charge-
dominated regime occurs at much larger charges in panel (d). The reason for
this difference is that the small-$s$ attractive force between droplets with
the same charge is weaker than the attractive force for droplets with opposite
charges. Moreover, the electrostatic force for droplets with the same charge
is repulsive at large interfacial distances. Since the force attracts at small
separations, it must change sign at a non-dimensional critical distance
$s_{\rm c}$ that depends on $\lambda$ and $\lambda_{q}$. Therefore, the
droplet interfacial separation at the saddle point must remain smaller than
this critical distance, $s_{\rm c}$, as the charge increases. For the
parameters used here, $s_{\rm c}\approx 0.018$. This means that the collision
efficiency for droplets with the same polarity becomes independent of the Kn
number only when ${\rm Kn}<0.018$. There is no such condition for droplets
charged with opposite charges, because the saddle point continues to move down
as the amount of charge increases (Section III.1).

As the magnitude of the excess charge increases, a bifurcation occurs above
the collision sphere. This bifurcation gives rise to two new saddle points,
S$2$ and S$3$. The saddle-point S$2$ results from the balance between
electrostatic repulsion and differential settling. The stable manifold of S$2$
is $13$-dimensional. Its unstable manifold is two-dimensional, therefore it
does not form a separatrix for the collision dynamics. Nevertheless, the
collision efficiency appears to vanish already before the bifurcation for
reasons discussed below. In addition, droplets which approach the fixed point
along the $R_{3}-$ axis spend a long time in each others vicinity. The stable
manifold of S$3$ is $12$-dimensional, and its unstable manifold has dimension
three. This saddle point is not approached by trajectories starting at large
separations for the Stokes numbers considered here. However, at much larger
Stokes numbers, trajectories could overshoot the saddle S$2$ and approach the
fixed point S$3$.

We determined the charge $q_{\rm c}$ at which the bifurcation occurs by
numerically solving Eq. (12) for the interfacial separation $s$ and the common
settling speed $v_{\rm s}$. We determined the location of the bifurcation by
computing the value of the droplet charge for which Eq. (12) had exactly one
solution for $R_{3}>0$. The locations of the two saddles on the $R_{3}$-axis,
for $R_{3}>0$, are shown in panel (e). Panel (f) shows the relative dynamics
in the $R_{1}$-$R_{3}$-plane before bifurcation. It is qualitatively similar
to panel (b). At weak charges, it does not matter whether the droplets have
the same or opposite polarities. For both cases, a saddle point below results
due to the attractive force at close approach.

Panel (g) shows the trajectories in the $R_{1}$-$R_{3}$-plane after the
bifurcation. The flow in the vicinity of S$2$ confirms that this fixed point
is a saddle point. The second saddle point above, S$3$, and the saddle point
S$1$ below the collision sphere are not shown, because these fixed points are
not approached by trajectories approaching from afar. While the unstable
manifolds of the saddle point S$2$ do not form separatrices that prevent
collisions from above, at the droplet sizes and initial conditions considered,
we do not observe any collisions for charges above the bifurcation values. The
bifurcation values for the charges are shown in panel (d) as dashed vertical
lines. However, we also observe that $\mathscr{E}$ tends to become vanishingly
small well before the bifurcation occurs. The reason is that even before the
saddle point appears, the dynamics in the $R_{3}-$direction slows down in its
vicinity. These slow trajectories are then deflected in the unstable $R_{1}$
direction so that only trajectories approaching infinitesimally close to the
$R_{3}$ axis can collide. We speculate that the qualitative form of the
collision rate before the bifurcation can be understood by studying the normal
form of the bifurcation for our system.

##  IV Discussion

The amount of excess charge required to reach the Coulomb-dominated collision
regime depends on the polarity, upon whether the two droplets have excess
charges of the same sign, or not. For droplets carrying the same charges,
larger amounts of charge are required to qualitatively change the collision
dynamics, about $q\sim 10^{4}\,e$ for $20\,\mu$m-droplets. For droplets that
are charged with opposite polarities, much smaller charges are required, about
$10^{3}\,e$ for $20\,\mu$m-droplets. In this case, the critical charge is much
closer to the estimate of Davis (1965) than that of Tinsley and Zhou (2014).
However, these earlier estimates did not consider how the hydrodynamic and
non-continuum forces vary at small interfacial distances. While Davis (1965)
included hydrodynamic forces valid up to $10^{-3}a_{2}$ (where $a_{2}$ is the
radius of the larger droplet), he did not account for non-continuum effects.
As a consequence the results for the collision efficiency depend on an
arbitrary cutoff assumed to define a collision Dhanasekaran  _et al._ (2021).

Tinsley and Zhou (2014) state that at least $10^{4}\,e$ are required for
electrical forces to significantly affect the collision efficiency of
oppositely charged droplets with radii of the order of $20$ microns, and that
for similarly charged droplets this threshold is even higher. Our results
indicate, by contrast, that the threshold is much smaller, about $10^{3}\,e$
for oppositely charged droplets, and about $10^{4}\,e$ for droplets with the
same polarity. These differences are not surprising, given that early
calculations Paluch (1970); Schlamp  _et al._ (1976, 1979); Tinsley and Zhou
(2014) used highly idealised models that neither account for how the
hydrodynamic forces change as the droplets come close, nor how hydrodynamics
breaks down below the mean free path, causing their predictions to fail Abbott
(1974).

Eq. (14) shows how the critical charge $q_{c}$ depends upon the parameters of
the problem. We see that $|q_{1}q_{2}|\sim a_{1}^{5}$. This means that the
critical charge $q_{c}$ decreases rapidly as the droplet size decreases.
Keeping all parameters the same as in Fig. 2, but with $a_{2}=4\,\mu$m and
$a_{2}=5\,\mu$m and Kn $=0.001$, we see that a charge of only $30\,e$ is
enough to lead to a transition for oppositely charged droplets. Note, however,
that spatial diffusion (not considered here) must affect the relative dynamics
of such small droplets.

The model includes non-continuum effects due to the breakdown of hydrodynamics
at small interfacial distances for both radial and tangential resistance
functions. Ref. Li Sing How (2021) demonstrated that at small Stokes numbers,
the collision efficiencies match those obtained for a St = 0 model from Ref.
Dhanasekaran  _et al._ (2021), which did not include non-continuum corrections
to the tangential resistance functions. However, the non-continuum corrections
to tangential resistance functions might become important at larger St
numbers, or when the radius ratio of the droplets deviates significantly from
unity Li Sing How  _et al._ (2021).

Our model has five non-dimensional parameters: the Stokes number St (particle
inertia), the Coulomb number Cu (charge), the ratio $\lambda$ of droplet
radii, the charge ratio $\lambda_{q}$, and the Knudsen number Kn (effect of
mean free path). We have not yet systematically studied the effect of changing
St, $\lambda$, and $\lambda_{q}$. Instead we chose $\lambda=0.8$,
$\lambda_{q}={\pm}1$, and two values of the Stokes number, ${\rm St}=0.42$,
and $3.96$. Since the droplet-size distribution in clouds at the onset of
collisional growth is narrow Grabowski and Wang (2013), it is natural to
consider radius ratios close to unity, as we did here. As a next step, it is
of interest to study the effect of the parameter $\lambda_{q}$. While the
charge distribution of cloud droplets is not known in general, extreme values
of $\lambda_{q}$ can significantly enhance collision rates Khain  _et al._
(2004). Therefore it is important to study how the collision efficiency
changes as the charge ratio $\lambda_{q}$ varies.

Let us briefly comment on the limitations of the model. First, we did not
consider van-der-Waals forces, short-range attractive dipole forces that
change the collision dynamics at short interfacial distances $s$, of the order
of the London length $\lambda_{L}$ Zhang and Davis (1991). Rother  _et al._
(2022) computed the collision rate of droplets settling in still air including
van-der-Waals forces and Maxwell slip (the first-order correction to the non-
continuum effects). They found that van-der-Waals forces cause a significant
increase in the collision rate for droplets smaller than $10\,\mu$m in radius,
but that the effect is small for larger droplets. This is consistent with the
findings of Dhanasekaran  _et al._ (2021), who stated that van-der-Waals
forces do not matter, compared to non-continuum effects, for water droplets
with radii larger than $10\,\mu$m. In order to understand how van-der-Waals
forces affect the collision dynamics of charged droplets, we included the van-
der-Waals force in our simulations, using Ho and Higuchi (1968); Zhang and
Davis (1991)

$\displaystyle F_{\rm vdW}^{(i)}=-{A_{\rm
H}}\begin{cases}\frac{{a_{1}a_{2}}}{a_{1}+a_{2}}\frac{\lambda_{L}(\lambda_{L}+7.0768\pi
s)}{s^{2}(\lambda_{L}+2\pi(1.7692)s)^{2}}&\quad\mbox{for}\quad
s<0.25\,\lambda_{L}/\pi\,,\\\
\frac{{a_{1}a_{2}}}{a_{1}+a_{2}}\left(\frac{4.9\lambda_{L}}{60\pi
s^{3}}-\frac{6.51\lambda_{L}^{2}}{360\pi^{2}s^{4}}+\frac{2.36\lambda_{L}^{3}}{1680\pi^{3}s^{5}}\right)&\quad\mbox{for}\quad
s\geq 0.25\,\lambda_{L}/\pi\,.\\\ \end{cases}$ (16)

Here $A_{H}$ is the Hamaker constant, with typical value $A_{H}=5\times
10^{-20}$ J for water droplets in air. Two new non-dimensional parameters come
with Eq. (16): $\lambda_{L}/{\overline{a}}$ quantifies the interfacial
separation below which van-der-Waals forces become important, while
$A_{H}/(\tilde{m}V_{0}^{2})$ quantifies the strength of van-der-Waals forces
with respect to the relative kinetic energy at large separations.

We found that including van-der-Waals forces in our simulations leads to a
small increase of about 15% in the collision efficiency for Kn$=0.01$ and
droplets of sizes $16\,\mu$m and $20\,\mu$m, but does not cause qualitative
changes to Figs. 2(a,d). This is consistent with the findings of Refs.
Dhanasekaran  _et al._ (2021); Rother  _et al._ (2022). In the limit of small
charges and small Kn, by contrast, van-der-Waals forces may change of the
collision dynamics qualitatively. We expect van-der-Waals force gives rise to
a saddle point below the collision sphere, even for neutral droplets. In this
case, the stable manifold of this new fixed point may determine the collision
efficiency at small values of Kn.

Second, our model uses the Stokes approximation for the hydrodynamic forces
and torques. The effects of convective and unsteady fluid inertia were
neglected. Klett and Davis (1973) explained that convective fluid inertia
makes a significant difference to the collision efficiency of droplets with
radius ratios close to unity, because the droplet dynamics in the Stokes
approximation is degenerate when $a_{1}=a_{2}$ (see also Ref. Candelier and
Mehlig (2016)). Magnusson  _et al._ (2022) considered the leading-order
effects of the particle Reynolds number and Strouhal number for droplets with
${\rm Re}_{p}\approx 0.1$ and ${\rm Sl}\approx 0.1$, at large droplet
separations. Here we do not consider these effects because we do not know how
to account for the combined effects of the Reynolds and Strouhal effects at
small interfacial distances.

Third, we assumed that the droplets remain spherical before they collide.
Droplet deformation during collisions is controlled by the Weber number We,
defined as the ratio of the droplet kinetic energy upon contact to its surface
energy Qian and Law (1997). When the Weber number is much smaller than unity,
as for a $20\,\mu$m droplet settling at its terminal velocity (We $\sim
10^{-6}$), droplet deformation is negligible. At much larger collision speeds
the Weber number is larger. In this case, and the droplets are expected to
deform. This changes not only the electrostatic forces, but also the
hydrodynamic forces and their non-continuum regularisations, effects not
considered here.

Fourth, we assumed that the droplets are perfect conductors. The time scale
$\tau$ of charge redistribution on a water droplet depends on the ratio of the
permittivity of water $\varepsilon_{\rm w}$ to its electrical conductivity
$\sigma$, $\tau=\varepsilon_{\rm w}/\sigma$ Haus and Melcher (1989). With the
conservative estimate $\sigma\sim 5\times 10^{-6}$ S/m for pure water, and
$\varepsilon_{w}\sim 7\times 10^{-10}$ C2/Nm Rumble (2022), one finds
$\tau\sim 1.4\times 10^{-4}$ s. This time is much shorter than the particle-
relaxation time for relative dynamics,
$\tau_{\gamma}=|\gamma_{2}^{-1}-\gamma_{1}^{-1}|$. For a pair of $16\,\mu$m-
and $20\,\mu$m-droplets, one finds $\tau/\tau_{\gamma}\sim 0.08$. The small
but non-zero charge-relaxation may still have a small effect on the
electrostatic forces which we neglected.

Fifth, we considered special initial conditions. We assumed that the droplets
are initially so far apart that they settle independently with their
respective Stokes settling speeds, and we set their initial angular velocities
to zero. As a result, the separatrices shown in Fig. 2 correspond to
intersections between the invariant manifolds and the one-parameter family of
curves determined by the initial conditions.

Sixth, we considered droplets settling in still air. What is the effect of
turbulence on the collision efficiency of charged droplets? Saffman and Turner
(1956) explained how turbulent strains increase collision rates of cloud
droplets, neglecting their interactions. Dhanasekaran  _et al._ (2021)
analysed how hydrodynamic interactions and their regularisation due to the
breakdown of hydrodynamics affect the collision dynamics of droplets settling
in a steady straining flow. They found much smaller collision rates than
Saffman & Turner, and that hydrodynamic interactions lead to an intricate
dependence of the collision rate on the straining flow. This behaviour is
explained by a sequence of bifurcations, both bifurcations of equilibria and
grazing bifurcations Dubey  _et al._ (2022). It remains an open question how
electrostatic interactions change this picture. Moreover, the close approach
of charged droplets in turbulence was measured in Refs. Lu  _et al._ (2010a,
b); Lu and Shaw (2015), and compared with a theory valid for large droplet
separations – similar to the model described by Magnusson  _et al._ (2022). In
order to quantitatively describe the relative dynamics of nearby droplets, we
intend to generalise the model to arbitrary linear flows. In the presence of a
background flow, convective-inertia effects may give rise to inertial lift
forces Candelier  _et al._ (2019).

##  V Conclusions

We analysed the collision dynamics of weakly charged, $\mu$m-sized water
droplets settling in still air by numerical integration of a model that
incorporates particle inertia, hydrodynamic interactions in the Stokes
approximation, non-continuum effects, and electrical forces. We observed two
distinct collision regimes. For small charges, the collision dynamics is
dominated by short-range hydrodynamic interactions. In this regime it is
crucial to take into account how these interactions are regularised below the
mean free path of air. As a consequence, the collision efficiency depends on
the Knudsen number, Kn, in this regime. At large charges, in contrast, the
collision efficiency does not depend on Kn, because the separatrix between
colliding and non-colliding trajectories is the stable manifold of a saddle
point far from the collision sphere, consistent with the conclusions drawn in
Magnusson et al. Magnusson  _et al._ (2022).

The two regimes shown in Fig. 1(b) occur both for droplet pairs with excess
charges of the same and of different polarities. This may be surprising at
first sight, considering that equal point charges repel each other, while
opposite point charges attract. The reason is that induced charges cause the
electrostatic force to always be attractive, at small enough separations
Lekner (2012). Ref. Khain  _et al._ (2004) describes how induced electrical
charges enhance the collision rate between a charged droplet and neutral
droplet, where the effect of induced charges is apparent at large separations
(where their model applies). For micron-sized droplets of similar sizes, the
electrostatic force changes sign at an interfacial distance much smaller than
the droplet radius. But since this distance is of the order of the mean free
path, where the hydrodynamic approximation breaks down, the induced charges
have a significant effect on the collision efficiency, as we demonstrated in
this paper.

We found that the saddle point S$1$ (Fig. 2) below the collision sphere exists
even at small charges. But when the interfacial distance at the saddle-point
$s<{\rm Kn}$, its stable manifolds depend on the Knudsen number. As a result
the collision efficiency depends on Kn, at small charges. The crossover
between the two limits occurs when the saddle point moves further away from
the collision sphere than the mean free path. The critical charge where this
qualitative change in the collision dynamics occurs is much lower than that
stated by Tinsley and Zhou (2014). While our prediction matches that of Davis
(1965), they did not include non-continuum effects at small separations, so
that their result depends upon the cutoff used to define a collision. We
observed that for $20\,\mu$m-droplets, the transition between the two regimes
happens at $\sim 10^{3}\,e$ for droplets charged with opposite polarities, and
at $\sim 10^{4}\,e$ for droplets charged with the same polarities. These
charge magnitudes are comparable to those observed for thunder-cloud droplets
Takahashi (1973). Note, however, that in thunder clouds, electric fields
affect the droplet dynamics. Here we did not consider this effect.

These conclusions are based on a model with five non-dimensional parameters,
and we have only analysed a small part of the parameter space. Including van-
der-Waals forces give rise to two additional dimensionless parameters Zhang
and Davis (1991); Rother  _et al._ (2022), one quantifying the interfacial
separation below which the van-der-Waals force is important, and another
quantifying its strength. While previous studies Dhanasekaran  _et al._
(2021); Rother  _et al._ (2022) suggest that including van-der-Waals force
does not give rise to a qualitative changes in the dynamics for the parameters
typical of cloud droplets with radii larger than $10\,\mu$m, we speculate that
the qualitative dynamics are different when the Kn number becomes vanishingly
small. The reason is that the attractive van-der-Waal force gives rise to a
saddle point below, which determines the relative dynamics for vanishingly
small Kn numbers.

Finally, we assumed here that the droplets settle in quiescent air. We did not
consider the effects of turbulent flow. Turbulence can substantially increase
the rate at which droplets collide Saffman and Turner (1956); Dhanasekaran
_et al._ (2021), because turbulent strains bring similar-sized droplets
together. More generally, turbulence increases spatial clustering of droplets,
as measured by the pair-correlation function (see Ref. Gustavsson and Mehlig
(2016) for a review). Lu et al. Lu  _et al._ (2010a, b); Lu and Shaw (2015)
measured how electrical charges change the pair-correlation function, and
developed models to explain the observed effects – at separations much larger
than the droplet radii where hydrodynamic interactions can be neglected. It
would be of interest to include in the model the mechanisms described here, in
order to determine how charges affect the pair-correlation function at smaller
separations. Similarly, the relative-velocity statistics of charged droplets
in turbulence was investigated in Ref. Alipchenkov  _et al._ (2004),
neglecting hydrodynamic interactions, but our results indicate that these
interactions (and their regularisation by non-continuum effects) cannot be
neglected at small separations. In summary, the phase-space picture presented
here shows that fairly small charges can have a significant effect on the
relative dynamics of nearby droplets. How this affects the pair-correlation
function, relative-velocity statistics, and collision efficiencies in
turbulence remains an open question.

Acknowledgements. The research of AD and BM was supported by grants from
Vetenskapsrådet (grant nos. 2017-3865 and 2021-4452 and from the Knut and
Alice Wallenberg Foundation (grant no. 2014.0048). The research of GB was
supported by the National Science Foundation under grant no. CBET-1605195. BM
acknowledges support from a Mary Shepard B. Upson Visiting Professorship with
the Sibley School of Mechanical and Aerospace Engineering at Cornell.

## References

  * Pruppacher and Klett (1997) H. R. Pruppacher and J. D. Klett,  _Microphysics of clouds and precipitation, 2nd edition_ (Kluwer Academic Publishers, Dordrecht, The Nederlands, 1997) 954p.
  * Shaw (2003) R. A. Shaw, “Particle-turbulence interactions in atmospheric clouds,” Annu. Rev. Fluid Mech. 35, 183–227 (2003).
  * Devenish  _et al._ (2012) B. J. Devenish, P. Bartello, J. L. Brenguier, L. R. Collins, W. W. Grabowski, R. H.A. Ijzermans, S. P. Malinowski, M. W. Reeks, J. C. Vassilicos, L. P. Wang, and Z. Warhaft, “Droplet growth in warm turbulent clouds,” Q. J. R. Meteorol. Soc. 138, 1401–1429 (2012).
  * Berry and Reinhardt (1974) E. X. Berry and R. L. Reinhardt, “An analysis of cloud drop growth by collection: Part I. Double distributions,” Journal of Atmospheric Sciences 31, 1814 – 1824 (1974).
  * Grabowski and Wang (2013) W. W. Grabowski and L.-P. Wang, “Growth of cloud droplets in a turbulent environment,” Annu. Rev. Fluid Mech. 45, 293–324 (2013).
  * Saffman and Turner (1956) P. G. Saffman and J. S. Turner, “On the collision of drops in turbulent clouds,” J. Fluid Mech. 1, 16–30 (1956).
  * Kostinski and Shaw (2005) A. B. Kostinski and R. A. Shaw, “Fluctuations and luck in droplet growth by coalescence,” Bulletin of the American Meteorological Society 86, 235–244 (2005).
  * Kim and Karrila (1991) S. Kim and S. J. Karrila,  _Microhydrodynamics: principles and selected applications_ (Butterworth-Heinemann, Boston, 1991).
  * Sundararajakumar and Koch (1996) R. R. Sundararajakumar and D. L. Koch, “Non-continuum lubrication flows between particles colliding in a gas,” J. Fluid Mech. 313, 283–308 (1996).
  * Dhanasekaran  _et al._ (2021) J. Dhanasekaran, A. Roy, and D. L. Koch, “Collision rate of bidisperse spheres settling in a compressional non-continuum gas flow,” J. Fluid Mech. 910, A10 (2021).
  * Dubey  _et al._ (2022) A. Dubey, K. Gustavsson, G. P. Bewley, and B. Mehlig, “Bifurcations in droplet collisions,” Phys. Rev. Fluids 7, 064401 (2022).
  * Sundaram and Collins (1997) S. Sundaram and L. R. Collins, “Collision statistics in an isotropic particle-laden turbulent suspension,” J. Fluid. Mech. 335, 75–109 (1997).
  * Falkovich  _et al._ (2002) G. Falkovich, I. Fouxon, and M. G. Stepanov, “Acceleration of rain initiation by cloud turbulence,” Nature 419, 151–154 (2002).
  * Gustavsson and Mehlig (2016) K. Gustavsson and B. Mehlig, “Statistical models for spatial patterns of heavy particles in turbulence,” Adv. Phys. 65, 1 (2016).
  * Klett and Davis (1973) J. D. Klett and M.H. Davis, “Theoretical collision efficiencies of cloud droplets at small Reynolds numbers,” J. Atmos. Sci. 30, 107–117 (1973).
  * Candelier (2007) Fabien Candelier, “Relative motion of two identical bubbles, rising in a fluid at rest, at small Reynolds numbers,” Comptes Rendus Mecanique 335, 732–736 (2007).
  * Rother  _et al._ (2022) M.A. Rother, J.K. Stark, and R.H. Davis, “Gravitational collision efficiencies of small viscous drops at finite Stokes numbers and low Reynolds numbers,” International Journal of Multiphase Flow 146, 103876 (2022).
  * Harrison  _et al._ (2015) R. G. Harrison, K. A. Nicoll, and M. H. P. Ambaum, “On the microphysical effects of observed cloud edge charging,” Q. J. R. Meteorol. Soc. 141, 2690–2699 (2015).
  * Davis (1965) M. H. Davis, “The effect of electric charges and fields on the collision of very small cloud drops,” Proceedings of the International Conference on Cloud Physics, Tokyo and Sapporo. (1965).
  * Tinsley and Zhou (2014) B.A. Tinsley and L. Zhou, “Comments on ”Effect of electric charge on collisions between cloud droplets”,” Journal of Applied Meteorology and Climatology 53, 1317–1320 (2014).
  * Paluch (1970) I. Paluch, “Theoretical collision efficiencies of charged cloud droplets,” Journal of Geophysical Research (1896-1977) 75, 1633–1640 (1970).
  * Schlamp  _et al._ (1976) R. J. Schlamp, S. N. Grover, H. R. Pruppacher, and A. E. Hamielec, “Numerical Investigation of the Effect of Electric Charges and Vertical External Electric Fields on the Collision Efficiency of Cloud Drops.” Journal of the Atmospheric Sciences 33, 1747–1755 (1976).
  * Schlamp  _et al._ (1979) R. J. Schlamp, S. N. Grover, H. R. Pruppacher, and A. E. Hamielec, “A Numerical Investigation of the Effect of Electric Charges and Vertical External Electric Fields an the Collision Efficiency of Cloud Drops: Part II,” Journal of the Atmospheric Sciences 36, 339–349 (1979).
  * Abbott (1974) C. E. Abbott, “Experimental cloud droplet collection efficiencies,” Journal of Geophysical Research 79, 3098–3100 (1974).
  * Magnusson  _et al._ (2022) G. Magnusson, A. Dubey, R. Kearney, G. P. Bewley, and B. Mehlig, “Collisions of micron-sized charged water droplets in still air,” Phys. Rev. Fluids 7, 043601 (2022).
  * Lekner (2012) J. Lekner, “Electrostatics of two charged conducting spheres,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 468, 2829–2848 (2012).
  * Khain  _et al._ (2004) A. Khain, V. Arkhipov, M. Pinsky, Y. Feldman, and Ya Ryabov, “Rain enhancement and fog elimination by seeding with charged droplets. Part i: Theory and numerical simulations,” Journal of Applied Meteorology 43, 1513–1529 (2004).
  * Fletcher (2013) N. H. Fletcher, “Effect of Electric Charge on Collisions between Cloud Droplets,” Journal of Applied Meteorology and Climatology 52, 517–520 (2013).
  * Hidy and Brock (1965) G. M. Hidy and J. R. Brock, “Some remarks about the coagulation of aerosol particles by Brownian motion,” J. Coll. Science 20, 477–491 (1965).
  * Lu  _et al._ (2010a) J. Lu, H. Nordsiek, E.W. Saw, and R. A. Shaw, “Clustering of charged inertial particles in turbulence,” Phys. Rev. Lett. 104, 184505 (2010a).
  * Lu  _et al._ (2010b) J. Lu, H. Nordsiek, and R. A. Shaw, “Clustering of settling charged particles in turbulence: theory and experiments,” New J. Phys. 12, 123030 (2010b).
  * Lu and Shaw (2015) J. Lu and R. A. Shaw, “Charged particle dynamics in turbulence: Theory and direct numerical simulations,” Phys. Fluids 27, 065111 (2015).
  * Jackson (1975) J. D. Jackson,  _Classical electrodynamics; 2nd ed._ (Wiley, New York, NY, 1975).
  * Jeffrey and Onishi (1984) D. J. Jeffrey and Y. Onishi, “Calculation of the resistance and mobility functions for two unequal rigid spheres in low-Reynolds-number flow,” Journal of Fluid Mechanics 139, 261–290 (1984).
  * Townsend (2018) A. K. Townsend, “Generating, from scratch, the near-field asymptotic forms of scalar resistance functions for two unequal rigid spheres in low-Reynolds-number flow,” arxiv:1802.08226 (2018).
  * Li Sing How  _et al._ (2021) M. Li Sing How, D. L. Koch, and L. R. Collins, “Non-continuum tangential lubrication gas flow between two spheres,” Journal of Fluid Mechanics 920, A2 (2021).
  * Li Sing How (2021) M. Li Sing How,  _Numerical simulation of near-contact motion and coalescence of inertial droplets in turbulence_ , Ph.D. thesis, Cornell (2021).
  * Bhatnagar  _et al._ (1954) P. L. Bhatnagar, E. P. Gross, and M. Krook, “A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems,” Phys. Rev. 94, 511–525 (1954).
  * Wiggins (2013) S. Wiggins,  _Chaotic Transport in Dynamical Systems_, Interdisciplinary Applied Mathematics (Springer New York, 2013).
  * Zhang and Davis (1991) X. Zhang and R. H. Davis, “The rate of collisions due to Brownian or gravitational motion of small drops,” Journal of Fluid Mechanics 230, 479–504 (1991).
  * Ho and Higuchi (1968) N. F. H. Ho and W. I. Higuchi, “Preferential Aggregation and Coalescence in Heterodispersed Systems,” Journal of Pharmaceutical Sciences 57, 436–442 (1968).
  * Candelier and Mehlig (2016) F. Candelier and B. Mehlig, “Settling of an asymmetric dumbbell in a quiescent fluid,” J. Fluid Mech. 802, 174–185 (2016).
  * Qian and Law (1997) J. Qian and C. K. Law, “Regimes of coalescence and separation in droplet collision,” J. Fluid Mech. 331, 59–80 (1997).
  * Haus and Melcher (1989) H. A. Haus and J. R. Melcher,  _Electromagnetic Fields and Energy_ (Prentice-Hall, 1989).
  * Rumble (2022) J. Rumble,  _CRC Handbook of Chemistry and Physics, 103rd Edition_ (Taylor & Francis, 2022).
  * Candelier  _et al._ (2019) F. Candelier, B. Mehlig, and J. Magnaudet, “Time-dependent lift and drag on a rigid body in a viscous steady linear flow,” Journal of Fluid Mechanics 864, 554–595 (2019).
  * Takahashi (1973) T. Takahashi, “Measurement of electric charge of cloud droplets, drizzle, and raindrops,” Reviews of Geophysics 11, 903–924 (1973).
  * Alipchenkov  _et al._ (2004) V. M. Alipchenkov, L. I. Zaichik, and O. F. Petrov, “Clustering of charged particles in isotropic turbulence,” High Temperature 42, 919–927 (2004).