Relative trajectories of two charged sedimenting particles in a Stokes flow

Motion of particles in a Stokes flow, such as sedimentation, [1–3] has applications including medical technology [4–6], microfluidics [7–11], swimming of microorganisms [12–14], deformation of vesicles [15], waste water treatment [16], marine snow [17], mantle plums [18] and motion within volcanic magma [19, 20]. There are many introductions to the physics of Stokes flows, introduced in [21], such as [22–30] and the included references.

There has recently been interest in bound states of particles in viscous flows. For example, formation of doublets and other bound states of particles in viscous flows have been explored for drops [31], pairs of magnetically active rollers near a repelling wall [32] and for pairs of identical rigid spheres in a background flow with walls [33]. There are also results about large scale spontaneous self-organization into ordered structures with many drops [34] or many rollers [35].

There is also a rich literature in bound states of sedimenting particles. A fundamental problem for the motion of a group of particles close to each other is understanding whether they stay together for a long time or disperse and what are the physical mechanisms responsible for keeping them together, see [36–45]. To make the notion of capturing precise, we define a capturing state of two particles as a configuration whose relative trajectories do not go to infinity in future time. The capturing set is then the set of all capturing states.

In [46] it is shown that there can exist a capturing set of a finite non-zero measure for two uncharged spherical particles of different radii and different masses settling under gravity in a Stokes flow. This set consists of neutrally stable periodic orbits. Because of the neutral stability, there are no basins of attraction, so that even a small perturbation can have a destabilizing effect. Further, because the capturing set has finite measure, if particles begin very distantly they come closer to each other but later move away, and are not trapped. This structure of the capturing set is often used to predict that trajectories of two captured particles are likely to be disturbed by the presence of other particles in a suspension, and therefore capture will probably have no significant effect on the dynamics of sedimenting suspensions e.g., [47–50].

In this paper we study pairs of charged sedimenting particles, in order to see if the charge can create large basins of attraction for bound states. Our interest in electrostatic forces originates from the simple observation that systems of charged particles settling in viscous fluids are common, and therefore stable doublets of charged particles could be potentially used in many practical applications. In a vacuum, electrostatic interactions are destabilizing. This fact is known as Earnshaw's Theorem. However, in [51] it was shown that even a very small charge can stabilize pairs of sedimenting particles; two charged spherical particles of different radii and different masses settling under gravity in a Stokes flow can have locally asymptotically stable stationary states with particle centers in line with gravity.

The stability result in [51] is local. To understand the practical significance of this local result, in this paper we will investigate how two charged point particles sedimenting in a Stokes flow can capture one another, depending on the characteristic parameters of the system: the ratios of the particle radii and masses. We will determine the regions in the phase space of these parameters where stable stationary states of a given interparticle distance can exist. We will determine the capturing set and its structure. We will also find a new class of stable stationary states with particle centers inclined with respect to gravity. We will show that the capturing set of a stable stationary state for two charged point particles settling in a Stokes flow is infinite in measure and consists of one or several basins of attraction to a stable stationary state. These findings open the path toward future experimental observation of these stable doublets. Such doublets could have relevance in dilute charged suspensions [52], for particles in flows at non-zero Reynolds number [53–55] and in plasma [56].

We start our investigation with the mathematical model of point-like particles in section 2. We examine the benchmark case of uncharged pairs of sedimenting point particles (as in [46]) in section 3. In section 4 we present the stability conditions by looking at the linearized dynamics near three kinds of stationary states: centers of particles aligned with gravity with larger particle up, centers of particles aligned with gravity with larger particle down and stationary states where centers of particles are inclined with respect to gravity. In section 5 we discuss generic examples of the relative motion. The vector field of the particle relative velocities and relative trajectories determined numerically on a grid of initial conditions for the case when the larger particle has a greater Stokes velocity. We find the capturing set in physical space of the relative positions using the Poincare-Bendixson theorem. The other case—when the smaller particle has a greater Stokes velocity—is examined in section 6. In section 7 we give the phase diagram for non-overlapping stable stationary states in terms of the ratio of Stokes velocities and ratio of particle radii. We conclude with a summary of our results.

We investigate the dynamics of pairs of charged spherical particles of different radii and masses settling under gravity in a viscous fluid, in the range of the Reynolds number much smaller than unity. A schematic of the system is shown in figure 1. Our goal is to construct the simplest possible analytical model. Therefore, we assume that the particles are point-like and use Coulomb force to model their electrostatic interactions. We also use the point-like model to describe the hydrodynamic interactions. It means that the two-particle mutual mobility is given by the Oseen tensor [26], taken at the relative position of the particles, and self-mobility follows from the Stokes law. The hydrodynamic point-like model is valid when the particles are distant. Further, the electrostatic point-like Coulomb force does not describe effects of, e.g., anisotropic distribution of mobile charges on the particle surface, nor electrostatic screening. Therefore, our approach is an approximation.

Zoom In Zoom Out Reset image size

The particle pair is labeled so that radius a₁ of particle 1 is less than or equal to the radius a₂ of particle 2. Via the Stokes equations, we will arrive at equations of motion in terms of the electrostatic and gravitational forces. These are then non-dimensionalized, parameterized and expressed in terms of coordinates.

2.1. Force diagram and dynamics in the point-force approximation

2.1.1. External forces

We start our diagram of forces with gravity. Let $\hat{{\bf{z}}}$ be a unit vector pointing anti-parallel to the constant gravitational field ${\boldsymbol{g}}$. Further, let m_i be the reduced mass of particle i = 1, 2. This means that the particle weight is corrected for the buoyancy force of the fluid. If the fluid has density ρ and the particle has radius a_i and mass M_i , ${m}_{i}={M}_{i}-\tfrac{4}{3}\pi {a}_{i}^{3}\rho$. The gravitational force on particle i is

$$\begin{eqnarray}&&{{\bf{f}}}_{m,i}=-{m}_{i}g\hat{{\bf{z}}}\end{eqnarray} \tag{ 1 }$$

In this paper, we will assume that m_i  > 0. The case of m_i  ≤ 0 can be studied in an analogous way as seen in the supplemental materials of [51]. Now we move on to the electrostatic field. Let ${{\bf{r}}}_{1}$ and ${{\bf{r}}}_{2}$ be the positions of the centers of particle 1 and 2, so that the relative position is

$$\begin{eqnarray}&&{\bf{d}}={{\bf{r}}}_{2}-{{\bf{r}}}_{1}\end{eqnarray} \tag{ 2 }$$

We now move on to the electrostatic forces. We denote the charge on particle i by q_i . Then the Coulomb electrostatic force which acts on particle i is

$$\begin{eqnarray}&&{{\bf{f}}}_{q,i}={\left(-1\right)}^{i}{{kq}}_{i}{q}_{j}\displaystyle \frac{{\bf{d}}}{| {\bf{d}}{| }^{3}}\end{eqnarray} \tag{ 3 }$$

where k is Coulomb's constant and $j=3-i$. The sum of the external forces acting on particle i is

