Clustering-induced self-propulsion of isotropic autophoretic particles

Abstract

Self-diffusiophoretic particles exploit local concentration gradients of a solute species in order to self-propel at the micron scale. While an isolated chemically- and geometrically-isotropic particle cannot swim, we show that it can achieve self-propulsion through interactions with other individually-non-motile particles by forming geometrically-anisotropic clusters via phoretic and hydrodynamic interactions. This result identifies a new route to symmetry-breaking for the concentration field and to self-propulsion, that is not based on an anisotropic design, but on the collective dynamics of identical and homogeneous active particles. Using full numerical simulations as well as theoretical modelling of the clustering process, the statistics of the propulsion properties are obtained for arbitrary initial arrangement of the particles. The robustness of these results to thermal noise, and more generally the effect of Brownian motion of the particles, is also discussed.

---

1 Introduction

Achieving self-propulsion through fluids at the micron scale entails many challenges associated with the dynamics of the surrounding fluid medium, including overcoming the dominant effect of viscous dissipation that precludes any inertia-related history effect, and breaking temporal and spatial symmetries in the flow forcing.^1,2 Swimming micro-organisms (e.g. bacteria or swimming algae) achieve self-propulsion using the propagation of waves along their cilia or flagella, in order to generate periodic yet non-reciprocal strokes.^3–5 These biological systems have inspired the experimental design of many artificial swimmers in the lab.^6–9 However, such mechanical systems present fundamental and practical limitations, such as moving parts or their reliance on an external, often macroscopic, actuation (e.g. unsteady magnetic field).

Catalytic colloids represent another category of artificial micro-swimmers that rely instead on self-generated gradients in the physico-chemical properties of their immediate environment.¹⁰ Self-diffusiophoresis is such an example, where a chemically-active colloid moves in response to a concentration gradient of a solute species that is either produced or consumed at its surface through a catalytic reaction.^11,12 Other phoretic systems rely on self-generated gradients of temperature (thermophoresis) or electric potential (electrophoresis) to self-propel.^13,14 Such autophoretic motion fundamentally requires two distinct physico-chemical properties of the particle. Short-range interactions between its surface and the surrounding solute generate an effective slip motion of the fluid in response to gradients in surface concentration, a property generally termed as phoretic mobility, .¹⁵ The second property, the chemical activity, , refers to the particle's ability to absorb or release solute through chemical reactions at their surface.¹⁶

A major advantage of such active colloids in comparison with macroscopically-actuated microswimmers lies in the local nature of the interaction phenomena at the heart of their propulsion, which allows for complex collective behaviour; in contrast, collective motion of magnetic swimmers is essentially driven by a macroscopic forcing that is uniform at the scale of the inter-swimmer distance. Experimentally, these active colloids have been observed to reach speeds ranging from a few to tens of diameters per second.^17–19 Promising experimental results and theoretical predictions have motivated research into possible biomedical and engineering applications such as cargo transport for targeted drug delivery^20,21 and micromachines.²²

These colloids are fundamentally out of equilibrium as they continuously convert physico-chemical energy of their environment into mechanical work. Hydrodynamic and chemical interactions of such self-propelled particles may lead to complex collective dynamics including clustering^23–25 or richer patterns,^26–29 which can often be linked to instabilities of homogeneous suspensions.^30,31 Individually non-motile particles can further achieve self-propulsion by forming stable chemically-inhomogeneous clusters,^32,33 a phenomenon which was recently also observed in experiments.^34,35

To achieve self-propulsion, autophoretic colloids must set the surrounding fluid into motion which fundamentally requires an asymmetric distribution in the solute concentration at their surface. So far, three different mechanisms have been identified to achieve such symmetry-breaking of the concentration field: (i) an asymmetric chemical patterning of the surface (e.g. Janus particles,^16,19) (ii) an asymmetric shape of the chemically-homogeneous colloid^36,37 and (iii) an instability resulting from the non-linear advective coupling between the solute dynamics and the flow motion.^38,39 The former two are fundamentally associated with the particle design, and are built into its architecture. The latter, in contrast, arises spontaneously from the destabilization of a non-motile isotropic state.

The purpose of the present work is to introduce and characterize a fourth route to symmetry-breaking of the concentration field and self-propulsion, wherein geometrically- and chemically-isotropic particles interacting with diffusive solutes do not propel on their own but instead gain locomotion by forming asymmetric clusters. All particles are identical here, which is a fundamentally different situation from the assembly of chemically-inhomogeneous molecules from a mixture of two different types of particles.^32,33 Identical isotropic phoretic particles which can attract (for and of opposite signs) each other because they generate radial concentration gradients which induce a phoretic drift on the other particles. This phoretic attraction combined with steric constraints enable only a discrete set of stable configurations that may display a geometric asymmetry (Fig. 1), which is a sufficient ingredient for self-propulsion of this assembly.³⁷

Fig. 1 Stable planar arrangements of N = 6 attracting phoretic particles. Colorbar shows the relative surface concentration. The asymmetry of the geometry sets whether the cluster translates (left) remains stationary (center) or rotates (right).

This collective self-propulsion is intimately linked to the exact geometry of the particle assembly, which itself results from the dynamic phoretic clustering of the particles. An essential goal of the present work is therefore to characterize the statistical properties of the particle arrangement arising from the clustering process, and therefore requires a careful modelling of this dynamics.

Note that self-propulsion and collective dynamics of chemically-isotropic and individually non-motile particles was also observed for colloidal particles trapped at a fluid–fluid interface;⁴⁰ then fluid motion resulted from the Marangoni stresses at the free surface rather than a direct hydrodynamic forcing by each particle as considered here.

The general problem for N particles is first presented in Section 2. In Section 3, the collective dynamics of two identical isotropic particles is considered in detail and their relative motion is computed analytically. Such analytic solutions are not available for larger number of particles and Section 4 focuses on the modeling of the clustering dynamics and cluster self-propulsion. A reduced-order model, validated with full numerical solution using a regularized Boundary Element Method (BEM), is then used in Section 5 to determine the statistics of formation of different clusters, their velocity and the resulting mean properties and their evolution with N. Section 6 analyses the effect of noise and Brownian motion on these results and conclusions are finally drawn in Section 7.

2 Collective dynamics of N isotropic phoretic particles

The dynamics of N identical spherical particles is considered in an unbounded fluid of viscosity η and density ρ. The surface of the particles interacts with a chemical solute suspended in the fluid phase; following the classical continuum framework,^15,16,41,42 this local interaction results in an effective slip velocity where C(x) is the surface solute concentration, n the local normal unit vector to the surface and the phoretic mobility, which is a signed characteristic of the solute–particle interaction. The second physico-chemical property of the particles' surface is its activity, namely its ability to alter the solute concentration (e.g. through catalytic reactions). For simplicity, this activity is modeled here as a fixed-flux emission () or absorption () of solute, but this framework could easily be extended to more complex chemical kinetics.^42–45 The solute diffuses in the fluid domain with diffusivity κ, and its far-field concentration is noted C_∞. The particles considered here are geometrically and chemically identical and homogeneous, i.e. and are constant and uniform for all spherical particles, which are also of the same radius.

The radius a of the particles is chosen as the characteristic length scale. The relative solute concentration is then defined as , and characteristic velocity and pressure are chosen as and , respectively. In the following, all quantities are non-dimensionalized using these characteristic scales.

We assume that the particles are small enough so that inertial and advective effects are negligible on the fluid and solute transports (i.e. the Reynolds number, and Péclet number, are both small). Thus, in this purely diffusive limit, the solute concentration field obeys the steady state diffusion equation,

∇²c = 0, (1)

with boundary conditionswhere is the dimensionless surface activity on the surface of each particle j. Since advection of the solute by the fluid is neglected here, the solute dynamics can be solved for independently and the fluid dynamics problem can then be obtained in a second step using the dimensionless Stokes flow equations,

∇²u = ∇p, and ∇·u = 0, (3)

with boundary conditions

ũ_j = M(I − nn)·∇c and u(r → ∞) → 0, (5)

