1. Introduction
Weather modification activities such as rain enhancement and fog elimination are essential for combating water crises in desert areas and increasing visibility on roads and airport runways. To facilitate these activities, one needs to accelerate the cloud microphysical processes, particularly the collision–coalescence process. Cloud seeding (i.e. seeding of hygroscopic agents in clouds) is one artificial method of enhancing the efficiency of droplet–droplet collisions by increasing the number density of large droplets in the size spectrum. The evolution of the drop size distribution (DSD) in warm clouds depends primarily on the rate of collisions between two drops. This collision–coalescence process sets the time to rain formation by shaping the DSD. Thus, a detailed study of the collision process can aid in analysing the growth mechanism of cloud droplets and improving the parameterizations of cloud microphysical processes. Better parameterizations, in turn, can significantly reduce the uncertainties involved in weather forecasting and climate models. In aerosol impactors, the knowledge of the collision process between charged particles can reduce errors in measuring the aerosol size distribution. The droplet volume fraction in a cloud is approximately 
$O(10^{-6})$ (see Grabowski & Wang Reference Grabowski and Wang2013). So, we consider that the system is dilute; thus, the pairwise interactions will serve our purpose. We assume the two spheres to be rigid, ignoring the role of interfacial mobility. This assumption is also valid for smaller water droplets (
${\rm radii} < 25\,\mathrm {\mu }$m) in air where the drop-to-medium viscosity ratio is high (
${\approx }10^2$). Here, we study the role of electrostatic forces in the collision rate of a pair of like-charged tiny droplets sedimenting due to gravity in quiescent air.
The rate equation for the particle number density when two species are present is
\begin{equation} -\frac{{\rm d} n_1}{{\rm d} t} ={-}\frac{{\rm d} n_2}{{\rm d} t} = K_{12}, \end{equation}
where 
$K_{12}$ is the collision rate of particles with radii 
$a_1$ and 
$a_2$ and respective number densities 
$n_1$ and 
$n_2$. The collision rate is difficult to predict theoretically, especially when the effects of interparticle interactions are not inconsequential. Smoluchowski (Reference Smoluchowski1917) derived the expression for the ideal collision rate 
$K_{12}^0$ for two non-interacting spheres settling under gravity in a still fluid and found that 
$K_{12}^0=n_1 n_2 [2{\rm \pi} (\rho _p-\rho _f)(a_1^2-a_2^2)g(a_1+a_2)^2]/(9\mu _f)$, where 
$\rho _p$ and 
$\rho _f$ are the particle and fluid densities, 
$g$ is the acceleration due to gravity and 
$\mu _f$ is the dynamic viscosity of the fluid. The interparticle interactions significantly alter the relative velocity between a particle pair at close separations, thus modifying the collision rate. The collision efficiency 
$E_{12}=K_{12}/K_{12}^0$, which is the ratio of the collision rate with interactions to that obtained ignoring interactions (the ideal collision rate), captures the effects of interparticle interactions on the collision rate. Davis (Reference Davis1984) and Melik & Fogler (Reference Melik and Fogler1984) used trajectory analyses for predicting the collision efficiency of two unequal-sized rigid spheres sedimenting due to gravity and interacting via continuum hydrodynamics and van der Waals forces. Zhang & Davis (Reference Zhang and Davis1991) and Rother, Stark & Davis (Reference Rother, Stark and Davis2022) performed similar calculations for differentially sedimenting viscous drops, where they quantified the effects of interfacial mobility on the collision rate.
We assume both the fluid and particle inertia are negligible, and thus inertia does not influence the collision dynamics. The Reynolds number 
$Re_p=[2\rho _f(\rho _p-\rho _f)g a_1^3]/(9\mu _f^2)$, capturing the fluid inertia, is defined here based on the terminal settling speed and radius of the larger particle. The inertial effects of two differentially settling spheres can be quantified by the Stokes number 
$St=[16\rho _p(\rho _p-\rho _f)g(a_1^2-a_2^2)(a_1a_2)^{3/2}]/[81\mu _f^2(a_1+a_2)^2]$ (see Davis Reference Davis1984). This Stokes number depends on particle sizes and the size ratio 
$\kappa =a_2/a_1 < 1$. The Péclet number (
$Pe$) measuring the relative importance of gravitational sedimentation to Brownian diffusion is given by 
$Pe=2{\rm \pi} (\rho _p-\rho _f)a_1^4 \kappa (1-\kappa ^2)g/(3k_B T)$, where 
$k_B=1.318 \times 10^{-23}$ J K
$^{-1}$ is the Boltzmann constant, and 
$T$ is the absolute temperature (see Zinchenko & Davis Reference Zinchenko and Davis1994). Let us calculate the values of 
$Re_p$, 
$St$ and 
$Pe$ for a water droplet in air with 
$a_1 = 10\,\mathrm {\mu }$m, 
$\rho _p \approx 10^3$ kg m
$^{-3}$, 
$\rho _f \approx 1$ kg m
$^{-3}$, 
$\mu _f \approx 1.8 \times 10^{-5}$ Pa s and 
$T=275$ K. We find 
$Re_p \approx 0.007$ (negligible fluid inertia); 
$St \approx 0.54$ for 
$\kappa = 0.3$ and 
$St \approx 0.03$ for 
$\kappa =0.99$; 
$Pe \approx 1504$ for 
$\kappa = 0.3$ and 
$Pe \approx 1085$ for 
$\kappa =0.99$. The above representative values of 
$St$ suggest that the particle inertia is negligible for nearly equal-sized droplets of radii less than 
$15\,\mathrm {\mu }$m. In warm clouds, condensation is the dominant growth mechanism for droplets of radii up to 
$15\,\mathrm {\mu }$m, thus leading to a nearly monodisperse size distribution. Therefore, the negligible particle inertia assumption is valid at the lower end of the ‘size gap’ of 15–40 
$\mathrm {\mu }$m droplets. The role of particle inertia in the collision rates of a particle pair interacting via continuum hydrodynamics has been studied by Davis (Reference Davis1984) (for rigid spheres) and Rother et al. (Reference Rother, Stark and Davis2022) (for spherical drops). Our non-inertial calculation will work as a reference calculation against which to compare future studies on collisions of charged spheres with inertial effects. For Brownian-dominated collisions 
$Pe \ll 1$ and gravity-dominated collisions 
$Pe \gg 1$. As the 
$Pe$ values are sufficiently large for droplets of radii more than 
$10\,\mathrm {\mu }$m, we carry out the current analysis for the pure gravitational settling case (i.e. Brownian diffusion is negligible). Wen & Batchelor (Reference Wen and Batchelor1985) solved an advection–diffusion equation using an asymptotic method to predict the collision efficiency in the 
$Pe \gg 1$ limit. Zinchenko & Davis (Reference Zinchenko and Davis1994) performed the collision calculations at arbitrary 
$Pe$ for droplets with continuum hydrodynamic interactions and van der Waals forces.
Interparticle interactions, especially hydrodynamic interactions, modulate the collision rates between particle pairs. At close separations, the continuum assumption of hydrodynamic interactions would not be valid, and the near-field non-continuum lubrication interactions become the dominant collision mechanism in media with long mean free paths, like air (see Sundararajakumar & Koch Reference Sundararajakumar and Koch1996). The Knudsen number 
$Kn=\lambda _0/a^*$, where 
$\lambda _0$ is the mean free path of the medium and 
$a^*=(a_1+a_2)/2$ is the average of the sphere radii 
$a_1$ and 
$a_2$, measures the strength of non-continuum effects. Previous studies have obtained collision rates due to non-continuum interactions for particles subject to differential sedimentation, uniaxial compressional flow (Dhanasekaran, Roy & Koch Reference Dhanasekaran, Roy and Koch2021a), simple shear flow (Patra, Koch & Roy Reference Patra, Koch and Roy2022), Brownian motion (Patra & Roy Reference Patra and Roy2022) and turbulent flow (Dhanasekaran, Roy & Koch Reference Dhanasekaran, Roy and Koch2021b). Surface deformations of droplet pairs in Stokes flow become significant when the lubrication pressure becomes comparable to the Laplace pressure. Droplet deformation becomes important when 
$3\mu _f V_{{rel}}a^*/2h^{*2} \sim 2\sigma /a^*$, where 
$V_{{rel}}=2 (\rho _p-\rho _f) (a^2_1-a^2_2) g/(9 \mu _f)$ is the relative velocity of two unequal-sized non-interacting droplets settling under gravity, 
$h^*$ is the gap thickness (in 
$\mathrm {\mu }$m) between the surfaces and 
$\sigma$ is the surface tension at the air–water interface. Gopinath & Koch (Reference Gopinath and Koch2002) showed that deformation becomes important when 
$h^*$ is approximately equal to 
$6.74 \times 10^{-5} a^{*2}$, where the average droplet radius is in 
$\mathrm {\mu }$m. However, we ignore the surface deformations of small water droplets moving in air. This is a reasonable assumption because van der Waals attraction and non-continuum effects will occur before deformation (see Dhanasekaran et al. Reference Dhanasekaran, Roy and Koch2021a). However, it is true that drops ultimately must deform to coalesce. Therefore, drops interacting through non-continuum hydrodynamics must begin to deform at small separations where the attraction force due to van der Waals or electrostatic interactions is large and drops have reached the point of no return for coalescence. For collision rate calculations, we are concerned about whether drops will come up to the separation, after which they must coalesce. So, we do not need to resolve the drop deformation to calculate the collision rate. In the present work, we will analyse the collision rate of non-deformable charged droplets sedimenting through a quiescent atmosphere while interacting through non-continuum hydrodynamics and electrostatic forces.
Cloud droplets can acquire electric charge through different mechanisms, such as diffusion of ions, convection charging, inductive charging, thermoelectric effects and contact potential effects (see Pruppacher & Klett Reference Pruppacher and Klett1997). Tinsley et al. (Reference Tinsley, Rohrbaugh, Hei and Beard2000) studied the effect of image charge forces on the collision between a charged aerosol particle (a point charge) and a charged conducting droplet and concluded that electrostatic forces result in a significant increase in collision efficiency as compared with that for an uncharged aerosol and droplet. They also reported that electrical effects accelerate the scavenging rate of charged aerosol particles even in non-thunderstorm clouds. Khain et al. (Reference Khain, Arkhipov, Pinsky, Feldman and Ryabov2004) obtained an approximate expression for the electrostatic force between two charged conducting spheres using the method of electrical images and showed that the electrostatic force between the charged droplets makes the collision process much more efficient. They suggested that droplet seeding in clouds can be an efficient tool for rain enhancement and fog elimination. However, none of these studies have considered appropriate electrostatic and hydrodynamic interactions between the droplets. Also, all studies have assumed the charged water droplets as conductors, which is not entirely correct because the dielectric constant of water is not infinite. In the present study, we will explore the role of finite dielectric constant on the collision dynamics of inertialess charged droplets sedimenting in still air.
Collision rate calculations of charged drops require information on the sizes and electrical charges of the interacting drop pairs. However, the size of a drop and the amount of charge it carries are not independent parameters. There are several field measurements on droplet sizes and charges in a weakly electrified cloud (see Twomey Reference Twomey1956; Krasnogorskaya Reference Krasnogorskaya1969; Colgate & Romero Reference Colgate and Romero1970). Almost all field observations reported a quadratic relationship between the droplet charge 
$q_1$ and droplet radius 
$a_1$ (i.e. 
$q_1 \propto a_1^2$). Various empirical fits of the measured droplet sizes and charges are available in the literature (see Pruppacher & Klett Reference Pruppacher and Klett1997). Rayleigh (Reference Rayleigh1882) derived the maximum charge carrying capacity of a drop by equating the repulsive electrostatic stress and the stabilizing surface tension stress and found 
$q_{{max}}$ to be 
$8{\rm \pi} (\epsilon _0\sigma _0a_1^3)^{1/2}$, where 
$\epsilon _0\approx 8.85 \times 10^{-12}$ C
$^2$ N
$^{-1}$ m
$^{-2}$ is the permittivity of free space or air, 
$\sigma _0 = 0.073$ N m
$^{-1}$ is the surface tension of water. If the amount of charge exceeds the Rayleigh limit charge, the drop will experience electromechanical instability and disintegrate into several smaller drops. Corona discharge occurs when the drop surface potential exceeds a breakdown value 
$V_b$. For drops with 
$a_1 < 100\,\mathrm {\mu }$m, 
$V_b \approx 0.327$ kV. Based on this condition, the maximum charge 
$q_{{max}}$ is determined to be 
$4{\rm \pi} \epsilon _0a_1V_b$, referred to as the Paschen limit. Figure 1(a) shows the variation of droplet charge with droplet radius for four different formulas, and it is evident from the figure that the charge of a cloud drop is much less than the Rayleigh or Paschen limits. In § 2, we will define a non-dimensional parameter 
$N_e$ measuring the relative strength of the electrostatic and gravitational forces. In figure 1(b), we show how 
$N_e$ varies with the droplet radius while charges on droplets vary according to the relations given in figure 1(a).
(a) The droplet charge as a function of radius from four different formulas. The top two lines show the maximum droplet charge based on the Rayleigh and Paschen limit. The blue line indicates the charge measured by Twomey (Reference Twomey1956), and the black line is an empirical fit given in Pruppacher & Klett (Reference Pruppacher and Klett1997). (b) The quantity 
$N_e$ as a function of droplet radius for charge values given in panel (a) when the size ratio is 
$0.5$.
The calculation of electrostatic interactions between two charged spheres has a long history. Maxwell (Reference Maxwell1873) calculated the electrostatic energy of two spheres of different sizes and charges, obtaining the expression for the interaction energy as an expansion in inverse powers of 
$r$ (centre-to-centre distance) until 
$O(r^{-22})$. Maxwell's calculation obtains the self-energy and Coulomb energy at leading order. At the next order, the mutual polarization term appears, and it is always negative irrespective of the sign of the charges. Russell (Reference Russell1909) derived the expressions for capacitance coefficients of two equal-sized charged conducting spheres at small separations. Davis (Reference Davis1964) calculated the electrostatic forces between two charged conductors in an imposed electric field by solving the potential field in a bispherical coordinate system. Following the work of Russell (Reference Russell1909), Lekner (Reference Lekner2012, Reference Lekner2016) provided an analytical expression for the electrostatic force near contact for arbitrary values of size ratio and charge ratio. He found that two like-charged conducting spheres always attract each other at short distances, except for those charge ratio values that the spheres would attain when they are brought into contact. The surface potentials of the spheres are equal in this exceptional case. Thus we would almost always encounter attractive near-field interactions when studying particle collisions in a polydisperse suspension of conducting spheres. Calculating the interaction force between charged dielectric spheres is more challenging than their perfectly conducting counterparts since it requires the additional calculation of the non-uniform electric potential inside the spheres. Feng (Reference Feng2000) used the Galerkin finite difference method to estimate the electrostatic interaction force between two charged dielectric spheres in contact. The theoretical calculation of the interaction forces between two charged dielectric spheres has received considerable attention recently (see Bichoutskaia et al. Reference Bichoutskaia, Boatwright, Khachatourian and Stace2010; Munirov & Filippov Reference Munirov and Filippov2013; Khachatourian et al. Reference Khachatourian, Chan, Stace and Bichoutskaia2014). Two like-charged dielectric spheres, like conducting spheres, almost always attract each other at close separations. However, the charge ratio–size ratio parameter space in which they repel each other forms a band region (in the case of conducting spheres, this repulsive region becomes a curve; see figure 4). The attraction between like-charged spheres occurs because of surface charge redistribution due to mutual polarization. In the current analysis, we utilize the work of Khachatourian et al. (Reference Khachatourian, Chan, Stace and Bichoutskaia2014), who calculated the electrostatic interaction force between two charged spheres of dielectric materials by using the general solution of the Laplace equation in bispherical coordinates.
We will study the collision efficiency of charged spheres of dielectric material settling in still air. In § 2, we will formulate the problem and outline the procedure for calculating the collision efficiency. Then, in § 3, we will calculate the collision efficiency of a pair of hydrodynamically interacting spheres due to the combined effect of gravity, electrostatic and van der Waals forces. Finally, in § 4, we will summarize our results and discuss their implications.
2. Problem formulation
2.1. Expression for the relative velocity of two particles
We consider a dilute polydisperse suspension of charged spherical particles settling under gravity in still air and interacting with each other through non-continuum hydrodynamics. In dilute systems, the probability of a third particle affecting the relative motion of two interacting particles is negligible. Thus we can restrict our analysis to binary interactions of particles with radii 
$a_1$ and 
$a_2$ as shown in figure 2(a). Particles are assumed inertialess due to their small sizes, and fluid motion is assumed to be sufficiently slow, satisfying the Stokes equations for creeping flow. Since the equations governing the flow field are linear, we can write the resultant relative velocity between a pair of spheres as a vector sum of the relative velocities caused by gravity, van der Waals and electrostatic forces (see Batchelor Reference Batchelor1982; Davis Reference Davis1984)
\begin{align} \boldsymbol{\hat{v}}_{12}(\boldsymbol{\hat{r}}) &= \boldsymbol{\hat{v}}^{(0)}_{12}\boldsymbol{\cdot} \left[L\frac{\boldsymbol{\hat{r}\hat{r}}}{\hat{r}^2} + M\left(\boldsymbol{\mathsf{I}}-\frac{\boldsymbol{\hat{r}\hat{r}}}{\hat{r}^2}\right)\right] \nonumber\\ &\quad +\frac{1}{6{\rm \pi} \mu_f }\left(\frac{1}{a_1} + \frac{1}{a_2}\right) \left[G\frac{\boldsymbol{\hat{r}\hat{r}}}{\hat{r}^2}+H\left(\boldsymbol{\mathsf{I}}- \frac{\boldsymbol{\hat{r}\hat{r}}}{\hat{r}^2}\right)\right]\boldsymbol{\cdot} \left(\boldsymbol{F}_{vdW}+\boldsymbol{F}_e\right), \end{align}
where 
$\boldsymbol {\hat {r}}$ is the vector from the centre of particle 2 to the centre of particle 1, 
$\hat {r} = |\boldsymbol {\hat {r}}|$, 
$\boldsymbol{\mathsf{I}}$ is the unit second-order tensor, 
$\boldsymbol {\hat {v}}^{(0)}_{12} = 2 (\rho _{p1}-\rho _f) (a^2_1-\gamma a^2_2) \boldsymbol {g}/(9 \mu _f)$ is the relative velocity between the particle pairs in the absence of interaction and 
$\boldsymbol {F}_{vdW}$ and 
$\boldsymbol {F}_e$ are, respectively, the van der Waals and electrostatic forces acting on the particles. Moreover, 
$\rho _{pj}$ is the density of the 
$j$th particle and 
$\gamma =(\rho _{p2}-\rho _f)/(\rho _{p1}-\rho _f)$ captures scenarios when the particles have different densities. The viscous drag on the two particles having a finite size difference can be different even if their masses are the same. Therefore, a pair of sedimenting spheres would have non-zero relative velocity if they were of unequal sizes, densities or both. The motivation for the current study is to analyse the settling dynamics of charged water droplets in air; thus, we assume 
$\rho _{p1}=\rho _{p2}=\rho _p\rightarrow \gamma =1$. Here, 
$L, M$ are axisymmetric mobility (representing the relative motion along the line of centres) functions for two unequal-sized spherical particles settling under gravity through a quiescent fluid, and 
$G, H$ are asymmetric mobility (representing the relative motion normal to the line of centres) functions for two spherical particles interacting hydrodynamically and moving because of central potentials. These mobility functions depend on the size ratio, 
$\kappa = a_2/a_1$, and non-dimensional centre-to-centre distance, 
$r=2\hat {r}/(a_1+a_2)$. For calculating continuum axisymmetric mobilities, we utilize the solution of the Stokes equations in a bispherical coordinate system (see Lin, Lee & Sather Reference Lin, Lee and Sather1970; Wang, Zinchenko & Davis Reference Wang, Zinchenko and Davis1994; Zinchenko & Davis Reference Zinchenko and Davis1994). We calculate continuum asymmetric mobilities using the twin multipole expansion technique developed by Jeffrey & Onishi (Reference Jeffrey and Onishi1984). Recently, Dhanasekaran et al. (Reference Dhanasekaran, Roy and Koch2021a) calculated the modifications of the axisymmetric mobilities due to non-continuum lubrication interactions, where they considered continuum lubrication interactions when the separation 
$\xi = r-2 > O(Kn)$ and non-continuum lubrication interactions when the 
$\xi \leq O(Kn)$. They utilized the work of Sundararajakumar & Koch (Reference Sundararajakumar and Koch1996) that provided the non-continuum lubrication force between particles colliding in a gaseous medium. In the present analysis, we use the uniformly valid solution of axisymmetric mobilities developed by Dhanasekaran et al. (Reference Dhanasekaran, Roy and Koch2021a). We expect that continuum breakdown will not strongly influence asymmetric mobilities because these mobilities remain finite at contact. Thus we consider continuum hydrodynamics for asymmetric motion at all separations. We also expect that the non-continuum interactions for the tangential motions will not significantly alter our collision rate calculations because axisymmetric and asymmetric relative motions of an inertialess particle pair are independent. However, one must consider non-continuum lubrication forces for asymmetric motions for calculating the collision rate of inertial particles in a gas (see Li Sing How, Koch & Collins Reference Li Sing How, Koch and Collins2021).
(a) Schematic of the coordinate system used in the analysis. Here, ‘1’ indicates the sphere with radius 
$a_1$ and charge 
$q_1$; ‘2’ indicates the sphere with radius 
$a_2$ and charge 
$q_2$. The sphere marked ‘3’ is the collision sphere of radius 
$a_1+a_2$. The electrostatic force 
$\boldsymbol {F}_e$ acts along the line joining the centres of the two spheres, and gravity acts along the negative 
$x_3$ direction. We use 
$\hat {e}_r$ and 
$\hat {e}_{\theta }$ as the unit vectors in the 
$r$ and 
$\theta$ directions, respectively. (b) Schematic of three trajectories: the red line is a collision trajectory, the golden line is the limiting trajectory, and the blue line is an open trajectory.
The van der Waals force 
$F_{vdW} = -{\rm d} \varPhi _{vdW}/{\rm d} r$, where 
$\varPhi _{vdW}$ is the van der Waals potential, always acts along the line of centres of the two spheres. Using a pairwise additivity theory, Hamaker (Reference Hamaker1937) derived an analytical expression for 
$\varPhi _{vdW}$. However, Hamaker's calculation ignored the electromagnetic retardation. One must consider the effects of retardation when the separation is comparable to or more than the London wavelength 
$\lambda _L (\approx 0.1\,\mathrm {\mu }{\rm m})$. In the present analysis, we use the work of Zinchenko & Davis (Reference Zinchenko and Davis1994), who derived the expression for the retarded van der Waals potential between two particles. The retarded 
$\varPhi _{vdW}$ is a function of 
$r$, 
$\kappa$, 
$A_H$ and 
$N_l$. Here, 
$A_H$ is the Hamaker constant and the non-dimensional quantity 
$N_l$ is the radius of the spheres scaled with 
$\lambda _L$ (i.e. 
$N_l = 2 {\rm \pi}(a_1+a_2)/\lambda _L = 2 {\rm \pi}a_1 (1+\kappa )/\lambda _L$). The values of 
$A_H$ for several common materials are available in the literature (see Friedlander Reference Friedlander2000).
The electrostatic interaction force between two charged conducting spheres has been studied extensively. Most importantly, the analytical expressions of the force for far and close interparticle separations are now well established (see Lekner Reference Lekner2012). Lekner (Reference Lekner2012) derived the near-field asymptotic expression for the force by approximating the capacitance coefficients for small separations. The study found that a pair of like-charged conducting spheres always attract each other at close separations except when the charge ratio 
$\beta = q_2/q_1$ obeys the following relation (shown as the yellow line in figure 4):
\begin{equation} \beta=\frac{\gamma+\psi\left(\dfrac{a_1}{a_1+a_2}\right)}{\gamma+\psi \left(\dfrac{a_2}{a_1+a_2}\right)}=\frac{\gamma+\psi\left(\dfrac{1}{1+\kappa}\right)}{\gamma+ \psi\left(\dfrac{\kappa}{1+\kappa}\right)}, \end{equation}
where 
$\gamma = 0.5772156649\cdots$ is the Euler constant, 
$\psi$ is the digamma function (i.e. 
$\psi (z)={\rm d}\ln \varGamma (z)/{\rm d} z$ is the logarithmic derivative of the gamma function 
$\varGamma (z)$) and 
$q_1$ and 
$q_2$ are the charges of the spheres with radii 
$a_1$ and 
$a_2$ respectively. We assume that the charges reside entirely on the surface of the particles; there are no volumetric charges.
The near-field results mentioned above rely on the asymptotic expansions for the capacitance coefficients to find the electrostatic energy for small separations and thus the force (Banerjee et al. Reference Banerjee, Peters, Brown and Song2021; Lekner Reference Lekner2021). The capacitance method will work only for the case of perfect conductors. To find the electrostatic force for small gaps for large but finite dielectric constant, one can adopt the lubrication analysis of Batchelor & O'Brien (Reference Batchelor and O'Brien1977) who calculated the effective thermal/electric conductivity of a suspension of densely packed highly conducting spheres. Khair (Reference Khair2013) did a similar lubrication analysis for perfect conductors, confirming the leading-order asymptotic expression of Lekner (Reference Lekner2012). For completeness, we formulate the lubrication analysis to derive the expression for the electrostatic force between two charged spheres separated by a small distance. We choose a cylindrical coordinate system (
$\rho, \varphi, z$) with the origin at the midpoint of the two interacting spheres, as shown in figure 3. The potential 
$\phi$ is invariant to rotation around the 
$z$-axis, and thus, Laplace's equation for 
$\phi$ reduces to
\begin{equation} \frac{\partial^2 \phi}{\partial z^2} + \frac{1}{\rho}\frac{\partial}{\partial \rho} \left(\rho\frac{\partial \phi}{\partial \rho}\right) = 0. \end{equation}
The boundary conditions for (2.3) are: 
$\phi = \phi _1$ at the surface 
$S_1$, 
$\phi = \phi _2$ at the surface 
$S_2$ and 
$\phi \rightarrow 0$ as the surface-to-surface distance approaches infinity. For small 
$\rho$, we approximate the spherical surfaces as paraboloids. Thus, the surface to surface distance between the spheres is given by
\begin{equation} h(\rho) = h_0 + \frac{\rho^2}{2a}, \end{equation}
where 
$h_0$ is half of the minimum gap thickness and 
$a=2a_1a_2/(a_1+a_2)$ is the reduced radius. The length scale across the gap (i.e. in the 
$z$ direction) over which 
$\phi$ varies is 
$h_0$, whereas the length scale along the gap for variations in 
$\phi$ is 
$L_{\rho }=O(\sqrt {ah_0})$ (see Khair Reference Khair2013). So, the latter is 
$O(\sqrt {h_0/a}) (\ll 1)$ smaller than the former, and hence at the leading order, (2.3) becomes 
$\partial ^2\phi /\partial z^2 = 0$. Hence, the variation of potential within the gap is given by
\begin{equation} \phi = \frac{\phi_1+\phi_2}{2} + \frac{\phi_1-\phi_2}{2}\frac{z}{h(\rho)}. \end{equation}
We consider that the particles carry unequal but fixed charges. The amount of charge carried by the particle 
$1$ can be expressed as
\begin{equation} q_1 ={-}\epsilon \oint\frac{\partial \phi}{\partial n}{\rm d} S_1, \end{equation}
with 
${\rm d} S_1=\rho \,{\rm d}\varphi \,{\rm d}\rho$ with 
$\varphi$ being the azimuthal coordinate. For sphere 1, 
$\partial \phi /\partial n=-\partial \phi /\partial z=-(\phi _1-\phi _2)/2h(\rho )$. Thus (2.6) becomes
\begin{equation} q_1 \sim \epsilon {\rm \pi}\int_0^R\frac{(\phi_1-\phi_2)\rho}{h_0 + (\rho^2/2a)}{\rm d}\rho. \end{equation}
The upper limit 
$R (\gg L_{\rho })$ of the above integration corresponds to the outer boundary of the lubrication region. Here, 
$q_2=-q_1$ to the leading order because 
$\partial \phi /\partial z$ is constant across the gap. This simply says that two almost touching particles subject to a fixed potential difference acquire charges of opposite signs but equal magnitudes at the leading order. Thus we can write
\begin{equation} q_1-q_2 \sim 2\epsilon {\rm \pi}\int_0^R\frac{(\phi_1-\phi_2)\rho}{h_0 + (\rho^2/2a)}{\rm d}\rho. \end{equation}
For small gap thickness, the force on each particle acts along the 
$z$ direction. Let us say that the force on particle 
$1$ is 
$\textbf { {F}}_1=F\hat {e}_z$, where 
$\hat {e}_z$ is the unit vector in the 
$z$ direction and
\begin{equation} F = \frac{1}{2}\epsilon \oint\left(\frac{\partial \phi}{\partial n}\right)^2 \boldsymbol{\hat{n}}_1\boldsymbol{\cdot}\boldsymbol{\hat{e}}_z \,{\rm d} S_1. \end{equation}
Here, 
$\boldsymbol {\hat {n}}_1$ is the outward unit normal to 
$S_1$ and at leading order 
$\boldsymbol {\hat {n}}_1\boldsymbol {\cdot }\boldsymbol {\hat {e}}_z=-1$. So, (2.9) becomes
\begin{equation} F \sim{-}\frac{\epsilon {\rm \pi}}{4} \int_0^R\frac{(\phi_1-\phi_2)^2\rho}{\left(h_0 + (\rho^2/2a)\right)^2}{\rm d}\rho. \end{equation}
By symmetry, the force on particle 
$2$ equals 
$-\textbf { {F}}_1$. For perfect conductors, the dielectric constant is infinite, and 
$\phi _1$ and 
$\phi _2$ are constants. The leading-order contribution to the force can be found by taking 
$R \rightarrow \infty$. Thus, the expression for 
$F$ reduces to
\begin{equation} F \sim{-}\frac{\epsilon {\rm \pi}}{4} (\phi_1-\phi_2)^2 \frac{a}{h_0}. \end{equation}
Similarly, we can simplify the expression for the charge difference between particles 
$1$ and 
$2$ as follows:
\begin{align} q_1-q_2 &\sim 2\epsilon {\rm \pi}(\phi_1-\phi_2) \int_0^R \frac{\rho}{h_0 + (\rho^2/2a)}{\rm d}\rho \nonumber\\ &= 2\epsilon {\rm \pi}(\phi_1-\phi_2) \int_0^a\frac{\rho\,{\rm d}\rho}{h_0 + (\rho^2/2a)} + \cdots \nonumber\\ &\sim 2\epsilon {\rm \pi}a(\phi_1-\phi_2)\ln[a/h_0] . \end{align}Thus, the potential difference between the spheres for the small gap thickness can be expressed in terms of the charge difference as
\begin{equation} \phi_1-\phi_2 \sim \frac{q_1-q_2}{2\epsilon {\rm \pi}a \ln(a/h_0)}. \end{equation}
Putting 
$\phi _1-\phi _2$ from (2.13) into (2.11), we have
\begin{equation} F \sim{-}\frac{1}{16 \epsilon {\rm \pi}a^2}\left(\frac{a}{h_0}\right)\frac{(q_1-q_2)^2}{[\ln(a/h_0)]^2}. \end{equation}Schematic of two-sphere geometry at close separation. Here, 
$k$ and 
$\alpha k$ are the relative permittivities of the fluid medium and sphere materials.
Equation (2.14) implies that the attraction force has an 
$O(\xi ^{-1}[\ln \xi ]^{-2})$ singularity, which approaches 
$-\infty$ as 
$\xi \rightarrow 0$. From this, we can conclude that two like-charged spherical conductors interacting through continuum hydrodynamics can come into contact in finite time. In practice, the perfect conductor approximation is sufficiently accurate for particles with very high dielectric constants, such as metal particles. However, this assumption would not be valid for cloud drops since the dielectric constant for water is finite (
$k\approx 80$).
In the present work, we consider the effects of finite dielectric constant on the collision rate. Depending on the size ratio and charge ratio values, a pair of like-charged spheres of dielectric materials attract one another at small surface-to-surface distances (see Bichoutskaia et al. Reference Bichoutskaia, Boatwright, Khachatourian and Stace2010; Munirov & Filippov Reference Munirov and Filippov2013; Khachatourian et al. Reference Khachatourian, Chan, Stace and Bichoutskaia2014). Here, we will utilize the work of Khachatourian et al. (Reference Khachatourian, Chan, Stace and Bichoutskaia2014), who calculated the electrostatic interaction force between two charged spheres of dielectric materials using bispherical coordinates. In Appendix A, we outline the procedure for computing the non-dimensional electrostatic force 
$f_e = |\boldsymbol {F}_e|/(q_1^2/4{\rm \pi} \epsilon _0a_1^2)$ for given values of 
$r$, 
$\kappa$, 
$\beta$, 
$k_1$ (dielectric constant of sphere 1), 
$k_2$ (dielectric constant of sphere 2) and 
$k_m$ (relative permittivity of the fluid medium). Figure 4(a) shows the contour plot of 
$\ln |f_e|$ in the 
$\beta - \kappa$ parameter space for two like-charged water drops in air when they almost touch each other. To illustrate the distinct feature in the case of dielectric, in figure 4(b), we further show the contour plot of 
$f_e$ for the same values of 
$k$, 
$\beta$, 
$k_1$, 
$k_2$ when 
$0 \leq \beta \leq 1$ and 
$0 \leq \kappa \leq 1$. The dark red region in figure 4(b) indicates the positive values of 
$f_e$ (i.e. repulsion). This band-like region becomes a curve given by (2.2) (the yellow line in figure 4b) in the case of perfectly conducting spheres (i.e. 
$k_1=k_2=\infty$). For a given size ratio, the attractive electrostatic force in near-field decreases with increasing the charge ratio, and it turns into a repulsive force for a certain range of 
$\beta$ values. With a further increase in 
$\beta$, the force again becomes attractive and increases monotonically with 
$\beta$. The inset in figure 4(b) shows this behaviour for 
$\kappa =0.5$. The electrostatic force between two like-charged dielectric spheres deviates significantly from that of a perfect conductor case when the dielectric constant is higher than that of water (see figure 5). Therefore, in the lubrication regime, collision dynamics of dielectric spherical pairs is expected to be altered from that of a pair of perfect conductors. The inset in figure 5 shows that the magnitude of the attraction force at close separation increases with increasing dielectric constant.
Variation of electrostatic force with size and charge ratio when the separation 
$\xi =10^{-6}$. Panel (a) shows the contour plot of the logarithm of the normalized electrostatic force (normalized by 
$q_1^2/4{\rm \pi} \epsilon _0a_1^2$) at contact in the charge ratio–size ratio parameter space when 
$k_1=k_2=80$. In (b), the dark red region indicates the repulsive forces (i.e. the force values are positive in this region). Everywhere else, the forces are attractive (i.e. negative). The yellow line is for a conducting spherical pair along which the forces are repulsive. The inset in (b) indicates the non-monotonic variation of the near-field electrostatic force with the charge ratio when the size ratio is 
$0.5$.
The non-dimensional electrostatic forces as a function of separation for different dielectric constant when 
$\kappa = 0.5$, 
$\beta =1.0$. This plot shows that the interaction forces between two dielectric spheres deviate more from the perfect conductor case when they come close to each other. The inset shows how the forces at contact vary with the dielectric constants.
We choose a spherical coordinate system 
$(r,\theta,\varphi )$ with the origin at the centre of sphere 2 with 
$\theta =0$ (i.e. 
$x_3$-axis) being the gravity axis. To non-dimensionalize the system, we consider 
$a^*$, 
$\hat {v}^{(0)}_{12}=|\boldsymbol {\hat {v}}^{(0)}_{12}|$ and 
$a^*/\hat {v}^{(0)}_{12}$ as the characteristic length, velocity and time scale of the problem. Thus the non-dimensional radial separation between the centres of the two spheres lies in the range 
$r=2$ (referred to as the collision sphere, indicated as sphere 3 in figure 2) to 
$\infty$ (where one sphere does not influence the other). We denote the dimensionless coordinates with an overbar, i.e. 
$\bar {x}_i=2x_i/(a_1+a_2)$ (
$i=1,2,3$). The size ratio 
$\kappa$, which can vary in the range 
$(0,1]$, captures the geometry of the two-sphere system. The dimensionless relative velocity 
$\boldsymbol {v} = \boldsymbol {\hat {v}}_{12}/\hat {v}^{(0)}_{12}$ can be written as 
$\boldsymbol {v} = v_r \hat {e}_r + v_{\theta } \hat {e}_{\theta }$, where
\begin{gather} v_r = \frac{{\rm d} r}{{\rm d} t} ={-} L\cos\theta - N_v G\frac{{\rm d} \varPhi_{vdW}}{{\rm d} r} + N_e G f_e, \end{gather}
\begin{gather}v_{\theta} = r\frac{{\rm d} \theta}{{\rm d} t} = M \sin\theta, \end{gather}
and 
$N_v$ and 
$N_e$ are dimensionless quantities that capture relative strengths of van der Waals and electrostatic forces to gravity
\begin{gather} N_v = \frac{3 A_H}{2{\rm \pi}\kappa\left(1-\kappa^2\right)(\rho_p-\rho_f){g} a_1^4}, \end{gather}
\begin{gather}N_e = \frac{3q_1^2}{16{\rm \pi}^2\epsilon_0(\rho_p-\rho_f)g\kappa(1-\kappa)a_1^5}. \end{gather}
In the definition of the non-dimensional parameter 
$N_e$ above, when we refer to gravity force in the ratio of electrostatic to gravity, we consider the force due to the sum of gravity and viscous forces that drives the relative motion of the two drops. The particles are inertialess (
$St=0$); thus, there is an instantaneous balance between viscous drag and external forces.
2.2. Expressions for the particle collision efficiency
We consider a dilute dispersion with different drop sizes and aim to determine the rate at which drops of radii 
$a_1$ and 
$a_2$ with number densities 
$n_1$ and 
$n_2$ collide with each other per unit volume. Mathematically, the collision rate 
$K_{12}$ is equal to the flux of pairs into the collision sphere of radius 
$\hat {r}=a_1+a_2$ and can be expressed in terms of the pair distribution function 
$P(\hat {r})$ and the relative velocity 
$\boldsymbol {\hat {v}}_{12}$ by
\begin{equation} K_{12} ={-}n_1n_2 \int_{(\hat{r}=a_1+a_2)\ {\rm and}\ \left(\boldsymbol{\hat{v}}_{12} \boldsymbol{\cdot}\textbf{{n}}<0\right)} \left(\boldsymbol{\hat{v}}_{12} \boldsymbol{\cdot}\textbf{{n}}\right) P\,{\rm d} A, \end{equation}
where 
$\textbf { {n}}$ is the outward unit normal at the collisional surface. For a dilute dispersion, the pair distribution function is governed by the quasi-steady Fokker–Planck equation for regions of space outside the contact surface
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\left(P\ \boldsymbol{\hat{v}}_{12}\right) = 0. \end{equation}
The absence of far-field correlations defines the upstream boundary condition: 
$P \rightarrow 1$ as 
$\hat {r} \rightarrow \infty$. For the calculation purpose, we take 
$\hat {r}=\hat {r}_{\infty }$, which is large but finite.
The relative motion of a pair of non-Brownian drops settling under gravity in a quiescent fluid is deterministic, and thus we can use a trajectory analysis to find the collision rate. Using (2.20) and the divergence theorem, the integral in (2.19) can be taken over the surface that encloses the volume occupied by all trajectories that originate at 
$\hat {r}=\hat {r}_{\infty }$ and terminate at 
$\hat {r}=a_1+a_2$. The flux through the cross-sectional area 
$A_c$ of this volume at 
$\hat {r}=\hat {r}_{\infty }$ determines the collision rate. As the motion due to gravity is symmetric about the gravity axis, the area 
$A_c$ is a circle. We label 
$A_c$ as the upstream interception area. Since 
$P=1$ and 
$\boldsymbol {\hat {v}}_{12}=\boldsymbol {\hat {v}}^{(0)}_{12}$ at 
$\hat {r}=\hat {r}_{\infty }$, the collision rate is
\begin{equation} K_{12} ={-}n_1n_2\int_{A_c}\left(\boldsymbol{\hat{v}}_{12}\boldsymbol{\cdot} \textbf{{n}}'\right)|_{\hat{r}_{\infty}} P\,{\rm d} A = n_1n_2\hat{v}^{(0)}_{12}{\rm \pi} x_{2c}^2, \end{equation}
where 
$\textbf { {n}}'$ is the unit outward normal vector at the area element of 
$A_c$ and 
$x_{2c}$ is the critical impact parameter (i.e. the radius of the circle at 
$\hat {r}=\hat {r}_{\infty }$). In other words, this critical impact parameter is the horizontal distance from the gravity axis for two widely separated drops that will graze each other. The pair trajectory of two grazing drops is called the limiting trajectory that separates the open and close trajectories (see figure 2b). Equation (2.21) bypasses the evaluation of the pair probability 
$P$. The collision rate 
$K_{12}^{0}$ without interparticle interactions and Brownian diffusion is given by the classical Smoluchowski model, where 
$x_{2c}=a_1+a_2$. Thus the ideal collision rate becomes
\begin{equation} K_{12}^{0} = n_1n_2\hat{v}^{(0)}_{12}{\rm \pi} \left(a_1+a_2\right)^2. \end{equation}
The collision efficiency 
$E_{12}$ is defined as the ratio of 
$K_{12}$ to 
$K_{12}^{0}$
\begin{equation} E_{12} = \frac{K_{12}}{K_{12}^{0}} = \frac{x_{2c}^2}{\left(a_1+a_2\right)^2} = \frac{1}{4}\bar{x}_{2c}^2, \end{equation}
where 
$\bar {x}_{2c}=2x_{2c}/(a_1+a_2)$ is the dimensionless critical impact parameter. For 
$\bar {x}_2<\bar {x}_{2c}$, the drops collide; for 
$\bar {x}_2>\bar {x}_{2c}$, the drops move past one another without colliding. The problem becomes reduced to determining 
$\bar {x}_{2c}$. We find the limiting trajectory by integrating the following dimensionless trajectory equation:
\begin{equation} \frac{{\rm d}\theta}{{\rm d} r}=\frac{v_{\theta}}{r v_r}= \frac{M\sin\theta}{r\left({-}L\cos\theta - N_v G \dfrac{{\rm d} \varPhi_{vdW}}{{\rm d} r} + N_e G f_e\right)}. \end{equation}
The closed-form analytical solution for the collision efficiency exists for 
$N_v=N_e=0$ in (2.24). The expression of 
$E_{12}$ for 
$N_v=N_e=0$ is given by (see Davis Reference Davis1984)
\begin{equation} E_{12}=\exp\left({-}2\int_2^{\infty}\frac{M-L}{r L}{\rm d} r\right). \end{equation}3. Results and discussion
We have discussed in § 1 that the interaction force between two like-charged spheres can be attractive at close separations depending on the size and charge ratio values. However, like-charged particles always repel each other for moderate to large separations. Here, we aim to quantify how electrostatic interactions affect the collision dynamics. Therefore, we initially consider a situation where two like-charged spheres settle under gravity and interact via non-continuum hydrodynamics and electrostatic forces but do not experience van der Waals forces. Later, we will expand the analysis and determine the collision rate for cases where non-continuum hydrodynamics, van der Waals and electrostatic forces act together. We will present all the results for droplets with 
$k_1=k_2=80$.
A collision trajectory is a path traced by an evolving satellite sphere (sphere 1) that starts far upstream and eventually collides with the test sphere (sphere 2) located at the origin. The relative radial velocity must be inward at all separations for a colliding trajectory (see the condition 
$v_r<0$ in (2.19)). Thus, from (2.15), we have 
$N_e<(L\cos \theta )/(G f_e)$ for 
$r \in [2,\infty )$. The radial relative velocity due to gravity (the first term on the right-hand side in (2.15)) is inward for 
$0\leq \theta <{\rm \pi} /2$, and it is maximum at 
$\theta =0$. Therefore, the condition of getting at least one collision trajectory becomes 
$N_e<(L/G)/ f_e$ for 
$r \in [2,\infty )$. For a given 
$\kappa$ and 
$\beta$, we found a critical 
$N_e$ value above which no collision trajectory exists. We denote this critical 
$N_e$ as 
$(N_e)_c$. To obtain 
$(N_e)_c$, we compute the quantity 
$(L/G)/ f_e$ as a function of separation and 
$(N_e)_c$ equal to the lowest positive value of 
$(L/G)/ f_e$. For a given 
$\kappa$ and 
$Kn$, 
$(N_e)_c$ decreases monotonically as the charge ratio 
$\beta$ increases (see figure 6a). In figure 6(b), we demonstrate the importance of 
$(N_e)_c$ by plotting 
$v_r$ with 
$\xi$ when 
$\theta =0$. As expected, 
$v_r$ is negative for all 
$\xi$ when 
$N_e<(N_e)_c$, and thus a pair of drops starting far apart will eventually collide in this condition. For 
$N_e=(N_e)_c$, 
$v_r$ is zero at one particular 
$\xi$ and negative for all other values of 
$\xi$. Beyond the critical 
$N_e$, there are some separation distances where the magnitude of the inward radial velocity due to gravity is less than the radial velocity due to electrostatic repulsion; hence, no collision will occur. Primarily gravity governs the relative motion at large separations since the electrostatic force decreases rapidly with increasing 
$r$.
(a) The computed values of the critical 
$N_e$ as a function of charge ratio 
$\beta$ for different 
$\kappa$ when 
$Kn=10^{-2}$ and 
$N_v=0$. (b) The radial component of the relative velocity 
$v_r$ as a function of separation for 
$N_e$ less than, equal to and greater than 
$(N_e)_c$ when 
$Kn=10^{-2}$, 
$\kappa =0.5$, 
$\beta =10$ and 
$N_v=0$.
Electrostatic interaction forces will modify the relative trajectories of the pure differential sedimentation problem. To obtain relative drop trajectories, we set initial conditions on the collision sphere and perform backward integrations of (2.24) using a fourth-order Runge–Kutta method (the ‘ode45’ subroutine in MATLAB). Exactly at 
$r = 2$, 
$v_r = 0$ since the hydrodynamic mobilities 
$L$ and 
$G$ are zero at 
$r = 2$. However, in the case of non-continuum lubrication interactions, 
$v_r$ approaches zero slowly enough with decreasing 
$\xi$ that the pair would come into contact in finite time. To avoid the locations of zero relative radial velocity, we start the integrations with initial conditions on a sphere of radius 
$2+\zeta$, where 
$\zeta$ is a small separation from the contact surface. We will show converged results without too much computational load when 
$\zeta = 10^{-6}$. Because of symmetry about the gravity axis, it is sufficient to analyse the problem in 
$\theta \in [0,{\rm \pi} ]$.
As mentioned earlier, the non-continuum lubrication force allows collision to occur in finite time. For 
$N_e=0$, we found that two types of pair trajectories (open and colliding) exist (see figure 7a) and the separatrix (the limiting colliding trajectory) that separates these two types of trajectories touches the contact surface at 
$\theta ={\rm \pi} /2$. Therefore, the trajectories that depart from the contact surface with an initial 
$\theta$ in the range 
$[0,{\rm \pi} /2]$ and go to infinity are collision trajectories. It is noticeable from figure 7(a) that open trajectories are top–bottom symmetric for 
$N_e=0$. With increasing 
$Kn$, the non-dimensional critical impact parameter and hence the collision efficiency will increase. For 
$0< N_e<(N_e)_c$, the types of trajectories remain the same (i.e. open and colliding). However, unlike the 
$N_e=0$ case, the satellite sphere can collide on the rear side of the test sphere. Depending on the value of 
$N_e$, the 
$\theta$ location of the limiting trajectory on the contact surface varies between 
${\rm \pi} /2$ and 
${\rm \pi}$. We take a sufficiently large number of initial points on the collision sphere with 
$\theta$ in the range 
$[{\rm \pi} /2,{\rm \pi} ]$ and then determine the relative trajectories for all initial conditions by backward integrating (2.24). Out of these trajectories, the limiting colliding trajectory is the one whose distance from the axis of symmetry is maximum in the far field. We find that, as 
$N_e$ increases from 
$0$ to 
$(N_e)_c$, the initial 
$\theta$ locations of the limiting trajectories on the collision sphere gradually shift from 
${\rm \pi} /2$ to 
${\rm \pi}$. Figure 7(b) shows typical pair trajectories when 
$\kappa =0.5$, 
$\beta =10$, 
$Kn=10^{-2}$ and 
$N_e=10^{-1}$. For these parameter values, 
$(N_e)_c$ is approximately 
$1.77$. It is evident from figure 7(b) that the electrostatic interaction breaks the top–bottom symmetry of open trajectories. For 
$N_e > (N_e)_c$, gravity fails to overcome the electrostatic repulsion force for some relative separations, and thus no collision trajectory exists in this case. Figure 7(c) shows the pair trajectories for the same values of 
$\kappa$, 
$\beta$ and 
$Kn$ when 
$N_e=1.8 (>(N_e)_c)$. As expected, in this case, a thin layer appears at some distance from the collision sphere where 
$v_r$ is positive. The green line in figure 7(c) indicates this layer. A trajectory starting from far upstream will never be able to cross this layer, and therefore, two surfaces will not come into contact. An external non-hydrodynamic attractive force or particle inertia can overcome this repulsive region and bring particles up to a separation where electrostatic forces are attractive. Therefore, colliding trajectories can exist in those cases. We do not consider those scenarios in the present study. However, to prove our argument, we consider a hypothetical situation where trajectories starting with an initial 
$r$ less than the radial coordinate of this layer can eventually hit the contact sphere (see pink lines in figure 7c). As per the definition in (2.19), these trajectories do not contribute to the collision rate of settling drops with initially large separations.
Relative trajectories of two like-charged sedimenting spheres (gravity acts vertically downwards, i.e. along negative 
$x_3$-axis) for 
$Kn=10^{-2}$, 
$\kappa =0.5$, 
$\beta =10$ and 
$N_v=0$ when (a) 
$N_e=0$, (b) 
$N_e=0.1$ and (c) 
$N_e=1.8$. The value of 
$(N_e)_c$ for the above specified values of 
$Kn$, 
$\kappa$ and 
$\beta$ is approximately 
$1.77$. Collision trajectories exist for the first two cases. For 
$N_e > (N_e)_c$, collision is not possible.
Hydrodynamic interactions retard the relative velocity at small separations, thus making the collision process less efficient. The collision efficiency, defined as the ratio of the collision rate with interparticle interactions to the ideal collision rate, quantifies the effect of interactions. In this case, the collision efficiency depends on 
$Kn$ (measures the strength of non-continuum hydrodynamic interactions), 
$\kappa$ (describes the relative geometry of interacting spheres), 
$\beta$ (represents the relative electrostatic strength of the spheres) and 
$N_e$ (captures the electrostatic interactions). We span these parameters to obtain the important features of the collision dynamics. We present most of our results for the Knudsen number 
$Kn=10^{-2}$, a realistic value for water droplets in warm cumulus clouds. For a given 
$Kn$, 
$\kappa$ and 
$\beta$, we vary 
$N_e$ in the range 
$(0,(N_e)_c)$. Figure 8 shows the variation of the collision efficiency with 
$N_e$ for 
$Kn=10^{-2}$ and 
$\kappa =\beta =0.3,0.5$ when 
$N_v=0$. For 
$N_e$ greater than or equal to 
$(N_e)_c$, the collision efficiency 
$E_{12}$ becomes zero. As we decrease 
$N_e$ from 
$(N_e)_c$, 
$E_{12}$ increases and attains a maximum before it decreases again. With further decrease of 
$N_e$, 
$E_{12}$ asymptotes to the value corresponding to the collision efficiency due to non-continuum lubrication interactions of two unequal-sized uncharged spheres settling under gravity through a still fluid medium. The variation of 
$E_{12}$ with 
$N_e$ is qualitatively similar for the two values of 
$\kappa$. The 
$\kappa =\beta$ line lies in the positive band of the force contour (given in figure 4b) beyond 
$\kappa =\beta \approx 0.8$, and thus no collision will occur for 
$\kappa =\beta > 0.8$.
Collision efficiencies as a function of 
$N_e$ when 
$Kn=10^{-2}$, 
$\kappa = 0.3, 0.5$ and 
$\beta = \kappa$. The dashed lines indicate collision efficiencies due to non-continuum lubrication interactions when 
$N_v=0$.
We find that the electrostatic interactions can either impede or assist the collision process of a pair of like-charged droplets settling under gravity, and the collision dynamics strongly depends on the parameters 
$\kappa$, 
$\beta$ and 
$N_e$. Nonetheless, the collision efficiency will asymptote to the results corresponding to 
$N_e \rightarrow 0$ (the effects of gravity dominate) and 
$N_e \rightarrow \infty$ (electrostatic forces dominate) limits. In 
$N_e \rightarrow 0$ limit, 
$E_{12}$ will approach gravitational collision efficiency due to non-continuum hydrodynamics alone when 
$N_v=0$ and non-continuum hydrodynamics plus van der Waals interactions when 
$N_v \neq 0$. For 
$N_e \rightarrow \infty$, 
$E_{12}$ will approach zero. Figure 9 shows how collision efficiencies vary with 
$N_e$ when 
$\kappa = 0.3, 0.5, 0.7, 0.9$, 
$Kn = 10^{-2}$, 
$\beta = 0.1, 1, 10$ and 
$N_v = 10^{-3}$. The electrostatic attraction force at small interparticle separations increases with increasing 
$\kappa$ when the charge ratio is small (for example, 
$\beta = 0.1$). With increasing 
$\beta$, this trend reverses gradually and becomes the opposite when 
$\beta$ is high (for example, 
$\beta = 10$). Therefore, for small values of 
$\beta$, the collision efficiency is higher for higher 
$\kappa$ when 
$N_e$ is large. When 
$\beta$ is high, the variation of the collision efficiency with 
$\kappa$ is exactly opposite to that of the small 
$\beta$ cases. On the other hand, when gravity with non-continuum hydrodynamics plus van der Waals forces drive the dynamics (i.e. 
$N_e \rightarrow 0$ regime), the collision efficiency always increases with increasing 
$\kappa$. Since for small charge ratio values, the variations of collision efficiencies with size ratios are similar in the two extreme regimes (i.e. 
$N_e \rightarrow 0$ and 
$N_e \rightarrow \infty$), the curves for different 
$\kappa$ in figure 9(a) do not intersect with each other. The strength of electrostatic attraction forces is weak for small size ratio values (for example, 
$\kappa = 0.3, 0.5$), and thus the collision efficiencies, in these cases, are always smaller than the respective 
$N_e \rightarrow 0$ asymptotic values. However, for high size ratio values (for example, 
$\kappa = 0.7, 0.9$), the electrostatic attraction forces are so strong that the collision efficiencies for a range of values of 
$N_e$ overshoot the results corresponding to 
$N_e = 0$. For 
$\beta = 1.0$, the behaviours of 
$E_{12}$ with 
$\kappa$ are opposite in the two extreme regimes, and thus the curves corresponding to various 
$\kappa$ values intersect with each other in some intermediate 
$N_e$ values (see figure 9b). Also, in this case, collision efficiencies overshoot the 
$N_e=0$ results for lower 
$\kappa$ values because the magnitudes of electrostatic attraction forces are higher for smaller 
$\kappa$. For 
$\beta = 10$, the electrostatic attraction forces in the lubrication regime are so strong that the collision efficiencies overshoot the results corresponding to the non-continuum plus van der Waals case for all 
$\kappa$ values (see figure 9c). However, as in the previous case, different 
$\kappa$ curves intersect before approaching the respective asymptotes. One must note that the non-monotonic behaviour of the collision efficiency for different 
$\beta$ and a fixed 
$\kappa$ is consistent with our prior observation of the near-field variation of 
$f_e$ with 
$\beta$ (see the inset of figure 4b).
Collision efficiencies as a function of 
$N_e$ for 
$Kn=10^{-2}$, 
$\kappa = 0.3, 0.5, 0.7, 0.9$, 
$N_v = 10^{-3}$ when (a) 
$\beta = 0.1$, (b) 
$\beta = 1$ and (c) 
$\beta = 10$.
So far, we have analysed the role of size and charge ratios on the collision rate of two like-charged spheres settling in a non-continuum gas. Now, we will quantify how the strength of non-continuum effects and van der Waals attraction force influence the collision rate, and to do so, we will fix 
$\kappa$ to 
$0.5$ and 
$N_e$ to 
$1.0$. Figure 10 shows the variation of the collision efficiency with 
$Kn$ when 
$N_v=10^{-3}$ and 
$\beta =0.1, 1.0, 10$. Taking the effects of non-continuum lubrication interactions into account, Dhanasekaran et al. (Reference Dhanasekaran, Roy and Koch2021a) calculated the collision efficiency of two uncharged sedimenting drops, both with and without van der Waals force. The black line in figure 10 obtained using (2.25) represents the collision efficiency due to non-continuum effects alone. As 
$Kn$ decreases, 
$E_{12}$ due to non-continuum effects only decreases monotonically and will approach zero for 
$Kn \rightarrow 0$. The collision rate due to non-continuum hydrodynamics and van der Waals forces (NC
$+$vdW) decreases with decreasing 
$Kn$. Still, it approaches a non-zero value because van der Waals interactions between the particles drive the entire collision dynamics when 
$Kn$ approaches zero (see the red line in figure 10). Collision dynamics becomes more complicated when non-continuum effects, van der Waals and electrostatic interactions (NC
$+$vdW
$+$EF) act simultaneously. For small charge ratio values (for example, 
$\beta =0.1, 1.0$), collision efficiency slowly decreases with decreasing 
$Kn$, and asymptotes to the collision efficiency value that can be obtained by considering full continuum hydrodynamic interactions. Collision efficiencies for small 
$\beta$ values are less than that of NC
$+$vdW result. Interestingly, in this case, the collision efficiency is less than the pure non-continuum result when 
$Kn > 10^{-2}$ and 
$\beta$ is unity. Electrostatic forces for large charge ratio values (for example, 
$\beta =10$) dominate over all other effects, and thus collision efficiency becomes higher than that of the NC
$+$vdW result. Also, 
$E_{12}$ becomes almost independent of 
$Kn$ when 
$\beta$ is large.
Collision efficiencies as a function of 
$Kn$ for 
$\kappa = 0.5$, 
$N_e = 1.0$, 
$N_v = 10^{-3}$ when 
$\beta = 0.1, 1.0, 10$. The black and red lines indicate the collision efficiency for non-continuum (NC) effects alone and non-continuum plus van der Waals (NC
$+$vdW) interactions, respectively. Blue lines represent the collision efficiencies for different charge ratios when non-continuum hydrodynamics, van der Waals and electrostatic forces act together (NC
$+$vdW
$+$EF).
In figure 11, we show the variation of collision efficiencies with 
$N_v$ for 
$\kappa =0.5$, 
$Kn=10^{-2}$, 
$N_e=1.0$ when 
$\beta =0.1, 1.0, 10$. For 
$\beta =0.1$, 
$E_{12}$ decreases monotonically as the strength of the van der Waals force decreases. It asymptotes to the value corresponding to the collision efficiency with non-continuum effects plus electrostatic interactions (NC+EF). The collision efficiency for higher charge ratio values (
$\beta =1.0, 10$) is independent of the van der Waals force because electrostatic forces dominate the collision dynamics. As expected, collision efficiency for 
$\beta =10$ is much higher than that of 
$\beta =0.1$ or 
$10$. To show how electrostatic interactions between like-charged droplets influence the collision dynamics, we compare our findings with the collision efficiency of a pair of uncharged drops settling under gravity and interacting through the continuum and non-continuum hydrodynamic interactions and van der Waals forces (see the black and red lines in figure 11).
Collision efficiencies as a function of 
$N_v$ for 
$\kappa = 0.5$, 
$Kn=10^{-2}$, 
$N_e = 1.0$ when 
$\beta = 0.1, 1.0, 10$. The black and red lines indicate the collision efficiency with full continuum hydrodynamics plus van der Waals force and NC+vdW, respectively. Blue lines represent the collision efficiencies with NC+vdW+EF for different charge ratios.
Now, let us estimate the collision efficiency of a physical system with 
$a_1=10\,\mathrm {\mu }$m and 
$q_1=200e$ (
$e=1.602 \times 10^{-19}$ C is the elementary charge) as a function of 
$\kappa$ when the charge ratio 
$\beta = 5, 10$ (see figure 12). The mean free path of air 
$\lambda _0$ increases with altitude (Wallace & Hobbs Reference Wallace and Hobbs2006) and is approximately 
$0.1\,\mathrm {\mu }$m for warm cumulus clouds. Thus the Knudsen numbers as a function of size ratio are given by 
$Kn = 0.02/(1+\kappa )$. For water droplets in air, 
$A_H \approx 3.7 \times 10^{-20}$ J (see Friedlander Reference Friedlander2000). So, the 
$N_l$ and 
$N_v$ written in terms of 
$\kappa$ are given by 
$N_l = 200{\rm \pi} (1+\kappa )$ and 
$N_v = 1.8 \times 10^{-4}/\kappa (1-\kappa ^2)$. The dimensionless parameter 
$N_e$ varies with 
$\kappa$ according to the relation 
$N_e=8.58\times 10^{-4}/\kappa (1-\kappa )$. Given the numerical values of 
$\kappa$, we calculate the relevant dimensionless parameters from the above relations and then calculate the collision efficiencies. Collision efficiencies initially decrease with increasing 
$\kappa$ and then increase when 
$\kappa$ becomes close to unity. For a fixed 
$\kappa$, the collision efficiency is larger for 
$\beta = 10$ than for 
$\beta = 5$. As discussed in § 1, gravity-induced collisions without electrostatic interactions have extensive treatment in the literature (see Davis Reference Davis1984; Zhang & Davis Reference Zhang and Davis1991; Dhanasekaran et al. Reference Dhanasekaran, Roy and Koch2021a). We show these results in figure 12 to demonstrate that electrostatic interactions can enhance collision efficiency for higher charge ratios.
Collision efficiencies as a function of 
$\kappa$ for 
$a_1=10\,\mathrm {\mu }$m, 
$q_1=200 e$ and 
$\beta = 5, 10$. We considered the combined effects of non-continuum hydrodynamics (NC), van der Waals interactions (vdW) and electrostatic forces (EF) in performing these collision calculations. The results from previous studies for uncharged drops are included for reference.
4. Summary and conclusions
We have quantified the effects of electrostatic interactions on the collision dynamics of tiny charged particles sedimenting through a quiescent gaseous medium. We restricted our analysis to dilute dispersions, where pairwise interactions between the particles are dominant, and we can safely neglect the interactions between three or more particles. When the surface-to-surface distance between two interacting spheres is less than the mean free path of the gas, the lubrication interactions become non-continuum, resulting in a finite collision rate even in the absence of attractive potentials. We utilized the work of Dhanasekaran et al. (Reference Dhanasekaran, Roy and Koch2021a), who calculated the uniformly valid axisymmetric mobility functions that capture non-continuum lubrication interactions at close separations and full continuum hydrodynamic interactions at moderate to large separations. The charged water droplets in atmospheric clouds might not behave like perfect conductors, mainly when the sphere separation is less than 
$0.01 a_1$ (see Pruppacher & Klett Reference Pruppacher and Klett1997). Therefore, to accurately capture the dynamics, we considered the electrostatic force between two spheres of dielectric materials. We utilized the work of Khachatourian et al. (Reference Khachatourian, Chan, Stace and Bichoutskaia2014), who computed the electrostatic interaction force between two like-charged dielectric spheres. Using trajectory analysis, we calculated the collision efficiencies of two differentially sedimenting spheres interacting through the continuum and non-continuum hydrodynamics, van der Waals and electrostatic forces for a wide range of size ratio and charge ratio values. Finally, we compared our findings with the existing studies on gravity-induced collisions of uncharged spheres.
In § 2, we formulated the problem for inertialess systems. Relevant mobility functions that accurately represent the hydrodynamic and interparticle interactions have been used to express the relative velocity between two inertialess spheres due to the combined effects of gravity, van der Waals attractions and electrostatic forces. To obtain the non-dimensional relative velocity equations, we used the mean radius of the two spheres, the relative velocity without interactions and the mean radius divided by the relative velocity without interactions as the characteristic length, velocity and time scales. The parameters 
$N_v$ and 
$N_e$ in the non-dimensional pair trajectory equation capture the relative strength of van der Waals and electrostatic forces to gravity. With the relative velocity and the pair probability at contact as the integrand, we have expressed the non-dimensional collision rate equation as a surface integral over the collision sphere. We have determined the collision efficiency by calculating the critical impact parameter using trajectory analysis.
We found that there exists critical 
$N_e$ for given values of the Knudsen number, charge ratio and size ratio. No collision occurs for 
$N_e > (N_e)_c$. The electrostatic interaction force alters the pair trajectories and the collision efficiencies due to pure differential sedimentation. Figure 7 provides a typical map of pair trajectories in the absence of van der Waals forces for three different cases: 
$N_e = 0$, 
$0 < N_e \leq (N_e)_c$ and 
$N_e > (N_e)_c$. We have presented the typical variation of the collision efficiency with non-continuum lubrication interactions as a function of 
$N_e$ when 
$N_v = 0$ (see figure 8). The interplay between different forces made the dynamics even more complicated when gravity, non-continuum hydrodynamics, van der Waals and electrostatic forces act together. We have shown that polydispersity in terms of particle sizes and charges has a significant influence on the collision dynamics. As 
$N_e$ increases from zero, we observed two different types of dynamics depending on the size and charge ratio values. In one case, the collision efficiency increases to a maximum and then decreases to zero. In the other case, the collision efficiency decreases monotonically with increasing 
$N_e$ and asymptotes to zero as 
$N_e \rightarrow (N_e)_c$.
It is evident from this study that collisions of like-charged droplets can be much more efficient than collisions of uncharged ones. Thus the seeding of charged aerosol particles (i.e. the cloud seeding technology) that charge the droplets can increase the raindrop formation rate in clouds by accelerating the collision process. Nowadays, cloud seeding is an artificial method of rain enhancement in countries suffering from water shortages. Cloud seeding is also an efficient tool that can improve visibility and ensure safe operations of roads, airport runways and industrial facilities by clearing fog and low clouds. Better knowledge of the collision process of charged droplets in natural clouds might improve the DSD evolution prediction. The current work might also help in designing systems where differential sedimentation drives collisions between charged particles interacting through non-continuum hydrodynamics and van der Waals interactions. This work ignores the role of Brownian diffusion on the collision dynamics. Recently, Patra & Roy (Reference Patra and Roy2022) calculated the Brownian motion-induced (
$Pe=0$) collision rate with non-continuum hydrodynamics, van der Waals and electrostatic interactions. In many physical scenarios, both Brownian diffusion and gravitational settling would act in tandem during collisions of micron-sized particles. For the coupled problem, one needs to solve the advection–diffusion equation for the pair probability at arbitrary 
$Pe$. Zinchenko & Davis (Reference Zinchenko and Davis1994) calculated the collision efficiency at arbitrary 
$Pe$ in the continuum limit where van der Waals forces and/or interfacial mobility induce collisions. Thus, extending the current study to include the effects of Brownian diffusion would be logical. We plan to communicate this work in the future.
In the current study, we have assumed that the charges reside entirely on the surface of the particles; there are no volumetric charges. The particles are assumed to be perfect dielectrics. However, a cloud droplet is an aqueous solution, and ions would migrate towards the charged interface. Thus a question arises – would the surface charge change during the duration of droplet collisions? To answer this question, we will compare the drop interaction time with the other relevant time scales associated with the evolution of the surface charge. The time scale of drop encounter due to differential sedimentation would be 
$t_{{coll}} = a^*/V_{{rel}}$, 
$t_{{coll}} \approx 0.04$ s for water droplets in air with 
$a_1=10\,\mathrm {\mu }$m and 
$\kappa =0.99$. In an aqueous solution, there is substantial conductivity in the form of the electrical migration of ions. The diffusion time scale for the ions would be 
$t_{{diff}}=a_1^2/D$. For a typical diffusivity of ions 
$D=10^{-10}$ m
$^2$ s
$^{-1}$ for a 
$10\,\mathrm {\mu }$m drop, 
$t_{{diff}}=1$ s. Thus 
$t_{{coll}} \ll t_{{diff}}$. Another mechanism for altering the surface charge would be charge relaxation in the surrounding air medium due to conduction; the mobility of free ions in the suspending fluid medium can make the medium conductive. Thus, even when there is no corona discharge, the droplet charge can decrease due to the conductivity of the medium. In this case, the instantaneous droplet charge would be 
$q_0\exp (-t/\tau )$, where 
$q_0$ is the initial charge of the droplet, and 
$\tau =\epsilon _0/\sigma$ is the relaxation time with 
$\sigma$ being the conductivity of the fluid (Pruppacher & Klett Reference Pruppacher and Klett1997). The conductivity of clear air at sea level is approximately 
$2.3 \times 10^{-14}$ Sm
$^{-1}$, which gives 
$\tau \approx 6.5$ min. However, the conductivity inside warm clouds is substantially lower than the fair weather sea level conductivity because the concentration of gas-phase free ions in clouds is significantly lower than in cloud-free air. The conductivity inside a cloud can range from 
$1/40$ to 
$1/3$ of the fair weather sea level conductivity, which suggests that 
$\tau$ can vary approximately from 
$20$ min to 
$2$ h (Pruppacher & Klett Reference Pruppacher and Klett1997; Tinsley et al. Reference Tinsley, Rohrbaugh, Hei and Beard2000). Thus 
$t_{{coll}} \ll \tau$, implying charge relaxation to be a much slower process than drop encounter. A possible dominant mechanism for the rearrangement of surface charge would be due to an electrokinetic flow inside charged drop 1 under the influence of the electric field from charged drop 2. When the two drops are sufficiently far apart, this slip velocity would be 
$u_{{slip}} \sim (1/\mu _w)\beta /a_1^3q_1^2/(16{\rm \pi} ^2\epsilon _0)$. This slip velocity would rearrange the surface charge on a convective time scale 
$t_{{rearr}} = 2{\rm \pi} a_1/u_{{slip}} \approx 32 a_1^4 {\rm \pi}^3 \epsilon _0 \mu _w /q_1^2/ \beta$. Thus 
$t_{{coll}}/t_{{rearr}} \sim (\mu _a/\mu _w)\beta N_e$. Here, 
$\mu _w$ and 
$\mu _a$ are the dynamic viscosities of water and air, respectively. The present study shows that collisions for sedimenting charged spheres are most efficient when 
$N_e=O(1)$; thus, a drop encounter event would happen on a shorter time scale than the rearrangement of the surface charge – 
$t_{{coll}} \ll t_{{rearr}}$. Future work might consider the interaction of falling charged water droplets, accounting for the interior electrokinetic flow.
Direct numerical simulations of homogeneous isotropic turbulence seeded with electrically charged, inertial point particles reveal that the electrostatic charge can significantly modify the particle radial relative velocity, radial distribution function and hence the collision rate (see Lu et al. Reference Lu, Nordsiek, Saw and Shaw2010; Lu & Shaw Reference Lu and Shaw2015). These studies assumed pairwise interactions, which is justifiable for dilute suspensions. However, one must account for the influence of all other particles while tracking a single particle in suspensions with high volume fractions. Yao & Capecelatro (Reference Yao and Capecelatro2018) have recently developed an efficient numerical scheme in an Eulerian–Lagrangian framework for accurately capturing the effects of Coulomb forces between charged inertial particles in an isotropic turbulent flow. We want to highlight that any point particle simulation of charged particles will always lead to repulsive interactions between like-charged particles. Our present analysis borrows from classical and more recent electrostatic studies to show that like-charged particles can attract each other near contact, thus fundamentally altering inferences regarding collision and evolution of particle size distributions. Collision rate calculations of like-charged, finite-size particles in a turbulent flow with the appropriate near-field form of the electrostatic forces would be an important subject for future studies. A theoretical study of Zhang, Basaran & Wham (Reference Zhang, Basaran and Wham1995) showed that an external electric field could significantly enhance the coalescence rate of uncharged droplets settling under gravity. Experiments in a cloud chamber by Wang et al. (Reference Wang, Li, Li, Zhang, Yang and Yu2020) demonstrated that a strong electric field can dominate the collision process. Recently, Thiruvenkadam et al. (Reference Thiruvenkadam, Patra, Kadaba and Roy2023) theoretically analysed the effects of an external electric on relative trajectories of a pair of uncharged conducting spheres, and they showed that electric-field-induced forces allow surface-to-surface contact in a finite time by overcoming continuum lubrication resistance. We plan to study the combined effects of gravity, non-continuum interactions, van der Waals forces and electric-field-induced forces on the collision dynamics of uncharged micron-sized droplets.
Funding
The authors acknowledge support from IIT Madras for its support of the ‘Laboratory for Atmospheric and Climate Sciences’ research initiative under the Institute of Eminence framework. P.P. would like to acknowledge financial support from the Prime Minister's Research Fellows (PMRF) scheme, Ministry of Education, Government of India (project no. SB22230184AMPMRF008746). D.L.K. acknowledges support from NSF grant 2206851.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Computation of the electrostatic force
The expression for the electrostatic force is given by
\begin{equation} \boldsymbol{F}_e ={-}\hat{e}_r\frac{1}{K}\sum_{n=0}^{\infty}f_n \left[\frac{n}{2}\varPhi_{1,n-1}\,{\rm e}^{\eta_1}-\left(n+\frac{1}{2}\right) \varPhi_{1,n}+\left(\frac{n+1}{2}\right)\varPhi_{1,n+1}\,{\rm e}^{-\eta_1}\right]\varPhi_{2,n}, \end{equation}
where 
$K=1/(4{\rm \pi} \epsilon _0)\approx 9 \times 10^9$ VmC
$^{-1}$ and 
$\varPhi _{1,n}$ and 
$\varPhi _{2,n}$ are the coefficients of the electrostatic potential generated by the charges residing on spheres 1 and 2, respectively. These coefficients are calculated from two recursive relations. One of them is given by
\begin{align} & \left[\left(n+\frac{1}{2}\right)\varPhi_{1,n}\cosh\eta_1-\frac{n}{2} \varPhi_{1,n-1}-\frac{n+1}{2}\varPhi_{1,n+1}\right]\left(k_m+k_1\right)+ \frac{\sinh\eta_1}{2}\left(k_m-k_1\right)\varPhi_{1,n} \nonumber\\ &\quad -\left[\left(n+\frac{1}{2}\right)\varPhi_{2,n}f_n\cosh\eta_1-\frac{n}{2} \varPhi_{2,n-1}f_{n-1}-\frac{n+1}{2}\varPhi_{2,n+1}f_{n+1}\right]\left(k_m-k_1\right) \nonumber\\ &\quad +f_n\frac{\sinh\eta_1}{2}\left(k_m-k_1\right)\varPhi_{2,n} = \frac{\sqrt{2}a K q_1 \exp\left({-\left(n+\dfrac{1}{2}\right)\eta_1}\right)}{a_1^2}, \end{align}
and the other one is obtained by interchanging subscripts 1 and 2 in (A2). The expressions for 
$\eta _1$, 
$\eta _2$, 
$a$ and 
$f_n$ involved in the above two recursive relations are given by
\begin{gather} \eta_1 = \cosh^{{-}1}\left[\frac{\left(1+\kappa\right)r}{4}+ \frac{\left(1-\kappa^2\right)}{\left(1+\kappa\right)r}\right], \end{gather}
\begin{gather}\eta_2 = \cosh^{{-}1}\left[\frac{\left(1+\kappa\right)r}{4\kappa}+ \frac{\left(1-\kappa^2\right)}{\kappa\left(1+\kappa\right)r}\right], \end{gather}
\begin{gather}a = a_1\sqrt{\frac{\left(1+\kappa\right)^2r^2}{16} + \frac{\left(1-\kappa^2\right)^2}{\left(1+\kappa\right)^2r^2}-\frac{\left(1+\kappa^2\right)}{2}}, \end{gather}
\begin{gather}f_n = \exp\left({-\left(n+\frac{1}{2}\right)\left(\eta_1+\eta_2\right)}\right). \end{gather}The detailed derivations for the electrostatic force and the recursive relations of potential coefficients are given in Khachatourian et al. (Reference Khachatourian, Chan, Stace and Bichoutskaia2014).
Open access