$$\begin{eqnarray}&&{{\bf{f}}}_{i}={{\bf{f}}}_{m,i}+{{\bf{f}}}_{q,i}.\end{eqnarray} \tag{ 4 }$$

2.1.2. Fluid forces and point force dynamics

We now move on to discuss the fluid forces. We assume that the particles are suspended in an infinite fluid with viscosity μ. We assume that Brownian motion, fluid compressibility and inertia are irrelevant and we describe the fluid flow by the Stokes equations.

$$\begin{eqnarray}&&\mu {{\rm{\nabla }}}^{2}{\bf{u}}-{\rm{\nabla }}p=0\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\bf{u}}=0\,\end{eqnarray} \tag{ 6 }$$

where p is the pressure and ${\bf{u}}$ is the velocity field of a fluid [22, 23]. Physically, modeling fluid interactions by equations (5) & (6) entails that the external forces on the particles in a Stokes fluid are in balance with the resistance forces exerted on the particles by the fluid so that inertia may be neglected.

In the point particle approximation, the hydrodynamic interactions between the particles are described by a linear dependence of the particle velocities ${\dot{{\bf{r}}}}_{i}$ on the external forces ${{\bf{f}}}_{j}$ acting on them. The mutual interaction is determined by the Oseen tensor [26]

$$\begin{eqnarray}&&{\bf{G}}({\bf{d}})=\displaystyle \frac{1}{| {\bf{d}}| }({\bf{I}}+\displaystyle \frac{{\bf{d}}\otimes {\bf{d}}}{| {\bf{d}}{| }^{2}})\end{eqnarray} \tag{ 7 }$$

where $| {\bf{d}}|$ is the length of the vector ${\bf{d}}$, ${\bf{I}}$ is the identity tensor and ⨂ is the tensor product. We assume the self-interaction is determined by the Stokes velocity for the particle given by the external forces the particle is experiencing. To sum up, our dynamic equation will be of the form

$$\begin{eqnarray}&&{\dot{{\bf{r}}}}_{i}=\displaystyle \displaystyle \displaystyle \displaystyle \displaystyle \frac{1}{8{\pi }{\mu }}{\bf{G}}({\bf{d}})\cdot {{\bf{f}}}_{j}+\displaystyle \displaystyle \displaystyle \displaystyle \displaystyle \frac{1}{6{\pi }{\mu }{a}_{i}}{{\bf{f}}}_{i}\end{eqnarray} \tag{ 8 }$$

where $j=3-i$. Notice the above dynamics are first order, which means that the velocities of the particles are given by their positions. This is the mathematical expression of the irrelevance of inertia.

2.2. Characteristic dimensions

In order to aid our analysis we will choose characteristic dimensions. In particular, we choose

$$\begin{eqnarray}&&L={a}_{1}+{a}_{2},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&V=\displaystyle \frac{{m}_{2}g}{6\pi \mu L}\end{eqnarray} \tag{ 10 }$$

as our characteristic length and velocity scale. The time scale is therefore T = L/V. We will use the non-dimensional separation vector

$$\begin{eqnarray}&&{\boldsymbol{\alpha }}=\displaystyle \frac{{\boldsymbol{d}}}{L}\end{eqnarray} \tag{ 11 }$$

and the time normalized by T to nondimensionalize (8). We further define three independent non-dimensional numbers

$$\begin{eqnarray}&&{\beta }=-\,\displaystyle \frac{{{kq}}_{1}{q}_{2}}{{L}^{2}{m}_{2}g}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\gamma =\displaystyle \frac{{a}_{1}}{{a}_{2}},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\delta =\displaystyle \frac{{m}_{1}}{{m}_{2}}\end{eqnarray} \tag{ 14 }$$

which parameterize our equation of motion. One can see that β is the ratio of characteristic electrostatic force to characteristic gravitational force, γ is the ratio of particle radii and δ is the ratio of reduced particle masses. Our assumption that m_i  > 0 entails that δ > 0 and our labeling convention entails γ ≤ 1. An important physical parameter is also the ratio δ/γ of the particle Stokes velocities.

2.3. Dimensionless equations of motion and their basic properties

Having set up our diagram of forces and nondimensionalization choices, starting from the point-force model (8) we arrive at the following non-dimensional dynamical equation [51]

$$\begin{eqnarray}&&\dot{{\boldsymbol{\alpha }}}=\displaystyle \frac{3}{2{\alpha }^{3}}\beta \,{ \mathcal G }\cdot {\boldsymbol{\alpha }}-\beta \displaystyle \frac{{\left(1+\gamma \right)}^{2}}{\gamma {\alpha }^{3}}{\boldsymbol{\alpha }}+\displaystyle \frac{3}{4}(1-\delta ){ \mathcal G }\cdot \hat{{\bf{z}}}-\displaystyle \frac{(\gamma -\delta )(1+\gamma )}{\gamma }\hat{{\bf{z}}}\end{eqnarray} \tag{ 15 }$$

where dot denotes the non-dimensional time derivative, $\alpha =| {\boldsymbol{\alpha }}|$ and ${ \mathcal G }$ is the non-dimensional Green tensor given by

$$\begin{eqnarray}&&{ \mathcal G }({\boldsymbol{\alpha }})=\displaystyle \frac{1}{{\alpha }}\left({\bf{I}}+\displaystyle \frac{{\boldsymbol{\alpha }}\otimes {\boldsymbol{\alpha }}}{{{\alpha }}^{2}}\right)\end{eqnarray} \tag{ 16 }$$

It is helpful to write equation (15) in terms of coordinates of ${\boldsymbol{\alpha }}$. To give a convenient geometry to work in, we will choose x and y axes of the coordinate system so that the particle centers are in the plane y = 0 (the direction of gravity has been already chosen along the z axis). Any orbit with an initial condition in this plane will never experience a force pointing out of this plane, so we will suppress the y-coordinate for the rest of the paper. Using these conventions, equation (15) transformed into Cartesian co-ordinates becomes

$$\begin{eqnarray}&&{\dot{\alpha }}_{x}=3\beta \displaystyle \frac{{\alpha }_{x}}{{\alpha }^{4}}-\beta \displaystyle \frac{{\left(1+\gamma \right)}^{2}}{\gamma }\displaystyle \frac{{\alpha }_{x}}{{\alpha }^{3}}+\displaystyle \frac{3(1-\delta )}{4}\displaystyle \frac{{\alpha }_{x}{\alpha }_{z}}{{\alpha }^{3}}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\dot{\alpha }}_{z}=3\beta \displaystyle \frac{{\alpha }_{z}}{{\alpha }^{4}}-\beta \displaystyle \frac{{\left(1+\gamma \right)}^{2}}{\gamma }\displaystyle \frac{{\alpha }_{z}}{{\alpha }^{3}}+\displaystyle \frac{3(1-\delta )}{4}\displaystyle \frac{{\alpha }_{x}^{2}+2{\alpha }_{z}^{2}}{{\alpha }^{3}}-\displaystyle \frac{(\gamma -\delta )(1+\gamma )}{\gamma }\end{eqnarray} \tag{ 18 }$$