where, is the dimensionless mobility, R_j is the position of the center of particle j, and U_j = Ṙ_j and Ω_j are its translation and rotation velocities. The latter quantities are determined uniquely by imposing that each particle j remains force- and torque-free, i.e. for each j,with σ = −pI + (∇u + ∇u^T), the stress tensor. F_j and T_j are the contact (e.g. steric) forces applied on particle j in a cluster. When the particles are independent, they are force- and torque-free, i.e.F_j = 0, T_j = 0 for each value of j. When the particles form a cluster, their overlap is prevented by steric and other short-ranged effects. Thus, in this case, the particles are individually not force- and torque-free, but rigidly-bound. However, the system of equations can be closed by imposing that the system as a whole is force- and torque-free . In the following, we assume that steric interactions between particles generate only central forces so that T_j = 0.

Indeed, for spherical and isotropic particles with negligible solute advection, a rotation of any particle around its center leaves the solute concentration unchanged. In the following, we therefore solely focus on the translation velocities of the particles, U_j (note that Ω_j may not be zero, but is set to satisfy the condition that each particle remains torque-free).

For a single isolated particle (N = 1), the concentration field is purely isotropic c(r) = A/r and does not generate any surface gradient (or slip velocity) along its surface; a single isotropic phoretic particle is therefore unable to self-propel despite its chemical activity (see Fig. 2a).

Fig. 2 Relative concentration field around (a) a single isotropic phoretic particle and (b) two identical isotropic phoretic particles (A = 1 and M = −1). In (a), the isotropic concentration field leads to no fluid motion nor propulsion. In (b), the asymmetric concentration field around each particle leads to a mutual attraction.