A few words can be said about what writing the dynamics in these coordinates entails. First of all, (17) entails that if the particles begin with their centers vertically aligned there are no forces pushing the particles off vertical. This means that, for instance, there is never a periodic orbit which intersects the z-axis. Looking at the signs, we see that (17) is anti-symmetric in α_x and α_z while (18) is symmetric in α_x . Because the dynamics in this coordinate system are given by analytic functions, the pole and the zeros are always isolated points. Further, away from the pole at the origin the dynamics are differentiable, so that orbits never intersect. Finally, the stationary states are given when the LHS of (17) & (18) are both zero, which provides limits on the count of stationary states. When examining the dynamics it is useful to also write the equations in polar co-ordinates and θ the angle from the x-axis. We use polar rather than cylindrical coordinates because there are no forces in the y-direction. These defined so that

$$\begin{eqnarray}&&{\alpha }_{x}=\alpha \cos (\theta )\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{\alpha }_{z}=\alpha \sin (\theta )\end{eqnarray} \tag{ 20 }$$

where α > 0. This definition entails that the positive x-axis is the ray such that θ = 0, the positive z-axis is $\theta =\tfrac{\pi }{2}$, the negative x-axis is θ = π and the negative z-axis is $\theta =\tfrac{3\pi }{2}$. We could also use $\psi =\tfrac{\pi }{2}-\theta$, the angle from the z-axis. Applying this definition to equation (15), we derive

$$\begin{eqnarray}&&\dot{\alpha }=\displaystyle \frac{\beta }{{\alpha }^{2}}\left(\displaystyle \frac{3}{\alpha }-\displaystyle \frac{{\left(1+\gamma \right)}^{2}}{\gamma }\right)+\left(\displaystyle \frac{3(1-\delta )}{2\alpha }-\displaystyle \frac{(\gamma -\delta )(1+\gamma )}{\gamma }\right)\sin (\theta )\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\alpha \dot{\theta }=\left(\displaystyle \frac{3(1-\delta )}{4\alpha }-\displaystyle \frac{(\gamma -\delta )(1+\gamma )}{\gamma }\right)\cos (\theta )\end{eqnarray} \tag{ 22 }$$

Writing the dynamics in these coordinates also has a few obvious implications worth making explicit. First of all, the stationary states are given when the LHS of (21) & (22) are both zero. Therefore, equation (22) means that any stationary state must have either $\cos (\theta )=0$ (i.e., the line of the particle centers must be parallel with gravity) or have radial coordinate α equal to

$$\begin{eqnarray}&&{\alpha }^{\dagger }=\displaystyle \frac{3(1-\delta )\gamma }{4(\gamma -\delta )(1+\gamma )}\end{eqnarray} \tag{ 23 }$$

Because α > 0, this second option exists only if

$$\begin{eqnarray}&&(1-\delta )(\gamma -\delta )\gt 0\end{eqnarray} \tag{ 24 }$$

We begin our enumerative strategy by considering the case when at least one particle is uncharged and there is no electrostatic interactions between them. Our analysis will demonstrate the importance of non-stable stationary states in establishing the qualitative global dynamics as well as allow us to discuss later the limit of $\beta \to 0$. We will also compare the output of this model to the classical results [46].

We start by giving the dynamical equations. If there is at least one uncharged particle, then β = 0 and the equations of motion (21) & (22) become

$$\begin{eqnarray}&&\dot{\alpha }=\left(\displaystyle \frac{3(1-\delta )}{2\alpha }-\displaystyle \frac{(\gamma -\delta )(1+\gamma )}{\gamma }\right)\sin (\theta )\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\alpha \dot{\theta }=\left(\displaystyle \frac{3(1-\delta )}{4\alpha }-\displaystyle \frac{(\gamma -\delta )(1+\gamma )}{\gamma }\right)\cos (\theta )\end{eqnarray} \tag{ 26 }$$

There are a few trivial cases we will quickly deal with. If δ = 1 then a stationary state only exists if γ = 1 also. In this case, the particles are totally identical and there is no relative motion no matter where the particles begin. Similarly, if δ = γ, either δ = γ = 1 as before or no finite separation can be a stationary state. Either way, this would eliminate any interesting behavior such as bounded orbits. Moreover, only if (24) holds can there be stationary α.

Finally we come to the non-trivial cases, those in which the above inequality is satisfied. Solving for the stationary states gives four solutions: two horizontal and two vertical. The vertical cases haves ${\theta }^{* }=\tfrac{\pi }{2}$ or $\tfrac{3\pi }{2}$ and give

$$\begin{eqnarray}&&{\alpha }^{* }=\displaystyle \frac{3(1-\delta )\gamma }{2(\gamma -\delta )(1+\gamma )}\end{eqnarray} \tag{ 27 }$$

The horizontal cases have ${\theta }^{\dagger }=0$ or π and lie at a distance (23) from the origin. A few points can be established about these stationary states just from equations (23) and (27). First we note that ${\alpha }^{* }=2{\alpha }^{\dagger }$. Next, if δ > 1, then the right hand side of (27) is an increasing function of γ in our range of 0 < γ ≤ 1. Therefore we have that ${\alpha }^{* }\lt 1$ in this case. Stationary states of particles with different radii or different masses can be feasible, i.e., non-overlapping and non-touching, with ${\alpha }^{* }\gt 1$, only if δ < γ < 1.

We illustrate these parameter space results in figure 2(a). The solid line is δ/γ = 1, when the particles have equal Stokes velocities. The short dashed line is equal reduced masses δ = 1. Below this line and above the solid line stationary states do not exist (region D), while above this line stationary states exist but are infeasible (region E). As we have already shown, there is no stationary states on the solid line and short dashed lines except at δ = γ = 1 where every separation is a stationary state.

Zoom In Zoom Out Reset image size

The dash-dot line is α^† = 1, so that above this line but below the solid line (region C) all stationary states are feasible. On the long dashed line is ${\alpha }^{* }=1$, so that above this line but below the dash-dot line (region B) vertical stationary states are feasible but the horizontal ones are not feasible. Below dashed line (region A) all stationary states are not feasible. Properties of the regions are outlined in table 1.

Table 1. Properties of stationary states of uncharged particles.

	Region	Do stationary	Are vertical	Are horizontal
		states exist?	stationary states	stationary states
			feasible?	feasible?
(A)	$\displaystyle \frac{\delta }{\gamma }\leqslant \displaystyle \frac{2\gamma -1}{2-\gamma }$	yes	no	no
(B)	$\displaystyle \frac{2\gamma -1}{2-\gamma }\lt \displaystyle \frac{\delta }{\gamma }\leqslant \displaystyle \frac{4\gamma +1}{4+\gamma }$	yes	yes	no
(C)	$\displaystyle \frac{4\gamma +1}{4+\gamma }\lt \displaystyle \frac{\delta }{\gamma }\lt 1$	yes	yes	yes
(D)	$1\leqslant \displaystyle \frac{\delta }{\gamma }\leqslant \displaystyle \frac{1}{\gamma }\,{\rm{and}}\,\gamma \ne 1$	no	N/A	N/A
(E)	$\displaystyle \frac{\delta }{\gamma }\gt \displaystyle \frac{1}{\gamma }$	yes	no	no