However, the symmetry of the concentration field is broken by the presence of another identical particle (Fig. 2b). The surface concentration gradients lead to a mutual phoretic attraction (or repulsion, depending on the relative sign of their surfaces' properties, A and M) between the particles along their axis of symmetry, e_z. The resulting phoretic interactions can be computed and characterized analytically for two particles, and this is the focus of the next section.

3 Phoretic clustering of two identical particles

3.1 Intuitive model

To understand symmetry-breaking and the resulting clustering, we first take a look at the far-field limit when the distance between the particles d is much greater than their radius (a/d ≪ 1). In this limit, in addition to its own chemical field, particle i is exposed at leading order to the source field A/r_j of particle j ≠ i, which creates a locally uniform concentration gradientwith d = |R₂ − R₁|. In response to this gradient, particles experience equal and opposite phoretic drift velocities along their axis of symmetry¹⁵with e_z = (R₁ − R₂)/d. Phoretic attraction (resp. repulsion) therefore arises for AM = −1 (resp. AM = 1). We shall henceforth focus exclusively on the attractive case (AM = −1) which will be responsible for clustering. When the far-field assumption does not hold, this velocity must be corrected to account for confinement effects in both the diffusion and hydrodynamic problems. It is yet amenable to exact analytical calculation using classical results for Laplace's and Stokes' equations in bispherical coordinates.^37,46,47

3.2 Analytical solution for arbitrary separation

The clustering dynamics of two particles is solved analytically for arbitrary distance d using bispherical coordinates (τ, η, ϕ) defined from the cylindrical coordinate system (ρ, ϕ, z) asA constant value of τ corresponds to a sphere of radius α/sinh|τ| with center located at α/tanh τ. Hence, the unit vectors, e_τ and e_μ are respectively, normal and tangential to a spherical surface. Here both spheres have the same unit radius, and correspond to τ = ±τ₀. The positive constants α and τ₀ are determined so that

d = 2α cosh τ₀, α = sinh τ₀. (11)

The general solution to Laplace's equation for concentration c, that decays in the far-field, is obtained within this framework⁴⁶ aswhere L_n(μ) is the Legendre polynomial of degree n. Enforcing the surface flux boundary condition, eqn (2),determines the coefficients C_n uniquely.³⁷ The surface slip velocity is then determined from the concentration distribution on the surface,³⁷ and the reciprocal theorem for Stokes flow is used to compute the axial velocities ∓U of particles 1 and 2:^37,48where σ* is the hydrodynamic stress tensor of the auxiliary flow field corresponding to the solution of Stokes equations around rigid spheres and translating with velocities ±U*e_z, in response to an outer force ∓F*e_z. This auxiliary flow field can be expressed aswith the streamfunction ψ* given by^37,46The no-slip boundary conditions on the spheres' surface write asand appropriate projections on the Legendre polynomials provide an explicit determination for β_n and γ_n. The shear stress distribution on the surface is then determined as³⁷and the hydrodynamic force on each sphere is obtained directly as⁴⁶

Note that, due to action–reaction symmetry, this mutual attraction/repulsion velocity depends only on the distance d between the particles. Its value is hence computed numerically at various separation, d_c = d − 2a and is plotted in Fig. 3.

Fig. 3 Clustering velocity of a catalytic particle of unit radius. The solution converges to 1/d² decay, as expected in the far-field (d_c → ∞); in the case with repulsion velocity, C = 35, δ = 25 and d* = 1.95 were chosen in eqn (22).

The phoretic slip model,¹⁵eqn (5), at the heart of the present modeling likely breaks down when d_c ∼ λ with λ the thickness of the interaction layer where solute–particle interactions are significant. Yet, the phoretic slip model remains valid everywhere except within a small lubrication region between the particles which is only expected to introduce a subdominant correction to the present predictions.

The clustering velocity converges to a unit value when the particles are in contact. Hydrodynamic (lubrication) and phoretic effects are thus insufficient to avoid particle contact in finite time;⁴⁹ short-ranged interaction forces must therefore be added to account for intermolecular repulsive forces between the interacting surfaces preventing particle overlap (i.e. steric effect). To preserve the previous velocity formulation, this repulsive effect is therefore included as an additional repulsive velocity to mimic this behaviour, in the form of a smoothed step-function. The complete expression for the clustering velocity is hence given by

U_a(d) = U(d) − C(1 − tanh(δ(d − d*))), (22)

where C, δ and d* are chosen appropriately so that the particles have negligible separation after clustering (d_c ≪ a). This clustering velocity is a function solely of the particles' distance, shown in Fig. 3, and the dynamics of the two particles can therefore be described by an over-damped deterministic Langevin dynamics within an interaction potential:andThis potential is a measure of the stability of the two-particle cluster, and the system evolves at each instant down the direction of its steepest gradient.

For two particles, hydrodynamic and phoretic interactions introduce an asymmetry in the system that enables the motion of individual particles. Yet, the center of mass of the arrangement remains stationary for identical particles, since the final cluster (i.e. dimer) is front–back symmetric. For identical particles, the emergence of global self-propulsion is therefore restricted to N > 2.

4 Modeling the clustering and collective motion of N particles

While an analytical solution was amenable for N = 2 particles, this is not the case for larger N. Analytical methods, based for example on the Method of Reflections (MoR), can still be used to provide asymptotic approximations in particular for large separation distances between the particles (see Appendix B). Such approximations provide valuable qualitative and quantitative insight, but they fail to provide a description of the near-field interactions. These could nevertheless be studied using asymptotics through lubrication analysis.⁴⁹ Yet, in order to obtain a complete description of both the chemical and hydrodynamic problems for arbitrary arrangements, it is necessary to resort to numerical simulations of the Laplace and Stokes problems. In problems involving fluid–solid interactions, commonly employed numerical techniques include finite element and immersed boundary methods⁵⁰ and their variations.^51,52

Boundary element methods (BEM), based on classical boundary integral representations of the flow field, have also found an important niche in solving Stokes flow problems.⁵³ However, their implementation for phoretic problems has only been recent.^54,55 Boundary element methods solve the Laplace (1) and Stokes (3) equations outside a set of rigid particles by using the fundamental integral representation of the solutions to these harmonic and bi-harmonic equations, in terms of their values and normal gradients on the bounding surfaces alone. Such methods are therefore particularly well-suited for phoretic problems in Stokes flows since the coupling of the concentration and velocity fields occurs only on the particles' surfaces.

Because they represent the solution as the superposition of fundamental singularities (e.g. sources and stokeslet), classical boundary element methods for Laplace and Stokes problems involve singular kernels which require a separate analytical treatment of the singular contributions and precise quadrature techniques. An alternative is to regularize the singularities, for example using the method of regularized stokeslets for Stokes flows.^56–58

4.1 Regularised boundary element method

The concentration as well as the flow fields, being harmonic functions, can be evaluated at any point in the domain using boundary integral representations. BEMs are numerical approximations to these surface integrals obtained by discretising the surface into a collection of elements. The surface concentration field is evaluated by distributing regularized sinks and source dipoles on the surface of the particle while the surface flow velocities are computed by distributing regularized stokeslets and stresslets. These regularised singularities are collocated at discrete points on the surface (called nodes) which are element vertices.

We follow the computational framework of regularised BEM developed for phoretic problems by Montenegro-Johnson et al.⁵⁵ to generate mesh and quadrature routines for evaluation of surface integrals. The surface of each phoretic particle is discretised into 3072 piecewise-quadratic triangular elements (1538 linear nodes). Surface unknowns are discretized as linear functions over each element for accurate computation of the concentration field and traction on the surface. The method begins with the regularized boundary integral equation for the concentration field in response to the flux forcing on the active particles⁵⁵

where x is any point on the surface of the spherical particles with local normal n(x) and x₀ is the point where concentration is to be calculated. The regularised blob function and associated kernels are given byandwith r = |x − x₀| and r_ε² = r² + ε². G^ε and K^ε_j represent a regularised sink and source dipole respectively. The essence of the regularized boundary method is to express simply the left-hand side of eqn (25) in terms of the local concentration c(x₀) so as to obtain an integral equation to solve for c(x₀) on the boundary of the particles. When x₀ is located on the surface of a particle, the left-hand side of eqn (25) can be expressed as (see Appendix A)with κ the mean local curvature of the particle surface (here κ = 1 for spherical particles). It should be noted that the above result is true for surfaces of arbitrary shape and that the effect of the regularization introduces two different O(ε) corrections to the classical c(x₀)/2 result obtained for singular boundary integral methods. These corrections are respectively proportional to the concentration and normal flux, which can be easily implemented in classical boundary integral methods frameworks.† A similar correction should be implemented when x₀ is not exactly on the surface. It can however be demonstrated that these corrections scale at leading order as O(ε⁴/d⁴) and O(ε³/d²) (here d is the distance to the integration surface) and are therefore negligible except in the immediate vicinity of the surface.

Once the surface concentration field has been obtained, the surface slip velocity, ũ, is computed from the particle's phoretic mobility property. The boundary integral formulation of the Stokes flow problem and force- and torque-free conditions on each particle is then solved for the flow traction and translation and rotation velocities U_j and Ω_j of the different particles using^55,56

where f is the surface traction and the flow velocity at the surface, u = ũ + U + Ω × x. Here, the Green's functions S^ε and T^ε represent a regularised stokeslet and its associated stress tensor, respectively:If one is only interested in the resulting motion of the particles, such a boundary integral approach is particularly well-suited as no bulk quantity needs to be computed and remeshing the domain at each time step is unnecessary. The integral in the left-hand side of eqn (29) can be treated similarly as for the Laplace problem. Only if necessary (e.g. for plotting purposes), the surface distribution of concentration and velocity may be used to compute their bulk distribution, using the fundamental integral representations.

For simplicity, we restrict the following analysis to two-dimensional clustering and self-propulsion of the spherical particles. However, we note that the solute diffusion and hydrodynamic problem remain three-dimensional and unbounded. A variable mesh size is implemented along both polar and azimuthal directions of each sphere such that the mesh is refined within and near the plane of clustering. Further, in stable clusters, particles are arranged on a regular hexagonal lattice. Additional mesh refinement is then performed in the clustering plane near the positions of particles' contact on hexagonal lattices. This meshing provides an efficient framework to compute the swimming velocity of clustered particles (see Fig. 4).

Fig. 4 An example of the mesh used for regularised BEM computations (coarsened to 1/4th the total number of elements used so as to aid visualisation). The surface concentration distribution on the 5-particle cluster obtained from the simulations are also shown.

The regularization parameter ε must be chosen small enough to approach the true solution but large enough to avoid a singular behaviour of the integral equation; a value of 0.005 and 0.01 is typically used. This regularized boundary element method was shown to provide accurate computations for near- and far-field dynamics of various phoretic problems.⁵⁵ We validated here the method on the two-particle problem for which an analytical solution was obtained. Using the above mentioned meshing, an accuracy of about 1% was obtained on the swimming velocity of the cluster even for particles near contact with d_c = 0.01 and d_c = 0.02 (about 10 times the local element size in the adapted mesh).

In the following, two different types of problems are considered and solved numerically. In the first one, the detailed kinematics of the N-particle system is marched in time using an adaptive two-step Adam–Bashforth method, where the time step is determined on the smallest-separation between the particles. To prevent particle overlap, a soft repulsive potential is introduced akin to that leading to the corrective relative velocity for the two-particle system. A relatively coarse and uniform mesh (512 nodes per sphere) is used to compute the velocities of the particles, for which computing the complete set of trajectories for 6 particles initially distributed within a circle of radius R_max = 10 requires about 12 hours using Matlab on a desktop computer. This computational cost is prohibitive for large number of simulations (see Section 4.3), hence providing a clear motivation for a reduced-order model of the clustering dynamics. These accurate BEM simulations nevertheless provide the required reference for validation of such model.

In the second class of problems, regularized BEM are used to obtain the translation and rotation velocities of a cluster of phoretic particles (i.e. frozen relative positions). Numerically, the particle assembly is considered as rigid, and the resulting velocity of these stable clusters is computed with inter-particle separation, d_c = 0.01 and d_c = 0.02; a Taylor series expansion of the global swimming velocity of the cluster in terms of contact distance around d_c = 0 is then used to extrapolate the true self-propulsion velocities of clusters when particles are in actual contact. To achieve sufficient accuracy on the global motion of the cluster (in particular resolving properly the flow field around the particles in the lubrication regions), a finer mesh is used with 1538 nodes and refinement near the regions of contact as described previously. Computation of the cluster velocity for N = 13 yet requires ≈28 GB of allocated memory, but only two computations are performed for each cluster shape.

4.2 Clustering vs. self-propulsion

The numerical methods outlined in the previous section are now used to compute the dynamics of N identical spherical particles of unit radius, with chemical properties resulting in a mutual phoretic attraction i.e. AM = −1. For simplicity, their initial arrangements (and hence dynamics) are restricted to two dimensions, the solute diffusion and fluid motion remaining fully three-dimensional. The present approach and formalism could nevertheless be applied directly to 3D motion. Particles are initially arranged randomly in a plane; under the effect of phoretic attractions, complex non-equilibrium clustering dynamics are observed leading to the formation of a stable rigid planar assembly, held together by the balance between phoretic attraction and short-ranged repulsive forces. These clusters are stable in the classical sense that any slight two-dimensional perturbation would return them back to their original configuration. Using optimal packing arguments (attracting phoretic particles tend to minimize their relative distance), particles in such stable clusters are found to arrange on a regular hexagonal lattice where the particles have no relative degree of freedom left i.e. each particle cannot move without separating from its adjacent particle in the lattice.

An example of such dynamics is shown on Fig. 5. Initially, the asymmetric concentration field generated around each particle by its neighbours drives its motion. The general direction of each particle tends to cluster them around the geometric center of the assembly (referred to in the following simply as “center of mass”, although no inertial process is involved here). Initially and up until the particles form a rigid cluster, the velocity of the center of mass of the particles is noticeably at least one order of magnitude smaller than the velocity of individual particles (see Fig. 5): the center of mass of the system remains effectively stationary during that phase. This is somewhat unsurprising given the opposite (attractive) velocities induced by the particles on each other as observed in the 2-sphere case.

Fig. 5 (top) Clustering dynamics of 5 isotropic particles (A = 1, M = −1). The velocity of each particle is shown by black arrows and the color indicates the concentration field. The red cross indicates the position of the center of mass. (bottom) Ratio of the magnitude of the center of mass velocity to the average individual velocity magnitude of the 5 particles. A stable cluster is formed when there is no relative motion between the center of mass and individual particles (green shade). The vertical dashed lines correspond to the three snapshots above.

Yet, once a stable cluster is formed and particles experience no relative motion, the global velocity is interestingly non-zero (albeit small): particles swim as a cluster. This collective self-propulsion arises from the geometric asymmetry of the formed cluster which generates an asymmetric concentration field and phoretic forcing on the assembly, a phenomenon that was recently characterized for rigid systems of homogeneous properties and various shapes.^36,37

This transition clearly decomposes the dynamics of N phoretic particles into two different regimes, that differ fundamentally in the relative magnitude of the mean and relative velocities of the particles. In the clustering phase (Phase I), particles are individually force-free hydrodynamically. In the self-propulsion phase (Phase II), they are rigidly-bound by the balance of phoretic attraction and repulsive internal forces that prevent their overlap, and the hydrodynamic force on each particle is now non-zero. This profoundly modifies each particle's hydrodynamic signature. The results of the numerical simulation of the full Laplace and Stokes problems show that (i) Phase I is characterized by a global velocity (i.e. center of mass velocity) that is negligible in front of individual motion and (ii) this global velocity is significantly enhanced in Phase (II) when the particles are rigidly-bound.

These two features are confirmed by computing the asymptotic expansion of the particles' velocity during each phase using the method of reflections⁵⁹ (see Appendix B). When the particles are hydrodynamically force-free individually (i.e. not in contact, Phase I), the translation velocity of particle j is obtained as

with d_jk = |R_j − R_k| the interparticle distance and e_jk = (R_k − R_j)/d_jk the unit vector directed from particle j to particle k. The velocity of the center of mass is then

Considering now the final clustered configuration (Phase II), the translation and rotation velocities of the cluster, U_cluster and Ω_cluster are obtained as

with .

These results are valid up to O(ε⁶) corrections with ε ∼ a/d the typical radius-to-distance ratio. At leading order in ε, these results indeed confirm that

(i) in the clustering phase (Phase I), particles exhibit O(ε²) individual velocities, while their center of mass velocity is O(ε⁵),

(ii) in the self-propulsion phase (Phase II), the global cluster velocity is O(ε³) and there is no relative motion.

The analysis presented in Appendix B also demonstrates the physical origin of this increased global motion as a consequence of the non-uniform hydrodynamic resistance experienced by the different particles within the cluster, illustrated by the right-hand-side of eqn (34). At leading order, the phoretic forcing is pairwise and symmetric and therefore averages to zero. Differences in hydrodynamic resistance between two particles however introduces an imbalance that must be compensated by a global translation and/or rotation.

This decomposition of the dynamics of the particles into two successive and physically distinct phases is of further importance. Self-propulsion is essentially restricted to Phase II where the particles are clustered together. But, since the self-propulsion of the particles depend on the geometry of the cluster formed as an outcome of the clustering phase (Phase I), the self-propulsion velocity of the assembly is indirectly, but entirely determined by the detailed kinematic process of phoretic clustering.

4.3 A reduced-order model of the clustering dynamics: pairwise interactions

Characterizing the collective self-propulsion properties of isotropic particles therefore critically relies on a good understanding and modelling of the clustering phase. The clustering process conditions the final shape and therefore collective propulsion properties. Besides the maximum velocity achievable for N particles, prominent self-propulsion properties such as the mean and most probable velocities require determining the probability of formation of a given cluster configuration through the phoretic attraction of N particles that are initially randomly-distributed. The approach chosen here is a Monte-Carlo framework where, for each value of N, a large number of clustering simulations are run starting from random initial particle arrangements, and the probability to obtain a given cluster shape (indexed by q) is computed from the number of runs leading to that particular cluster shape. Due to their attractive nature, phoretic interactions of the particles tend to maximize their packing in their clustered configuration. Restricting here our analysis to two-dimensional clusters, particles are arranged on a regular hexagonal lattice in their final configuration. For N ≤ 5 particles, a single stable shape is obtained (see Fig. 4), while for N > 5, the number of final distinct cluster configurations is finite but increases exponentially with N; for example: N = 6, N = 8 and N = 10 lead to 3, 9 and 35 distinct configurations, respectively.

Boundary element methods are well suited to compute the detailed dynamics accurately, but are prohibitively expensive for running thousands of simulations of the full temporal dynamics as needed for obtaining the probabilities of different cluster shapes starting from random initial positions (see Section 4.1). Motivated by the distinct features of the clustering and propelling phases identified in the previous section, our approach is therefore to split the problem into two distinct parts: (i) determine the distribution of probabilities for the different cluster shapes using a reduced-order model of the clustering dynamics, and (ii) compute the exact propulsion velocity of each final cluster using regularized BEM.

The reduced-order model of the clustering phase required for the first part must be sufficiently accurate, both in the far-field and near-field limits, so that the final configuration is the same as that predicted by BEM, yet be sufficiently inexpensive computationally to be able to run a large number of simulations for each N. At leading order (i.e. far-field approximation), the relative velocity of the particles is determined as the superposition of symmetric interactions between pairs of particles. However, this method fails to provide accurate solutions when the separations between particles are small (d_c < a). The method of reflections detailed in Appendix B, which provides iterative approximations of increasing accuracy, is an appealing alternative to full simulations as it can capture the multi-particle dynamics, but it is fundamentally restricted to distances greater than the size of the particle (a/d < 1) for numbers of reflections small enough for practical implementation.

Turning back to the two-particle case, an exact solution for the Laplace and Stokes problems was obtained in Section 3 which includes a full description of confinement and lubrication effects. We exploit this solution here to account for far- and near-field dynamics properly, by superimposing the pairwise phoretic drift velocities between the particles whilst retaining the full exact solution in eqn (22)U_a(d_jk) (including the repulsive interaction that prevents particle overlap). Note that a fundamental restriction of this superposition assumption is its inability to predict any collective propulsion: the velocity of the center of mass of the arrangement is identically zero by definition. Yet, we observed in the full dynamics that such center of mass motion is negligible during the clustering phase, so that this constraint has little physical implication in practice. It should further be emphasized that this superposition of pairwise solutions is not derived from first principles, but rather on its practical ability to capture correctly the leading order interactions both for near- and far-field limits, and to match the BEM predictions with good accuracy.

The accuracy of this reduced-order interaction model is therefore checked in detail against BEM simulations in for N = 3 and N = 6 particles. Fig. 6 illustrates the clustering dynamics of N = 3 particles starting from arbitrary initial conditions and quantitatively compares the rate of change in relative distance at each time as predicted by the reduced-order model to that obtained using full BEM simulations. The model is observed to predict the clustering velocities with excellent accuracy except for the final stages of clustering where particles are already in contact but for two of them, and the final shape is already fully determined. A second validation of the ability of the reduced-order model to predict the final configuration correctly is proposed in Fig. 7 where the evolution in time of three different initial arrangements of six particles (N = 6) are plotted.

Fig. 6 Top: From right to left, snapshots of the phoretic clustering in a system of N = 3 particles obtained using BEM. For visualization purposes, (a) is zoomed out; the coloured lines represent the distance between centers of different pairs of particles. Bottom: The rate of change of separation (ḋ_c) between particles due to phoretic attraction predicted by the reduced model (dashed lines) and BEM (solid lines). The positions of the particles at each instant are set to that of the full BEM simulations. To avoid spurious repulsive velocities arising from slight differences in the equilibrium distance in the reduced model and BEM results and the stiffness of the repulsive model, the repulsive velocities in the reduced-order model are removed once two particles are in contact. The stages of clustering (a) to (e) are also shown (dashed-dotted lines).

Fig. 7 Comparison of the positions of the particles predicted by BEM (in black) and reduced model (in red) at various instances. The three different initial distributions of the particles, which differ only in the position of a single particle (in blue), lead to three different cluster shape using BEM, all of which are correctly captured by the reduced-order model.

These results confirm that this model is sufficient to capture the final shape of the cluster (which is the only piece of information that it will be used for): slight discrepancies in absolute positions of the particles arise only in the final stages of clustering, and do not affect their relative arrangement (i.e. they only correspond to a slight translation or rotation of the whole assembly). Hence, the reduced-order model is chosen in the following as a suitable substitute for BEM for modelling the clustering phase (Phase I) as it presents a good compromise between the accuracy of its prediction of the final shape and efficient computations (a few seconds are necessary for a single run, compared to tens of hours with BEM).

In practice, the velocity of particle j is computed in this reduced-order model as

with U_a(d) given in eqn (22). One should note that this reduced model is formally equivalent to the motion of N particles down the steepest gradient of an interaction potential :with

It should be stressed here that this formal equivalence only holds rigorously for the reduced-order model, and not for the full problem (e.g. in the latter, phoretic interactions lead to a net global motion). Within the reduced-order model framework, particles arrange to minimize which is equal to zero for infinitely-distant particles. As such, can be understood as a measure of the cohesive interactions within and stability of a particular cluster shape. Noting its value for cluster shape q, a Boltzmann distribution would therefore be expected if the cluster formation was purely stochastic and driven solely by the stability of the final shape, with the probability to obtain cluster q equal to

and σ² characterizes the background noise in the system: σ = 0 leads to for the most stable cluster (i.e. that with minimum ), and for all others, while σ = ∞ results in all cluster shapes having the same probability.

5 Phoretic clustering and self-propulsion

Using the reduced-order clustering model described in the previous section (Phase I), a large number of simulations are performed, with random initial spatial arrangements of particles in order to obtain the probability of formation of the different cluster configurations through phoretic attractions. In a second step, the self-propulsion velocity of each cluster (Phase II) is computed accurately using BEM and hence, the self-propulsion statistics are obtained for a given system of N particles (i.e. maximum, mean and most probable velocities).

5.1 Probabilities of formation of stable rigid clusters

To study the complete evolution of clustering of an N-particle system from zero interaction potential to a final cluster potential , the particles have to ideally begin from an infinite separation. For the practical purpose of simulations, the particles are initially distributed randomly within a large disk of radius R_max. It is ensured that R_max is sufficiently large that it does not significantly affect the probability statistics (Fig. 8). For N = 6, which is the smallest value of N for which multiple stable configurations are obtained, probability values converge for R_max ≳ 20 (see Fig. 8), which corresponds to a density (area fraction of particles in the clustering plane) of 1.5%. For all N ≤ 12, it is observed that this area fraction of 1.5% is sufficient to give accurate probability statistics.

Fig. 8 Stable 6-particle and 7-particle clusters and their probabilities obtained from 2000 independent trials. Bar graph shows probabilities when particles begin at different random locations within a circle of radius R_max for each trial (histogram shown for R_max = 20 for 6-particle and R_max = 22 for 7-particle system respectively). Measurements are made for different R_max to show its influence. Range bars represent 2 standard deviation limits of probabilities computed from 50 such events.

The number of independent trials needed to determine the probabilities of configurations accurately increases with N. We found that for N = 6 and N = 7, 2000 distinct runs are sufficient and the resulting probability distributions and accuracy are shown in Fig. 8. Clusters are labelled by their effective interaction potential , i.e. the most stable cluster is on the left.

Interestingly, Fig. 8 shows that the most stable cluster (i.e. that with least effective potential in the reduced-order model) does not have the highest probability as one would expect in a stochastic process. The phoretic clustering described by eqn (37) (i.e. down-gradient of the interaction potential) does not lead to an absolute but to a local minimum, as it depends on the detailed route followed in the configuration phase space during the clustering process. Less “stable” cluster shapes minimize the interaction potential only locally but may be wider attractors in the (2N − 3)-dimensional phase space that characterizes the clustering motion. Because of the intimate link of cluster shape and velocity, this is expected to hold profound consequences on the collective self-propulsion properties of the N particles.

Once formed, a cluster cannot transition to another configuration without additional forcing (i.e. “energy” input). Even though the differences in the final potential between configurations are relatively small in comparison to and (i.e. the change in potential from the initially-dispersed configuration to the final clustered shape), switching from one configuration to another requires overcoming a much larger potential barrier . This is illustrated in Fig. 9 in the case of N = 6 particles: changing the position of a single particle around the rest of the cluster allows the arrangement to describe all three possible stable configurations (identified as minima in effective potential ).

Fig. 9 Evolution of the assembly potential while rolling a single particle around the cluster from an initial stable configuration A of potential in order to sample all configurations. Clusters A, B and C represent local minima of the effective interaction potential . represents the minimum gain in effective potential required to transition from cluster A to cluster B, while the actual potential difference between clusters is just . is the change in effective potential between the initial stage and final clustered configuration. The results for 180° ≤ θ ≤ 360° are obtained by symmetry.

While the difference in effective potential of the three stable configurations is small ({p, q} ⊂ {A, B, C}), transition from one configuration to another requires reaching intermediate stages representing a potential barrier typically an order of magnitude larger, . Note also that the difference in effective potential between the pre-clustering state (ideally, ) and final configuration is yet an order of magnitude larger .

5.2 Propulsion velocities of rigid clusters

Geometrically-asymmetric phoretic systems self-propel with a global translation and/or rotation velocity determined only by the asymmetry of their geometric shape;^37,60 in particular, this velocity is independent of the particles' size. As expected from symmetry arguments, configurations with a single mirror plane of symmetry cannot undergo any rotation, while configurations which are antisymmetric about a plane (left unchanged by a 180° rotation) do not translate. Moreover, configurations with rotational symmetry neither translate nor rotate (see Fig. 8 for example).

As explained, in this clustered arrangement, the relative positions of the particles are fixed in a regular hexagonal lattice. For each N-particle system, the regularised Boundary Element Method detailed in Section 4.1 is used to compute the exact translational and angular velocities of propulsion (Phase II) of all possible cluster configurations. The fastest propelling and rotating clusters, as well as the most probable ones are shown in Fig. 10 for 9 ≤ N ≤ 12. By defining the normalised effective potential for shape q for fixed N as , we observe that the fastest propelling (translation or rotation) clusters are some of the least stable clusters with respect to the normalised effective potential . This is the result of their large geometric eccentricity, responsible for their larger propulsion velocity through the larger concentration gradients at their surface it creates. In contrast, the most probable clusters are characterized by their more compact and roughly symmetric arrangement around their geometric center (see Fig. 10).

Fig. 10 Fastest translating (top), fastest rotating (center) and most probable clusters (bottom) for N = 9 to 12. The relative surface concentration is shown in color for each cluster as well as its propulsion, stability and probability characteristics.

5.3 Clustering-induced self-propulsion properties

Combining the results of the previous sections, the statistical properties of the collective self-propulsion of N isotropic particles resulting from their phoretic clustering are obtained here. Fig. 11 shows the mean and maximum velocities of N-particle clusters. In particular, for a given N, noting the probability and U_q the velocity magnitude (or equivalently, rotational velocity magnitude Ω_q) of cluster shape q, the mean velocity is defined aswith the sum carried out on all possible cluster shapes.

Fig. 11 Magnitudes of maximum (dashed) and mean (solid) translational and rotational velocities of clusters obtained using BEM. The grey shading represents the mean velocity obtained through an equilibrium distribution, eqn (38) with the noise level σ², increasing from zero to infinity as the intensity of shading gets darker.

Because of the symmetry of the only existing final cluster, no self-propulsion is observed for N ≤ 4. Note that for all N > 5, at least one non-propelling cluster (having rotational symmetry) is found; and therefore the minimum velocity for all N > 5 is strictly zero. The maximum propulsion velocity is observed to increase steadily with N as the larger number of particles allows for more eccentric shapes and larger phoretic forcing. However, the mean velocity U_mean(N) is observed to saturate for larger values of N. Insight on this result as well as other statistical properties of importance such as the most probable velocity and the variance in propulsion velocities of large N-particle clusters is obtained by studying the probability distribution of the propulsion velocities (Fig. 12).

Fig. 12 Probability distribution of (top) translational velocities U and (bottom) rotational velocities Ω of N-particle clusters. The distribution mean is indicated by a dashed red line for each N.

The probability distributions of translational velocities of clusters with N ≥ 10 show a prominent peak near the mean value of the distribution (≈2 × 10⁻³) indicating that clusters with this velocity are also the most likely to form. The most probable angular velocity, however, is lower than the mean value (≈2.5 × 10⁻⁵). As the maximum velocity (translation and rotation) increases with N, the graph spreads indicating an increase in variance of velocities between the various configurations. However, the fastest propelling configurations, which are also some of the most eccentric and elongated ones, have very low probability so that such clusters are only seldom observed for large N, and hence their contribution to the mean velocity is minimal. For N ≥ 10, most of the contribution to the mean properties is brought by clusters with an intermediate velocity which remains relatively fixed with N, leading to the saturation in the mean velocity as N grows. This can be qualitatively understood from the most dominant features of the most probable clusters as depicted in Fig. 10: those clusters display more compact shapes with an asymmetry arising only from a small number of particles.

It was already emphasized above that the probability of formation of a given shape in the phoretic clustering process does not correlate with its effective potential for the reduced-order model, as would be expected in a classical system at thermodynamic equilibrium, for which the probability of formation of various configurations follows a Boltzmann distribution, eqn (38). Not unexpectedly, this has consequences for the collective propulsion properties: Fig. 11 shows that the mean velocity resulting from phoretic interactions is systematically larger than it would be under the sole constraint of minimizing the effective interaction potential, regardless of the importance of noise in the process.

6 Effect of noise

The previous sections focused on the clustering-induced collective dynamics of phoretic particles, in purely deterministic systems: for a given particles' arrangement at time t, the evolution of the particles' position in time is fully determined. Yet, in practice, all such systems which focus on microswimmers are subject to a variety of stochastic processes (including thermal noise) which continuously alter their motion. In the presence of background noise caused by temperature of the surrounding fluid, passive microscopic particles (or isolated isotropic phoretic swimmers) undergo Brownian motion characterized by zero mean displacement. Active particles however have a net displacement with continuous reorientation of the direction of propulsion, thus exhibiting diffusive behaviour in long time scales.^61,62 Commonly employed minimal models for these active Brownian particles describe the essential dynamics involving overdamped motion (Langevin dynamics) as well as their thermal reorientation.^26,62

Equivalently, in the case of isotropic phoretic particles, we describe the clustering process by its reduced-order model (see Section 4.3) rather than using the full BEM simulations because of its versatility and reduced computational cost, as well as to maintain consistency with the rest of the manuscript. Hence, the deterministic dynamics of the particles corresponds to a collective minimization of the effective interaction potential as defined in eqn (37). The purpose of this section is therefore to provide some insight on the robustness of these deterministic results, i.e. how the results of self-propulsion statistics obtained in the previous sections are modified in the presence of background noise on the kinematics of individual particles.

In the absence of inertia, the evolution of the position R^(b)_j(t) of particle j under the effect of background noise is given by the overdamped Langevin equation:

where its deterministic velocity, U_j(t) is given by eqn (37), and ξ_j(t) is a external Gaussian white noise, with zero mean and a variance, σ²I = 〈ξ_i(t)ξ_j(t′)〉 = 2Dδ(t − t′)δ_ijI, where D is the diffusivity of each particle in the fluid. Thus, the instantaneous displacement of particle j at any time t is

dR^(b)_j(t) = U_j(t)dt + dW_j, (41)

where W_j is a Weiner process with zero mean and variance σ²tδ_ijI. Discretizing eqn (41) by using Euler–Mayurama method gives

R^(b)_j(t + Δt) = R^(b)_j(t) + U_j(t)Δt + ΔW_j, (42)

ΔW_j is also a zero mean Weiner process with variance σ²Δtδ_ijI. Eqn (42) is solved using an adaptive time stepping method.

A system of N = 6 particles is the smallest system that exhibits multiple clustered configurations, and for which the introduction of noise is expected to potentially introduce significant modification in the collective dynamics; we shall henceforth consider N = 6 as an example to illustrate the effect of noise on the clustering statistics and resulting propulsion properties. A direct result of the component of randomness in position of particle in eqn (42), is the formation of multiple configurations from the same initial conditions of the system. For each set of initial conditions (typically a few hundreds), eqn (42) is used to run about 100 simulations. If the strength of noise satisfies , background fluctuations are sufficient to break the formed clusters and redistribute the particles far apart. Hence, we restrict ourselves to the case where σ² is much smaller than the absolute cluster potential, so that phoretic effects are still dominant and clustering occurs. In the following, two different behaviours are observed under the effect of noise, depending on its relative magnitude (σ²) and the effective potential barrier from configurations p to q (see Fig. 9).

6.1 Low noise

Low noise conditions are characterised by , which is not sufficient to change the configuration of a stable cluster once formed. However, it can influence the clustering dynamics by rearranging the particles as the system moves down the steepest gradient in interaction potential, as seen in Fig. 13.

Fig. 13 Formation of different clusters from the same initial positions of particles in the presence of noise of strength σ² = 0.01. Variation in potential of system with time during clustering for the cases shown.

Its influence is particularly important when particles are far apart, i.e. when differences in are small and of order σ². For this reason, one may wrongly presume that noise would change the probability statistics for the formation of cluster configurations. Instead, the probability statistics remain identical to that in the absence of noise (see Fig. 14).

Fig. 14 Comparison between probabilities of formation of 6 particle cluster configurations in the absence of noise and with low noise shows no significant difference in the statistics. The clustering was performed with R_max = 16 in both cases. In the case of low noise (σ² = 0.01), probabilities are computed from 200 random initial arrangements of particles with 100 trials per arrangement.

This result can be explained as follows: although a low background noise is expected to enable the system to explore neighboring routes in the configuration space to minimize the system's potential, , during the clustering phase, it is only effective in the initial stages where system has low ; in the later stages, the magnitude of the noise becomes too small in front of the deterministic velocity (determined by ) to significantly alter the particles' trajectories. However, a lower potential during the initial stages does not ensure a low potential of the cluster formed eventually (as observed in Fig. 13). Thus, the noise just effectively redistributes the initial arrangement of the particles. But this randomness is already taken into account in the probability of formation of different clusters by considering a large number of random initial positions. As a result the collective propulsion statistics in Fig. 11 remain unmodified in low-noise conditions.

6.2 High noise

In the presence of sufficiently large noise, , a cluster, once formed, continuously transitions from one configuration to another without dislocating fully since the noise intensity remains much smaller than . This situation is in stark contrast with the low-noise dynamics for two main reasons: (i) it is not possible anymore to define a fixed cluster shape to the arrangement of the particles which continuously evolves in time and transitions between all the available configurations, and (ii) the clustering dynamics does not play a significant role anymore. Indeed, since the particles' arrangement can be reconfigured, which cluster shape was reached in the first place is not relevant. Instead, how much time the particles spend in a particular configuration is now the significant information. This can be characterized as the probability to obtain a given configuration at any time, which is the result of the equilibrium between the background noise and the phoretic attraction. Not surprisingly, this probability now follows a Boltzmann distribution (Fig. 15). However, since multiple intermediate configurations (with the same potential) exist, the probability of a particular stable configuration cannot be defined.

Fig. 15 Formation of different cluster configurations by rearrangement of particles over time due to high noise, σ² = 0.1. Also shown is the probability distribution of various cluster configurations which follow a Boltzmann distribution. Since the average separation between particles are larger due to the continuous high noise, takes larger values.

7 Conclusions

The results presented in this work therefore establish a fourth route to self-propulsion of active colloids. Unlike previously-identified strategies which relied either on an asymmetric design of the particle or the non-linear convective transport of solutes by the phoretic flows, self-propulsion of individually non-motile yet active particles is achieved here by non-symmetric interactions between multiple particles. Under the effect of attractive phoretic attractions and steric constraints, particles form geometrically-asymmetric clusters that are able to maintain the asymmetric concentration fields required for propulsion. The self-propulsion velocities are much lower (by at least an order of magnitude) than typical relative velocities between particles during particle clustering, but lead nonetheless to a net migration of the particles.

Even though the governing dynamics of the particles (in the absence of external noise) are completely deterministic, multiple particle arrangements can be reached depending on the detailed dynamics of the cluster formation and balance between phoretic attraction, hydrodynamic interactions and steric repulsion. Slightly different initial positions of the different particles may therefore lead to fundamentally different propulsion characteristics, thereby introducing an inherent stochasticity in the system. Focusing for simplicity on two-dimensional arrangements of the particles (the hydrodynamics and solute diffusion are three dimensional), the probability of formation of the different cluster shapes was determined using Monte-Carlo simulations of the phoretic clustering for N particles with initially-random positions. These simulations were performed using a reduced-order model of the clustering process obtained by superimposing the velocity induced on a given particle by each of its neighbors as in a two-particle system, for which an analytical solution was first obtained.

Using these probabilities and the velocity of the different clusters computed numerically using a regularized boundary element method, the statistics of the collective self-propulsion were obtained, namely the mean, maximum and most probable velocities, for both translation and rotation. The maximum attainable velocity was found to increase with N: for larger clusters, a greater degree of asymmetry can be achieved leading to larger velocities. The mean propulsion velocity on the other hand was found to saturate for increasing N: the most eccentric clusters with large velocity also become less and less probable as N increases, while the most probable clusters exhibit a more compact geometry with velocities close to the mean value.

The particles are continuously in a state of non-equilibrium. The lack of external noise to relax the system creates probability statistics of various configurations which are strikingly different from equilibrium statistics. This difference directly affects the velocity statistics and we observe a larger mean velocity compared to that expected from a purely equilibrium process. Yet thermal noise and Brownian motion can become important for smaller particles, and their impact on the present result was tested by introducing a Brownian component in the motion of the active particles in addition to their deterministic phoretic clustering. Two different situations are observed depending on the magnitude of the noise with respect to the phoretic forces maintaining the cohesion of the formed clusters. The cluster probability and velocity statistics are found remarkably robust to low noise amplitude (i.e. thermal fluctuations that are unable to break or reconfigure a given cluster): such noise levels essentially redistribute the already random position of the particles in their initial state (before clustering). In contrast, higher levels of thermal noise can lead to cluster break-up or reconfiguration. In the latter, a given cluster undergoes continuous transition between different configurations, and the probability to find a given shape follows an equilibrium probability distribution. In that case, the group of N particles continuously reconfigure leading to changes in its translational and rotational velocity. This provides it with a self-propelled motion which is reminiscent in some regards to the run-and-tumble behaviour of bacteria where self-propulsion in a given configuration is followed by a geometric reorganization that modifies its swimming direction and velocity magnitude.³⁴ Although the process is quite different here, in particular because the relative length of a run is significantly smaller than for bacteria, it is expected to enhance transport by diffusion.

The central goal of the present work is to characterize the emergence of the self-propulsion from phoretic and hydrodynamic interactions of individually-non-motile isotropic particles. We demonstrated this property for the simpler configuration of two-dimensional motion of the particles. Yet, the different tools developed both analytically and numerically directly apply to three-dimensional and the present results could easily be extended to generic 3D clusters and associated propulsion properties, the main difficulty residing in the classification of the different cluster shape. We finally remark that the formalism developed in this manuscript, both numerical and analytical, could also be extended to address the collective dynamics of mixtures of particles of different chemical properties (e.g. mixture of active and inert, or isotropic and Janus particles^32,35).

Conflicts of interest

There are no conflicts to declare.

Appendix

A Effect of a finite regularization parameter

We present briefly here the derivation of eqn (28) which provides the leading order correction introduced by a finite regularization parameter ε in the regularized boundary element method framework. The general goal in BEM is to express the following integralwithin terms of the local value of c(x₀) and its gradients, when x₀ is on the surface of the particle. Here, V_f is the fluid volume outside the particle. To perform this computation, we consider a local reference frame centred on x₀ and such that n(x)₀ = e_z. The fluid domain outside the particle is locally defined as z + (κ_xx² + κ_yy²)/2 ≥ 0 with κ_x and κ_y are the two principal curvatures (this implicitly assumes that |κ|ε ≪ 1 with κ the mean curvature).

First, c(x) is expanded in Taylor series around x₀:

and using a local spherical polar coordinate system with μ = cos θ and μ_c(r,ϕ) = −(κ_x cos² ϕ + κ_y sin² ϕ)r/2 the parametric description of the local particle surfaceThe development above includes the first two leading orders in ε: the O(ε⁰) contribution is the classical result for a singular boundary integral methods and the O(ε) accounts for the finite size of the regularization in comparison with (i) the local curvature of the surface and (ii) the typical length scale for the variations of c.

B Method of reflections

The method of reflections is a classical technique to construct an asymptotic solution of Laplace or Stokes problem around several or many bodies using an iterative approach.⁵⁹ This asymptotic development is made in terms of ε = 1/d, i.e. the ratio of the particles' radius to typical distance among them. Taking the concentration problem as an example, the solution c⁰_j is first obtained for each particle j in the absence of all the other particles. These solutions are then corrected at each stage (i.e. each reflection) near each particle k ≠ j to cancel the extra contribution to the boundary conditions introduced at the surface of this particle by the previous reflections.

B.1 Diffusion problem. The Laplace problem for the concentration around all the particles of unit radii can be formulated in full as follows

∇²c = 0, (47)

n_j·∇c = −A for r_j = 1, (48)

c → 0 for r → ∞. (49)

As for the main text of the manuscript, r_j = r − R_j is the position relative to particle j, n_j its outward-pointing normal, d_jk = |R_j − R_k| the relative distance of particles j and k, and e_jk = (R_k − R_j)/d_jk is the unit vector point from particle j to particle k.

The solution of this problem for an isolated particle l is simply c⁰_l = A/r_l. Expanding this solution near particle k ≠ l (i.e. for small r_k),

where L_p(μ) is the Legendre polynomial of degree p.

The first reflection solution c¹_k therefore satisfies ∇²c¹_k = 0, decays at infinity, and at the surface of particle k

The unique solution is obtained asNote that this expansion includes no source term (p = 0) since c⁰_l satisfies ∇²c⁰_l = 0 in the vicinity of particle k.

The leading order correction introduced by the second reflection on particle j is a response to the uniform gradient G¹_j, created near particle j by the source dipole contribution in each c¹_k (k ≠ j) in the above equation, namely

and the leading order contribution to the second reflection writesUsing these results, the concentration c^s_j and slip velocity ũ_j = M(I − n_jn_j)·∇c^s_j at the surface of particle j using two reflections, are given by

c^s_j = C^j₀ + C^j₁·n_j + C^j₂ : (n_jn_j) + ⋯, (55)

ũ_j = C^j₁ + 2C^j₂·n_j + ⋯, (56)

with

B.2 Hydrodynamic problem for free particles. We now turn to the iterative solution of the Stokes problem for individually force-free particles, namely

η∇²u = ∇p, ∇·u = 0, u(∞) = 0 (59)

∀j, u = U_j + Ω_j × r_j + ũ_j on r_j = 1, (60)

Following the general framework of the method of reflections, the swimming problem for a single isolated particle k is solved for the velocity, rotation rate and stresslet of this particle as

S⁰_k = −10π〈n_kũ_k + ũ_kn_k〉 = −8πMC^k₂, (64)

The leading order flow field generated by particle j is then

The development above retains only the source and force dipoles as fundamental singularities. Force quadrupoles, which have the same far-field decay as that of the source dipole (1/r³), are nevertheless neglected because their intensity is determined by C^j₃ which is of subdominant order. The result of the first reflection is the swimming and rotation velocities of force- and torque-free particles in an ambient flow field, and are determined directly using Faxen's laws for spherical particles. As a result,Using the results for the Laplace problem (Section B.1), the velocity of force-free particle j is obtained asand the velocity of the center of mass, , is

For free (i.e. non-touching) particles, individual particles' velocities scale as 1/d² while their center of mass moves with an 1/d⁵ velocity.

B.3 Hydrodynamic problem for rigid clusters. When the particles are rigidly-bound in a cluster, they are not individually force-free, but, the total force and torque on the cluster are zero. Assuming that the non-hydrodynamic interaction forces between the particles are central, the net torque on particle j around its center of mass must vanish. Hence,

Within the framework of the method of reflections, solving the hydrodynamic problem for particle k, leads to

Particle k is not force-free anymore, and its flow signature must now also include a contribution from an additional Stokeslet and source dipole associated with F_k:From Faxen's law, the correction to the velocity of particle j introduced by the first reflection is thereforeUsing the result for the concentration problem and the rigid body kinematics, U_j = U_CM + Ω_CM × R_j,withThe translation and rotation velocity of the cluster are then obtained by imposing that the cluster as a whole is force- and torque-free, i.e.. The resulting leading order velocity of the center of mass, U_CM, now scales as 1/d³.

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 714027 to SM) and the Engineering and Physical Sciences Research Council of the United Kingdom (grant EP/R041555/1 to TDMJ). The authors would like to acknowledge valuable discussions with Eric Lauga.