With the locations of the stationary states in mind, we now turn to characterizing the orbits in the large. We can do this by solving exactly a few special cases of the equations of motion. If our initial condition involves α = α* but θ is non-vertical, then the system evolves along the curve given by

$$\begin{eqnarray}&&\dot{\theta }=-\displaystyle \frac{{\left(\gamma -\delta \right)}^{2}{\left(1+\gamma \right)}^{2}}{3{\gamma }^{2}(1-\delta )}\cos (\theta )\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\theta +\displaystyle \frac{\pi }{2}=2{\tan }^{-1}\left[{k}_{1}\exp \left(-\displaystyle \frac{{\left(\gamma -\delta \right)}^{2}{\left(1+\gamma \right)}^{2}}{3{\gamma }^{2}(1-\delta )}t\right)\right]\end{eqnarray} \tag{ 29 }$$

where k₁ is from the initial condition. This solution forms two heteroclinic curves, so designated because each connect two stationary states in infinite time.

The other pair of heteroclinic curves are vertical, ${\theta }^{* }=\tfrac{\pi }{2}$ or $\tfrac{3\pi }{2}$. These heteroclinic are symmetric. In this case the equations of motion are

$$\begin{eqnarray}&&\dot{\alpha }=\pm \left(\displaystyle \frac{3(1-\delta )}{2\alpha }-\displaystyle \frac{(\gamma -\delta )(1+\gamma )}{\gamma }\right)\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&\pm t={k}_{2}-\displaystyle \frac{{\gamma }}{({\gamma }-{\delta })(1+{\gamma })}{\alpha }-\displaystyle \frac{3(1-{\delta }){{\gamma }}^{2}}{2{\left({\gamma }-{\delta }\right)}^{2}{\left(1+{\gamma }\right)}^{2}}{\rm{ln}}\left|{\alpha }-\displaystyle \frac{3(1-{\delta }){\gamma }}{2({\gamma }-{\delta })(1+{\gamma })}\right|\end{eqnarray} \tag{ 31 }$$

where k₂ is from the initial conditions. This form of the solution makes it clear that a particle approaches θ* in infinite time.

Collectively, the curves sketched out by these four heteroclinic orbits make two half circles. The interior of the curves contain two neutrally stable stationary states at (23). Finally, the two vertical stationary states at (27) are not stable but rather saddle points.

We can now fully characterize the orbit traced out by any initial position. By the Poincare-Bendixson theorem every point on the interior of those closed heteroclinic curves will move in periodic motion: moving toward a stable fixed point is ruled out by linear stability analysis, the heteroclinics have all been found and wandering to infinity is impossible because they cannot touch the closed curves which initially bound them. Finally, every point not vertical and exterior to the circle with radius α* will go to infinity, by similar reasoning.

We complete the analysis by returning to figure 2(a) in order to characterize the qualitative dynamics as a function of the parameters. Regions (A)–(C) & (E) have qualitatively similar dynamics: in all cases there is periodic behavior. Region (D) has no stationary states and therefore all orbits are unbounded. Region (C) has the special property that some of the periodic orbits found are non-overlapping. This case is the most physically interesting and the corresponding relative trajectories are illustrated in figure 2(b). This completes the analysis of the qualitative behavior of uncharged point particles settling in a Stokes flow.

We will now briefly compare the results in this section to those in [46]. This seminal paper uses a precise model of hydrodynamic interactions, allowing for more realistic treatment of finitely sized particles. For reference, figure 2(a) can be compared directly to figure 5 in [46].

The essential similarity between the results of both approaches is that there appear such regions of the phase space γ & δ/γ that all the trajectories are unbounded, and such regions that some of the trajectories are unbounded, but the other ones are bounded. In case of the point-particle model, one should focus on the orbits (or their parts) outside the particle 'surface' determined by its radius. For the point-particle model, region (A) and regions (D) and (E) of figure 2(a) contain only unbounded feasible orbits. Regions (D) and (E), where the larger particle moves slower than the smaller one, are the same as the region of unbounded orbits in [46]. Region (A), in which the larger particle moves faster than the smaller one, is smaller than the region of unbounded orbits in [46].

In the region labeled (C) the feasible point-particle orbits are periodic or unbounded. Analogous (but lager) range was also found in [46] for the more precise hydrodynamic model. In the point-particle region (B) there are feasible unbounded orbits and feasible parts of periodic orbits, bounded owing to the non-overlapping vertical stationary states while the horizontal stationary states responsible for periodicity are overlapping. For the more precise hydrodynamic interactions used in [46], the analogous (but smaller) range exists. However, in this range there are no horizontal stationary states and the bounded orbits are not periodic. Rather, the bounded orbits are heteroclinics going from one vertical stationary state to the other. This comparison suggests that the existence of periodic or bounded orbits is predicted by the point-particle and the more precise hydrodynamic model in a qualitatively similar way when all or at least some of the stationary states are feasible.

Following the discussion of equation (22), we have divided all the stationary arrangements into three categories. In the first, the particle centers are aligned vertically and the larger particle is above the smaller particle, i.e. θ = π/2 (see figure 1 for illustration). We will call these larger up vertical stationary states. In the second case, the particle centers are also aligned vertically but the larger particle is below the smaller particle, therefore θ = 3π/2. We will call this second configurations larger down vertical stationary states. Finally, there are the inclined stationary states where θ takes on any other value.

The conditions for stability of vertical stationary states were found in [51]. Here we add a systematic analysis of vertical stationary states ¹ , both stable and not stable. In this section we also investigate a new class of stationary states: the inclined configurations, focusing on their stability.

4.1. Vertical larger up stationary states

We will first analyze properties of vertical stationary states with the larger particle above the smaller one, i.e. θ = π/2. In this case equation (21) becomes

$$\begin{eqnarray}&&{\alpha }^{3}\dot{\alpha }=3\beta -\beta \displaystyle \frac{{\left(1+\gamma \right)}^{2}}{\gamma }\alpha +\displaystyle \frac{3(1-\delta )}{2}{\alpha }^{2}-\displaystyle \frac{(\gamma -\delta )(1+\gamma )}{\gamma }{\alpha }^{3}\end{eqnarray} \tag{ 32 }$$

This equation gives rise to the following condition (equivalent to (14) in [51]) for a stationary state α = α* and $\theta =\tfrac{\pi }{2}$

$$\begin{eqnarray}&&0=6\gamma \beta -2\beta {\left(1+\gamma \right)}^{2}{\alpha }^{* }+3\gamma (1-\delta ){\alpha }^{* 2}-2(\gamma -\delta )(1+\gamma ){\alpha }^{* 3}\,\end{eqnarray} \tag{ 33 }$$

It follows from this condition that it is possible to choose a β to make any particular α* stationary. Therefore, there are stationary states that have no analogue with the uncharged case. We will return to the parameter space (γ, δ) behavior of vertical stationary states in a later section. For instance, in the range 1 > δ > γ, there are no stationary states for uncharged particles, but there are some if a charge is added. We now address linear stability of a stationary state $({\alpha }^{* },\tfrac{\pi }{2})$, given by equation (33). If we assume that $\alpha ={\alpha }^{* }+{\epsilon }_{r}$ and $\theta =\tfrac{\pi }{2}+{\epsilon }_{\theta }$, then the dynamics linear in epsilon are