Notes and references

1. E. M. Purcell, Am. J. Phys., 1977, 45, 3.
2. J. Happel and H. Brenner, Low Reynolds number hydrodynamics, Springer, 1965.
3. C. Brennen and H. Winnet, Annu. Rev. Fluid Mech., 1977, 9, 339–398.
4. E. Lauga and T. R. Powers, Rep. Prog. Phys., 2009, 72, 096601.
5. E. Lauga, Annu. Rev. Fluid Mech., 2016, 48, 105–130.
6. A. Najafi and R. Golestanian, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 69, 062901.
7. G. Grosjean, M. Hubert, G. Lagubeau and N. Vandewalle, Phys. Rev. E, 2016, 94, 021101.
8. A. Ghosh and P. Fisher, Nano Lett., 2009, 9, 3.
9. R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone and J. Bibette, Nature, 2005, 437, 3.
10. J. L. Moran and J. D. Posner, Annu. Rev. Fluid Mech., 2017, 49, 511–540.
11. F. Jülicher and J. Prost, Eur. Phys. J. E: Soft Matter Biol. Phys., 2009, 29, 27–36.
12. J. F. Brady, J. Fluid Mech., 2011, 667, 216–259.
13. V. Yadav, W. Duan, P. J. Butler and A. Sen, Annu. Rev. Biophys., 2015, 44, 77–100.
14. W. Duan, W. Wang, S. Das, V. Yadav, T. E. Mallouk and A. Sen, Annu. Rev. Anal. Chem., 2015, 8, 311–333.
15. J. L. Anderson, Annu. Rev. Fluid Mech., 1989, 21, 61–99.
16. R. Golestanian, T. B. Liverpool and A. Ajdari, New J. Phys., 2007, 9, 126.
17. G. Volpe, I. Buttinoni, D. Vogt, H.-J. Kümmerer and C. Bechinger, Soft Matter, 2011, 7, 8810–8815.
18. I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger and T. Speck, Phys. Rev. Lett., 2013, 110, 238301.
19. J. R. Howse, R. A. L. Jones, A. J. Ryan, T. Gough, R. Vafabakhsh and R. Golestanian, Phys. Rev. Lett., 2007, 99, 048102.
20. W. Wang, W. Duan, A. Sen and T. E. Mallouk, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 17744–17749.
21. M. N. Popescu, M. Tasinkevych and S. Dietrich, Europhys. Lett., 2011, 95, 28004.
22. J. Catchmark, S. Subramanian and A. Sen, Small, 2005, 1, 202–206.
23. I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert and L. Bocquet, Phys. Rev. Lett., 2012, 108, 268303.
24. J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine and P. M. Chaikin, Science, 2013, 339, 936–940.
25. W. Wang, W. Duan, S. Ahmed, T. E. Mallouk and A. Sen, Nano Today, 2013, 8, 531–554.
26. A. Zöttl and H. Stark, J. Phys.: Condens. Matter, 2016, 28, 253001.
27. P. Illien, R. Golestanian and A. Sen, Chem. Soc. Rev., 2017, 46, 1.
28. S. Saha, R. Golestanian and S. Ramaswamy, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 062316.
29. J. A. Cohen and R. Golestanian, Phys. Rev. Lett., 2014, 112, 068302.
30. R. Golestanian, Phys. Rev. Lett., 2012, 108, 038303.
31. B. Liebchen, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2017, 118, 268001.
32. R. Soto and R. Golestanian, Phys. Rev. Lett., 2014, 112, 068301.
33. R. Soto and R. Golestanian, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 052304.
34. M. S. Davies Wykes, J. Palacci, T. Adachi, L. Ristroph, X. Zhong, M. D. Ward, J. Zhang and M. J. Shelley, Soft Matter, 2016, 12, 4584–4589.
35. F. Schmidt, B. Liebchen, H. Löwen and G. Volpe, 2018, arXiv:1801.06868v1.
36. S. Shklyaev, J. F. Brady and U. M. Cordova-Figueroa, J. Fluid Mech., 2014, 748, 488–520.
37. S. Michelin and E. Lauga, Eur. Phys. J. E: Soft Matter Biol. Phys., 2015, 38, 7.
38. S. Michelin, E. Lauga and D. Bartolo, Phys. Fluids, 2013, 25, 061701.
39. Z. Izri, M. N. van der Linden, S. Michelin and O. Dauchot, Phys. Rev. Lett., 2014, 113, 248302.
40. A. Dominguez, P. Malgaretti, M. N. Popescu and S. Dietrich, Soft Matter, 2016, 12, 8398–8406.
41. R. Golestanian, T. B. Liverpool and A. Ajdari, Phys. Rev. Lett., 2005, 94, 220801.
42. S. Michelin and E. Lauga, J. Fluid Mech., 2014, 747, 572–604.
43. S. Ebbens, M.-H. Tu, J. R. Howse and R. Golestanian, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 020401(R).
44. N. Sharifi-Mood, J. Koplik and C. Maldarelli, Phys. Fluids, 2013, 25, 012001.
45. Y. Ibrahim, R. Golestanian and T. B. Liverpool, J. Fluid Mech., 2017, 828, 318–352.
46. M. Stimson and G. B. Jeffery, Proc. R. Soc. London, Ser. A, 1926, 111, 757.
47. S. Y. Reigh and R. Kapral, Soft Matter, 2015, 11, 3149.
48. H. A. Stone and A. D. T. Samuel, Phys. Rev. Lett., 1996, 77, 4102.
49. E. Yariv, Phys. Rev. Fluids, 2016, 1, 032101(R).
50. R. J. LeVeque and Z. Li, SIAM J. Sci. Comput., 1997, 18(3), 709–735.
51. G. J. Wagner, N. Moës, W. K. Liu and T. Belytschko, Int. J. Numer. Meth. Eng., 2001, 51, 293–313.
52. Z. Li and M.-C. Lai, J. Comput. Phys., 2001, 171, 822–842.
53. C. Pozrikidis, Boundary Integral and Singularity Methods for Linearized Viscous Flow, Cambridge University Press, 1992.
54. W. E. Uspal, M. N. Popescu, S. Dietrich and M. Tasinkevych, Soft Matter, 2014, 11, 434.
55. T. D. Montenegro-Johnson, S. Michelin and E. Lauga, Eur. Phys. J. E: Soft Matter Biol. Phys., 2015, 38, 139.
56. R. Cortez, SIAM J. Sci. Comput., 2001, 23(4), 1204–1225.
57. R. Cortez, L. Fauci and A. Medovikov, Phys. Fluids, 2005, 17, 031504.
58. D. J. Smith, Proc. R. Soc. London, Ser. A, 2009, 465, 3605–3626.
59. S. Kim and S. J. Karrila, Microhydrodynamics: Principles and Selected Applications, Butterworth-Heinemann, 1991.
60. S. Shklyaev, Europhys. Lett., 2015, 110, 54002.
61. G. Volpe, S. Gigan and G. Volpe, Am. J. Phys., 2014, 82, 659.
62. C. Bechinger, R. D. Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 1.

---

Footnote

† The latter was erroneously omitted in an earlier presentation of the method.⁵⁵

---