$$\begin{eqnarray}&&{\dot{\epsilon }}_{r}\approx -\displaystyle \frac{1}{\gamma {\alpha }^{* 3}}\left(\beta {\left(1+\gamma \right)}^{2}-3\gamma (1-\delta ){\alpha }^{* }+3(\gamma -\delta )(1+\gamma ){\alpha }^{* 2}\right){\epsilon }_{r}\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{\dot{\epsilon }}_{\theta }\approx -\displaystyle \frac{1}{4\gamma {\alpha }^{* 2}}\left(3\gamma (1-\delta )-4(\gamma -\delta )(1+\gamma ){\alpha }^{* }\right){\epsilon }_{\theta }\end{eqnarray} \tag{ 35 }$$

Because the matrix of coefficients is diagonal, the constants withing the parentheses must be positive in order for the system to be linearly stable. That is,

$$\begin{eqnarray}&&0\lt \beta {\left(1+\gamma \right)}^{2}-3\gamma (1-\delta ){\alpha }^{* }+3(\gamma -\delta )(1+\gamma ){\alpha }^{* 2}\,\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&0\lt 3\gamma (1-\delta )-4(\gamma -\delta )(1+\gamma ){\alpha }^{* }\end{eqnarray} \tag{ 37 }$$

It can be shown by Lypunov's method that equation (33) and inequalities (36) & (37) form necessary and sufficient conditions for a stationary stable state.

4.2. Vertical larger down stationary states

We will now discuss the stationary states with the larger particle below the smaller one. Similarly to equation (32), on the ray θ = 3π / 2, equation (21) becomes

$$\begin{eqnarray}&&{\alpha }^{3}\dot{\alpha }=3\beta -\beta \displaystyle \frac{{\left(1+\gamma \right)}^{2}}{\gamma }\alpha -\displaystyle \frac{3(1-\delta )}{2}{\alpha }^{2}+\displaystyle \frac{(\gamma -\delta )(1+\gamma )}{\gamma }{\alpha }^{3}\end{eqnarray} \tag{ 38 }$$

This gives rise to the following equations for the stationary state $({\alpha }^{* },\tfrac{3\pi }{2})$

$$\begin{eqnarray}&&0=6\gamma \beta -2\beta {\left(1+\gamma \right)}^{2}{\alpha }^{* }-3\gamma (1-\delta ){\alpha }^{* 2}+2(\gamma -\delta )(1+\gamma ){\alpha }^{* 3}\,\end{eqnarray} \tag{ 39 }$$

and for the local dynamics

$$\begin{eqnarray}&&{\dot{\epsilon }}_{r}\approx -\displaystyle \frac{1}{\gamma {\alpha }^{* 3}}\left(\beta {\left(1+\gamma \right)}^{2}+3\gamma (1-\delta ){\alpha }^{* }-3(\gamma -\delta )(1+\gamma ){\alpha }^{* 2}\right){\epsilon }_{r}\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{\dot{\epsilon }}_{\theta }\approx -\displaystyle \frac{1}{4\gamma {\alpha }^{* 2}}\left[-3\gamma (1-\delta )+4(\gamma -\delta )(1+\gamma ){\alpha }^{* }\right]{\epsilon }_{\theta }\end{eqnarray} \tag{ 41 }$$

In order to be linearly stable, the constants within the parentheses must be positive. That is,

$$\begin{eqnarray}&&0\lt \beta {\left(1+\gamma \right)}^{2}+3\gamma (1-\delta ){\alpha }^{* }-3(\gamma -\delta )(1+\gamma ){\alpha }^{* 2}\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&0\lt -3\gamma (1-\delta )+4(\gamma -\delta )(1+\gamma ){\alpha }^{* }\,\end{eqnarray} \tag{ 43 }$$

We have previously shown in [51] that there are no stable stationary states with the larger particle below the smaller one if we also require the particles do not overlap (α* > 1). However, conditions (42) & (43) will still aid in our understanding of the global dynamics because of the properties of non-stable stationary states.

4.3. Inclined stationary states

In the previous section, we showed how to count and classify the vertical stationary states of charged particles. We will now show how to do the same for non-vertical stationary states, which have a simpler form, given by equations (21), (22) & (23). This simplicity will allow us also to easily analyze typical behaviors in the parameter space (γ, δ).

We can discuss the stationary states of charged particles using figure 2(a) and relations in table 1. The same parameter space curves and regions are of interest but require new interpretations when conceived as relating to charged inclined stationary states. We start from considering two special cases. The short dashed line in figure 2(a) represents equal reduced masses, δ = 1. Equation (22) entails that there exists an inclined stationary state of charged particles only if further δ = γ = 1. Similarly, the solid line, the special case of equal Stokes velocities δ/γ = 1, has an inclined stationary state only if further δ = 1. That is to say, those lines are consistent with the existence of an inclined stationary state only on their intersection point. However, the interpretation of that point is different between the uncharged and the charged inclined case. In the uncharged case, all the relative positions of identical particles are stationary, while in the charged case only stationary states of identical particles must have separation distance ${\alpha }^{\dagger }=\tfrac{3\gamma }{{\left(1+\gamma \right)}^{2}}$ (and an arbitrary orientation). This never feasible distance is the same for all $\beta \ne 0$. Moving on to the general case of unequal Stokes velocities & unequal reduced masses, a general inclined stationary state must have

$$\begin{eqnarray}&&\sin ({\theta }^{\dagger })=\beta \displaystyle \frac{16(\gamma -\delta ){\left(1+\gamma \right)}^{2}(1+3\delta -3\gamma -\delta \gamma )}{9{\gamma }^{2}{\left(1-\delta \right)}^{3}}\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&{\alpha }^{\dagger }=\displaystyle \frac{3(1-\delta )\gamma }{4(\gamma -\delta )(1+\gamma )}\,\end{eqnarray} \tag{ 45 }$$

Equation (45), identical to (23), corresponds to ${\dot{\alpha }}_{\theta }=0$. Therefore, the RHS of (45) cannot be negative, and there are no inclined stationary states of charged particles in the region (D) defined in table 1 and shown in figure 2(a). Moreover, the distance between the particles in an inclined stationary state is independent of the charge/mass ratio β. However, unlike the uncharged case, the condition (44) for ${\dot{\alpha }}_{r}=0$ has a nontrivial form that does not allow for stationary states if charges are too large, because the sine of the stationary angle θ^† is proportional to charge/mass ratio β. In terms of phase diagram as in figure 2(a), equation (45) entails that all and only values of the parameters γ and δ in region (C), given in table 1 have feasible inclined stationary states of charged particles, providing that value of β is sufficiently small to satisfy (44). Horizontal stationary states, generic for the dynamics of uncharged particles, are exceptional for charged systems. One can see from equation (44) that there is a horizontal stationary state with the particle centers aligned perpendicular to gravity (i.e $\sin (\theta )=0$) if and only if either a) the particles are uncharged (β = 0) or b) the particles have Stokes velocity ratio given by the relation $\tfrac{\delta }{\gamma }=\tfrac{3\gamma -1}{\gamma (3-\gamma )}$ (and they overlap with ${\alpha }^{\dagger }=\tfrac{3\gamma }{{\left(1+\gamma \right)}^{2}}$). Because the curve traced by this relation is lower than the dash dot line representing the boundary of feasible inclined stationary states (except at their intersection δ = γ = 1) the angular part of an inclined feasible stable stationary is always less than π, which means that the larger particle is higher than the smaller one.

We now analyze the conditions for linear stability of inclined stationary states. Let $\alpha ={\alpha }^{\dagger }+{\epsilon }_{r}$ and $\theta ={\theta }^{\dagger }+{\epsilon }_{\theta }$ where _r and _θ are first order perturbations. One finds that to a first order

$$\begin{eqnarray}&&{\dot{\epsilon }}_{r}\approx -\displaystyle \frac{3\beta }{{\alpha }^{\dagger 4}}{\epsilon }_{r}+\displaystyle \frac{(\gamma -\delta )(1+\gamma )}{\gamma }\cos ({\theta }^{\dagger }){\epsilon }_{\theta }\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}&&{\dot{\epsilon }}_{\theta }\approx -\displaystyle \frac{(\gamma -\delta )(1+\gamma )}{{\alpha }^{\dagger 2}\gamma }\cos ({\theta }^{\dagger }){\epsilon }_{r}\end{eqnarray} \tag{ 47 }$$

The above linearized dynamics can be easily analyzed. A linear system of ODE with a matrix of constant coefficients is called stable if and only if the real parts of eigenvalues of the matrix are negative. Recall that the determinant is the product of the eigenvalues and the trace is the sum. Therefore, a necessary and sufficient condition for this stationary state to be linearly stable is the determinant to be positive and the trace negative. The determinant

$$\begin{eqnarray}&&\displaystyle \frac{{\left(\gamma -\delta \right)}^{2}{\left(1+\gamma \right)}^{2}}{{\alpha }^{\dagger 2}{\gamma }^{2}}{\cos }^{2}({\theta }^{\dagger })\end{eqnarray} \tag{ 48 }$$

is positive if and only if the particles have different Stokes velocities (recall that $\cos ({\theta }^{\dagger })\ne 0$ for inclined stationary states),

$$\begin{eqnarray}&&\displaystyle \frac{\delta }{\gamma }\ne 1\end{eqnarray} \tag{ 49 }$$

Further, the trace is equal to

$$\begin{eqnarray}&&-\displaystyle \frac{3\beta }{{\alpha }^{\dagger 4}},\end{eqnarray} \tag{ 50 }$$

which is negative if and only if

$$\begin{eqnarray}&&\beta \gt 0\end{eqnarray} \tag{ 51 }$$

Inequalities (49) & (51) entail that if the charges on the particles are opposed then the stationary states are stable whenever they exist: off of the solid and short dashed line in figure 2(a) and with β small enough that equation (44) has a solution. It is interesting to note that the particles having opposite charges, expressed in the inclined case by inequality (51), is also a necessary condition for a vertical stable steady state [51]. Therefore in the next sections we will focus on systems with oppositely charge particles (i.e. β > 0).

The local information derived in the previous section can be used to find the qualitative behavior of pairs of charged particles sedimenting in a Stokes flow. For instance, as just mentioned, we have shown that the local analysis entails that if there is to be a stable stationary state then the charges on the particles must be opposed. In this section we show some generic examples of the dynamics of charged sedimenting particle doublets with large capturing sets. It means that starting from a wide range of initial relative positions, both particles and will not separate from each other, and, on the contrary, they will decrease their distance for ever.

We choose four sets of parameters to demonstrate a variety of behaviors with the above capturing property. In this section we will choose parameters with β > 0 and δ/γ < 1, so that the particles have opposite charges and the larger particle has greater Stokes velocity when the particles are well separated. The dynamics for when δ/γ > 1 will be discussed in the next section. We organize by the count and arrangement of stable stationary states: one stable vertical stationary state, two stable vertical stationary states, two stable inclined stationary states and finally two stable inclined stationary states with one stable vertical stationary state.

Example orbits solving the nonlinear vector ordinary differential equation (15) are plotted in figure 3. Each orbit describes the relative motion of the larger particle 2 with the origin as the center of the smaller particle 1, and with the given initial relative positions of the larger particle. The orbits were calculated using a fourth order Runga-Kutta method with constant step size Δt = .1. We have decided to use (α_x , α_z ) space in this figure and this section for easy comparison with observations. One can see visually that the basins of attraction in our characteristic cases are large. We have also provided in figure 4 a visual way of demonstrating the local stability (or instability) of vertical stationary states. The solid curves are the vertical velocity ${\dot{\alpha }}_{z}$ of particle 2 relative to particle 1 when their line of centers is vertically aligned and the relative vertical position is α_z . On vertical axis, ${\dot{\alpha }}_{z}=\dot{\alpha }$ which can be obtained from (32) & (38). When the solid curve goes down through the horizontal dash-dot line (which corresponds to ${\dot{\alpha }}_{z}=0$) when α_z is increased, then separation between the particle centers is a stationary state which is vertically stable as given by inequalities (36) & (42). The two dashed lines capture the angular stability conditions (37) & (43) for the given parameters, separately for α_z  > 0 and α_z  < 0. In each of these ranges of α_z , if a stationary state is to the left of a dashed line, then it is horizontally stable. In the following we will demonstrate that using properties of the stationary states we are able to determine basic features of the dynamics and qualitatively describe the boundary of the basin of attraction in all of these generic cases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.1. One vertical stable stationary state

5.2. Two stable vertical stationary states

5.3. Inclined stable stationary states

5.4. Dynamics with both inclined and vertical stable stationary states

In all the cases of the previous section, when particle separation is large the larger particle moves in the direction of gravity relative to the smaller particle. This is because in all the cases examined the ratio of Stokes Velocities δ / γ < 1. It is also possible to have δ / γ ≥ 1, so that the larger particle moves against gravity relative to the smaller particle when the particle separation is large. We will call this case that of 'inverted Stokes velocity ratio'. The stability considerations from section 4 do not change in form for Inverted Stokes velocity ratios. We will see by example that the Poincare-Bendixson theorem used throughout the previous section still allows us to make qualitative conclusions about the global dynamics.

As our example, we choose parameters β = .293... δ = .82 & γ = .8. This corresponds to the range (D) of the phase space γ & δ/γ, see figure 2(a) and table 1. The corresponding orbits of the larger particle (with the label 2) relative to the smaller particle (with the label 1) are illustrated in figure 5(a). One can see that even though the stable stationary state is still 'larger up' (α_z > 0), the part of the capturing set when the larger particle is above the smaller particle is now bounded, whereas in figure 3 it was not bounded. The part of the capturing set with the larger particle below the smaller particle is now unbounded, whereas before it was bounded. This property of the capturing set corresponds to the inverted direction of the relative trajectories, which have now vertical components coming from $-\infty$ while in the previous section, they arrived from $+\infty$, in agreement with the inverted Stokes velocity ratio.

Zoom In Zoom Out Reset image size

We can now extract qualitative information about the dynamics of this system. The dynamics when the particle centers are aligned with gravity are displayed in figure 5(b). The stationary state at ${\alpha }_{z}^{* }=1.01...$ is stable and the stationary states at ${\alpha }_{z}^{* }=2.12$ and ${\alpha }_{z}^{* }=-0.63$ are saddle points. There are no inclined stable stationary states. The separatrix coming into the larger up saddle point must come in from infinity. This is due to the Poincare-Bendixson theorem: the only other place the separatrix could begin is the larger down saddle point, but this would require a closed curve of orbits without a steady state in the interior. The heteroclinics coming out of the larger down saddle point must go to the stable stationary state. This is also an application of the the Poincare-Bendixson theorem: the orbit cannot diverge without crossing over the heteroclinic coming into the larger up saddle point and cannot be closed because there is no inclined stationary state to be on a closed orbit's interior. Therefore it must end in a stationary state and the only choice is the sole stable one. The curve ending in the larger saddle point is the boundary of the capturing set. Almost all the orbits on the interior of these curves end up on the stable stationary state. All the orbits on its exterior diverge, which corresponds to separation of the particles.

In section 3, we gave a phase diagram of the potential stationary states of the pair of uncharged point-particles as a function of the ratio δ/γ of Stokes velocities and the ratio γ of particle radii. Further, in subsection 4.3, we demonstrated how to reinterpret this graph in the case of inclined stationary states of charged particles. In this section we will give the phase diagram appropriate to the case of 'larger up' vertical stationary states, including a heatmap with information about the separation of particle centers at the stable stationary state. We concentrate on such stationary states because we have already shown that 'larger down' vertical stationary states are never stable [51]. We will also compare the phase diagrams and heat maps for vertical and inclined stationary states.

In the remainder of this section we will consider particles of the opposite charges,

$$\begin{eqnarray}&&\beta \gt 0,\end{eqnarray} \tag{ 52 }$$

because only in this case there exist stable stationary states, as shown in [51]. We will also assume, again following [51], that at the stationary state the particles are non-overlapping (feasible), that is,

$$\begin{eqnarray}&&1\lt {\alpha }^{* }\end{eqnarray} \tag{ 53 }$$

This inequality along with equation (33) and inequalities (36) & (37) form necessary and sufficient conditions for a 'larger up' feasible locally asymptotically stable vertical stationary state.

To complete the list of the bounds imposed on the parameters, we remind that

$$\begin{eqnarray}&&0\lt \gamma \leqslant 1,\quad \delta \gt 0.\end{eqnarray} \tag{ 54 }$$

Our goal is to identify the characteristic regions in the parameter space of δ/γ and γ appropriate for feasible vertical 'larger up' stable stationary states of charged point particles, for certain values of β & α*. In fact, we will start by solving equation (33) for β as a function of α*

$$\begin{eqnarray}&&\beta =\displaystyle \frac{3\gamma (1-\delta )-2(\gamma -\delta )(1+\gamma ){\alpha }^{* }}{2{\left(1+\gamma \right)}^{2}{\alpha }^{* }-6\gamma }{\alpha }^{* 2}\end{eqnarray} \tag{ 55 }$$

The conditions (52) and (55) give the following bound on δ/γ as a function of α* and γ,

$$\begin{eqnarray}&&\displaystyle \frac{{\delta }}{{\gamma }}\gt \displaystyle \frac{2(1+{\gamma }){{\alpha }}^{\ast }-3}{2(1+{\gamma }){{\alpha }}^{\ast }-3{\gamma }}.\end{eqnarray} \tag{ 56 }$$

Because the right hand side of (56) is an increasing function of α*, the greatest lower bound for the ratio of Stokes velocities of oppositely charged particles at a feasible stationary state is reached at α* = 1. This bound is

$$\begin{eqnarray}&&\frac{{\delta }}{{\gamma }}\gt \displaystyle \frac{2{\gamma }-1}{2-{\gamma }}.\end{eqnarray} \tag{ 57 }$$

The bound $\delta /\gamma =2\gamma -1/2-\gamma$ is plotted by dashed line (green online) in figure 6(a). Below this line, in the range (A), there are no feasible stationary states.

Zoom In Zoom Out Reset image size

We now analyze angular stability of stationary states. We write inequality (37) as

$$\begin{eqnarray}&&\displaystyle \frac{\delta }{\gamma }\gt \displaystyle \frac{4(1+\gamma ){\alpha }^{* }-3}{4(1+\gamma ){\alpha }^{* }-3\gamma }\end{eqnarray} \tag{ 58 }$$

The minimum of the right hand side, reached at α* = 1, is the greatest lower bound in the phase space (δ/γ, γ) for existence of feasible stationary states stable against angular perturbations,

$$\begin{eqnarray}&&\displaystyle \frac{\delta }{\gamma }\gt \displaystyle \frac{4\gamma +1}{4+\gamma }\end{eqnarray} \tag{ 59 }$$

The bound is plotted by dash-dot line (blue online) in figure 6(a). Below this line, in the range (B), all feasible stationary states are angularly unstable. We now analyze radial stability of stationary states. We use equation (55) to eliminate β from (37). After a careful analysis of signs we obtain

$$\begin{eqnarray}&&\displaystyle \frac{\delta }{\gamma }\lt P(\gamma ,{\alpha }^{* })\end{eqnarray} \tag{ 60 }$$

where

$$\begin{eqnarray}&&P(\gamma ,{\alpha }^{* })=\displaystyle \frac{18\gamma +4{\left(1+\gamma \right)}^{3}{\alpha }^{* 2}-3(1+\gamma )(1+7\gamma ){\alpha }^{* }}{18{\gamma }^{2}+4{\left(1+\gamma \right)}^{3}{\alpha }^{* 2}-3\gamma (1+\gamma )(7+\gamma ){\alpha }^{* }}.\end{eqnarray} \tag{ 61 }$$

It leads to the following least upper bound in the phase space (δ/γ, γ) for existence of stationary states stable against radial perturbations,

$$\begin{eqnarray}&&\frac{{\delta }}{{\gamma }}\lt R({\gamma }),\end{eqnarray} \tag{ 62 }$$

where

$$\begin{eqnarray}R(\gamma )=\left\{\begin{array}{lll}1 & {\rm{for}} & \gamma \leqslant 2-\sqrt{3}\\ P(\gamma ,1) & {\rm{for}} & \gamma \gt 2-\sqrt{3}\end{array}\right.\end{eqnarray} \tag{ 63 }$$

The bound (62)–(63) is plotted in figure 6(a) as the sum of the boundaries between regions (C) and (H) & regions (G) and (H). Above it, in the range (H), all stationary states are radially unstable.

Basic properties of 'larger up' vertical stationary states in the regions (A)-(H) are outlined in table 2. Recall that the ratio of forces β > 0 has been eliminated using the stationary condition (33) and so remarks about stability should be interpreted with the β so derived. For δ/γ < 1 the requirement of the opposite charges (56) and the angular stability condition (58) give the upper bounds on the interparticle distance at a stationary state, ${\alpha }^{* }\lt 2M(\delta /\gamma ,\gamma )$ and ${\alpha }^{* }\lt M(\delta /\gamma ,\gamma )$, respectively, with

$$\begin{eqnarray}&&M(\delta /\gamma ,\gamma )=\displaystyle \frac{3\gamma (\tfrac{1}{\gamma }-\tfrac{\delta }{\gamma })}{4(1-\tfrac{\delta }{\gamma })(1+\gamma )}.\end{eqnarray} \tag{ 64 }$$

Table 2. Properties of 'larger up' non-overlapping vertical stationary states (each with a certain value of β > 0).

	Region	Properties of stationary states
(A)	$0\lt \displaystyle \frac{\delta }{\gamma }\leqslant \displaystyle \frac{2\gamma -1}{2-\gamma }$ or γ = 1	There are no stationary states with ${\alpha }^{* }\gt 1$.
(B)	$\displaystyle \frac{2\gamma -1}{2-\gamma }\lt \displaystyle \frac{\delta }{\gamma }\leqslant \displaystyle \frac{4\gamma +1}{4+\gamma }$	For any $1\lt {\alpha }^{* }\lt 2M(\delta /\gamma ,\gamma )$ there exists a stationary state. There are no stationary states with ${\alpha }^{* }\gt 2M(\delta /\gamma ,\gamma )$. Each stationary state is unstable with respect to angular perturbations.
(C)	$\displaystyle \frac{4\gamma +1}{4+\gamma }\lt \displaystyle \frac{\delta }{\gamma }\lt 1$	For any $1\lt {\alpha }^{* }\lt 2M(\delta /\gamma ,\gamma )$ there exists a stationary state. There are no stationary states with ${\alpha }^{* }\gt 2M(\delta /\gamma ,\gamma )$. For $1\lt {\alpha }^{* }\lt M(\delta /\gamma ,\gamma )$ there exist stable stationary states. For $M(\delta /\gamma ,\gamma )\leqslant {\alpha }^{* }\lt 2M(\delta /\gamma ,\gamma )$ stationary states are angularly unstable.
(G)	$1\leqslant \displaystyle \frac{\delta }{\gamma }\lt R(\gamma )$ and $2-\sqrt{3}\lt \gamma \lt 1$	For any ${\alpha }^{* }\gt 1$ there exists a stationary state. Each stationary state is stable against angular perturbations. There exist stable stationary states.
(H)	$R(\gamma )\leqslant \displaystyle \frac{\delta }{\gamma }$	For any ${\alpha }^{* }\gt 1$ there exists a stationary state. Each stationary states is radially unstable and angularly stable.

In general, for a given β, γ and δ/γ there might be several stable 'larger up' vertical stationary states with different values of the distance α*. We focus on the largest of them. We find the upper bound on the separation of particle centers, which is graphed in figure 6(b), where the regions of increasingly bright color correspond to regions of increasing separations. The white space corresponds to regions which do not have feasible stable stationary states. Analytically, we choose α_max as the least upper bound of the α* which satisfy inequalities (58) & (60) and we derive the following relation in region (C),

$$\begin{eqnarray}&&\displaystyle \frac{\delta }{\gamma }=\displaystyle \frac{3-4(1+\gamma ){\alpha }_{\max }}{3\gamma -4(1+\gamma ){\alpha }_{\max }}\end{eqnarray} \tag{ 65 }$$

In region (G) this relation is

$$\begin{eqnarray}&&\displaystyle \frac{\delta }{\gamma }=\max \left[1,P(\gamma ,{\alpha }_{\max })\right]\end{eqnarray} \tag{ 66 }$$

For comparison, we graph in figure 6(c) the separations of particles for inclined stationary states given by equation (45). As we previously noted, these are feasible only in region (C). One must take care during interpretation because what is plotted in this figure is the exact value rather than merely an upper bound. One can see that for inclined equilbria the lines of constant interparticle separation are monotonically functions of γ, that is that it takes an increasing ratio of Stokes velocities to balance a system with increasing reduced masses at the same separation distance.

We have shown that coexistence of hydrodynamic and electrostatic interactions between particles sedimenting in Stokes flows leads to the dynamics essentially different than in the absence of charge or in the absence of fluid. Using the point-particle model, we demonstrated analytically that charged particles can form stable doublets with basins of attraction in the space of the particle relative positions which are very large in comparison to particle radius. This result indicates that charged sedimenting particles can capture one another, even if the initial distance between them is large. The captured particles tend to a certain stable relative position, where the distance between the particle centers is larger than the sum of their radii, so that their surfaces are separated from each other by a fluid. In particular, there exist stable stationary configurations of charged particles separated by large distances, with large basins of attraction, within the range where the hydrodynamic interactions can be approximated as point-like.

Moreover, even if the ratio β of electrostatic to gravitational forces is very small, the dynamics of pairs of charged particles is both quantitatively and qualitatively different from the dynamics in the absence of any electrostatic interactions. The main qualitative difference is the structure of the capturing set in the space of relative positions. For charged particles, the capturing set consists of trajectories tending to a stationary state while for uncharged particles, it consists of neutrally stable periodic orbits.

The existence of a capturing set in the space of the relative positions is associated with the existence of stable stationary configurations of two charged sedimenting particles. However, stable stationary states are formed only for certain ranges of the ratio of particle radii γ and the ratio of Stokes velocities δ/γ. Therefore, we have determined the region in the parameter space of γ and δ/γ where stable stationary configurations exist for certain values of β, and for a certain range of values of the distance α* > 1 between the particles at the stable stationary state. In this way, we have shown that the capturing of charged particles takes place in a large region in the parameter space of γ and δ/γ. Interestingly, for some values of β, γ and δ/γ there exist multiple stable stationary states inside the capturing set.

We have found stable stationary states with the line-of-centers at an angle ψ inclined with respect to gravity. In the point-particle model, $\cos (\psi )$ is proportional to the ratio β of electrostatic to gravitational force while the particle-to-particle separation distance α^† at the inclined stationary state does not depend on β.

By analyzing examples of capturing sets, we have found that the basin of attraction of all stable stationary states has a boundary which is a surface of revolution of a trajectory which ends on a saddle point stationary state with the particle centers vertically aligned. It seems that the difference between vertical positions of the particles at this stationary state can be used as an estimate of 'the cross-section' of the capturing set. For example, in figure 3(b) and 3(d) the cross section is more than 10 times the sum of the particle radii ${a}_{1}+{a}_{2}$, and particle at the stable stationary configuration are separated by more than 4 (a₁ + a₂). For large capturing sets and large particle separations, the particle dynamics is well-approximated by the point-particle model. Therefore, the analysis presented here suggests that the existence of stable doublets could be confirmed by future experiments.

This work was supported in part by the National Science Centre under grant UMO-2018/31/B/ST8/03640.

All data that support the findings of this study are included within the article (and any supplementary files).