Unravelling the role of phoretic and hydrodynamic interactions in active colloidal suspensions

Active ﬂuids comprise a variety of systems composed of elements immersed in a ﬂuid environment which can convert some form of energy into directed motion; as such they are intrinsically out-of-equilibrium in the absence of any external forcing. A fundamental problem in the physics of active matter concerns the understanding of how the characteristics of the autonomous propulsion and agent-agent interactions determine the collective dynamics of the system. We study numerically suspensions of self-propelled diffusiophoretic colloids, in (quasi)-2 d conﬁgurations, accounting for both dynamically resolved solute-mediated phoretic interactions and solvent-mediated hydrodynamic interactions. Our results show that the system displays different scenarios at changing the colloid-solute afﬁnity and it develops a cluster phase in the chemoattractive case. We study the statistics of cluster sizes and cluster morphologies for different magnitudes of colloidal activity. Finally, we provide evidences that hydrodynamics plays a relevant role in the aggregation kinetics and cluster morphology, signiﬁcantly hindering the cluster growth.


Introduction
Collective behaviour is widespread in Nature: fish schools, insects swarms, bacterial colonies, plankton blooms are but a few instances of it. Collective phenomena in Active Matter are characterised by long-ranged correlations and large density fluctuations 1,2 , complex pattern-formation 3 , and non-equilibrium changes of state, such as a flocking 4-6 , clustering 7 , or mobility induced phase separation 8,9 . Answering key questions on how individuals interact and communicate goes even beyond fundamental goal of unravelling the physical mechanisms at the basis of self-organisation in living systems. It will help the design of micro-and nano-scale self-propelled objects [10][11][12][13][14][15][16][17] , with the aim of generating motion in miniaturised devices and developping biomimetic environments [18][19][20][21] . Despite most of these natural and artificial particles displace in a fluid medium, the role played by particle-motion induced hydrodynamic correlations has been essentially overlooked so far. Here we present a numerical study of a suspension of non-Brownian colloids which move responding an advection-diffusion-reaction equation which is integrated numerically by means of a finite difference scheme [32][33][34][35][36] . In Eq. (1), ρ f is the fluid density (assumed to be constant, since the flow regime is close to incompressible, the maximum Mach number being Ma ≈ 10 −2 ), P is the pressure field, ν is the kinematic viscosity. In Eq. (2), D is the solute diffusivity and Q C is a source field that represents the production of solute by the colloids and it is localised at particle surfaces, i.e. it has the form where the sum runs over all colloids (of radius is R), whose positions are X i (t). The function α(x (i) s ) is a prescribed activity profile over the surface (see Fig. 1 and the discussion in the next section). The local sink term −k d C models the degradation of products with rate k d (with associated characteristic screening length d = D k d ). The velocities attained in our simulations are such that the typical particle Reynolds Re = V p R/ν and Péclet Pe = V p R/D numbers are always smaller than 10 −1 (V p is the self-propulsion speed), thus making the advection terms in Eqs. (1)(2) negligible. The fluid is confined along the z-direction by two parallel walls, at which a no-slip boundary condition is imposed on the velocity field and a zero-flux condition applies for the equation for C. Periodic boundary conditions apply in the x and y directions. The boundary conditions at the particles surfaces need a deeper description that will be provided hereafter.

Particle-solute coupling: self-diffusiophoresis
According to the theory of colloidal phoresis 37 , a nonhomogeneous concentration field C around a solid particle causes a pressure imbalance that, in turn, induces a flow within the solidliquid interfacial layer. Such flow appears, on macroscopic length scales, as an effective slip boundary condition for the velocity. The interfacial layer thickness being much smaller than the typical particle size allows a lubrication theory analysis, which leads to the following effective slip velocity for the fluid velocity u at the particle surface Σ: wheren is the normal at the surface in x s and the phoretic mobility µ(x s ), x s ∈ Σ, contains the molecular details of the colloid-solute interaction. For uniform phoretic mobility µ(x s ) ≡ µ 0 , particles gain a net propulsion velocity V p ∼ −µ 0 ∇C, i.e. they perform positive diffusiophoresis for negative µ 0 ("chemoattractive" colloids) and negative diffusiophoresis for positive µ 0 ("chemorepulsive" colloids). Self-diffusiophoresis is realised by letting each particle produce solute C according to a prescribed activity profile α(x s ). Although we may deal, in principle, with arbitrarily patchy 38 (yet axisymmetric) active colloids, in the present study we will specialise to the case of Janus-like particles 39 , that produce solute at a constant rate per unit surface α 0 , homogeneously over only one hemisphere, i.e.: m being the particle characteristic unit vector (see the sketch in figure 1). An isolated free Janus SPC with the activity profile (5) and uniform phoretic mobility µ 0 performs a rectilinear motion with contast speed V p = |µ 0 |α 0 /(4D) 22,23,40 .

Particle dynamics
SPCs are described as solid spheres of radius R, mass M and moment of inertia I = 2 5 MR 2 . The boundary condition for the velocity field at the particle surfaces is implemented via the bounceback-on-links scheme for LB probability densities [41][42][43][44] . Details on the method and on how to introduce in it the concentrationdependent effective slip velocity (4) are given in the Appendix. The Lagrangian dynamics for the position X (i) and velocity V (i) of the centre of mass of the i-th SPC (i = 1, 2, . . . , N), and for its intrinsic orientationm (i) and angular velocity Ω (i) , is described by the following equations of motion: The total force and torque acting on the particle read: where ζ F and ζ T are friction matrices. F h and T h are the force and torque exerted by the fluid on the particle and depend on the global configuration of the velocity field u, thus including also hydrodynamic interactions; whereas, F p and T p depend on the solute field configuration C and mediate phoretic interactions. Equations (6) are solved numerically, time-marching first the positions and orientations vector by means of a forward-Euler scheme, and then the velocities and angular velocities by means of an implicit (backward-Euler) update.

Numerical details
We have performed numerical simulations of suspensions with N = 6400 SPCs of radius R = 2.5 lattice spacings (a value which is relatively small, such to allow simulations of many particle systems, but large enough to keep deviations from the expected physical behaviour, in terms, e.g., of the drag coefficient, below 10% 42,45 ), on lattices of 1024 × 1024 × 24 grid points (≈ 410 R × 410 R × 10 R) at fixed area fraction φ ≈ 0.12. The values, in lattice Boltzmann units (lbu), of the parameters used in the simulations are summarised in table 1. Particles are subjected to a gravity force F b strong enough to prevent them from leaving the bottom wall (the limit fall velocity being five times larger than the self-propulsion speed corresponding to the maximum phoretic mobility considered, i.e. . A soft-core particleparticle repulsion is introduced to prevent overlapping (further details on the treatment of close-to-contact particles are provided in the Appendix). The simulations are initialised with fluid velocity and concentration fields null everywhere and particles randomly distributed on the surface of the bottom wall, with random in-plane orientations and velocities equal to zero. We let each simulation run for 2×10 6 time steps, corresponding to approximately T run ≈ 1200 − 5800 τ, where the characteristic time τ is the time an isolated particle takes to displace by one radius, i.e. τ = R/V p , and it varies in the range τ ≈ 350 − 1700 lbu (depending on the value of |µ 0 |). It is worth noticing that, typical experimental values of size and propulsion speed of SPCs are, respectively, of the order of few µm and few µm/s 14,46,47 , such that τ ∼ 1 s. This means that our run lengths are comparable with the duration of experiments which last typically few tens of minutes 9,26 .

Dynamic scenarios controlled by the phoretic mobility.
A number of experimental and numerical/theoretical studies of self-propelled particles in (quasi)-2d have given indication of the emergence of clustering 9,[26][27][28][48][49][50] , however the nature of the mechanisms determining the formation of aggregates lacks a consensual agreement and seems to be strongly system-dependent (see also 16 for a recent review). In our simulations solvent and solute hydrodynamics is fully resolved, from the far field down to distances of the order of the particle size (below which it is regularised by the lubrication interaction). We deal with spherical particles, which rules out the possibility of alignment-induced collective motion; instead, chemical production and diffusion mediate an effective interaction, analogously to the experimental system studied in 26,28 . While in the experiments it was surmised that active colloids felt an attractive interaction, here we can tune the affinity of the particles for the solute via the phoretic mobility µ 0 , which can be regarded as an effective charge 51 , i.e. positive/negative values induce repulsive/attractive interactions, respectively. Indeed, while for µ 0 < 0 our simulations confirm the formation of clusters, for µ 0 > 0 such cluster phase disappears, with the average cluster size going to zero. Incidentally, let us remark that, in some respect, suspensions of SPCs may recall other systems of interacting microswimmers, as, for instance, attractive squirmers 52 ; there are however at least two major differences: for SPCs, unlike squirmers, the characteristic selfpropulsion speed is constant only for an isolated swimmer, but in general it depends on the concentration field; the second and probably most important one, from the point of view of collective dynamics, is that while in the squirmer case particle-particle interactions are frozen (i.e. dictated by the interaction potential once for all), in a SPCs suspension phoretic interactions are mediated by the solute field and are, therefore, dynamical, in the sense that they depend on the local (in space and time) field configuration. In other words, phoretic interactions are not pairwise additive but change as a function of the global dynamics, and as such, they give rise to a collective behaviour that is genuinely out-of-equilibrium. In what follows, the phoretic mobility will be expressed as µ ≡ µ 0 /|µ * |, where |µ * | is the absolute phoretic mobility for which an isolated particle of radius R would have unitary Péclet number. To address the impact of the chemical affinity on the collective dynamics quantitatively, we have performed a Voronoj tessellation analysis of the particle space configurations 53,54 * .
The bottom insets of Fig. 2 show the Voronoj diagrams both for the repulsive and cluster-forming regimes; as clearly visible to the naked eye, the geometries of the Voronoj cells for chemoattractive and chemorepellent active colloids are distinctively different. The standard deviation of the cell area distribu- normalised by the square of the mean value S ) turns out to be a good indicator to distinguish the two types of dynamics. In the top inset of Fig. 2 we plot σ S (t), as function of time, for two cases with different sign of the phoretic mobility. In the attractive case, µ < 0, cluster formation induces the appearance of very small (and large) cells and, hence, the surface fluctuations grow and eventually saturate at long times. For positive µ, instead, colloids repel each other and tend to reach an optimal covering of the space, implying that σ 2 S (t) attains a (lower) value which remains constant in time. Correspondingly, the dependence of σ 2 (SS) S (the time average of σ 2 S (t) over the steady state) on µ discriminates between the two regimes: it is high for negatively large µ, decreases as µ approaches zero and then stays low and constant for µ > 0. We will focus, in what follows, on the chemoattractive case, but before moving on we stress that the phase diagram for chemorepulsive self-phoretic colloids, as recently shown theoretically and numerically 55,56 , is indeed rather complex and deserves further investigation.

Cluster statistics and morphology.
We first characterise the cluster size distribution of chemoattractive SPCs at changing the colloid/solute coupling intensity, |µ|. We identify clusters according to a distance criterion. Two particles share a bond if their centers are at a distance equal to or less than a cutoff apart ‡ , and we define as clusters, groups of particles connected to each other through a bond. We compute probability density functions (PDFs) of cluster sizes over the steady state of each run. Fig. 3 shows such PDFs, which can be in all cases fitted to an exponential, P(n) ∝ e −n/n c , over a wide range of sizes n. The characteristic value n c and the mean size n = 1 N clus ∑ N clus i=1 n i increase linearly with |µ| (see inset of figure 3), hence with the velocity of an isolated particle, in agreement with experimental and numerical observations 9,26,50 . The global attractor for the

SPCs dynamics is a set
where C i is the region of plane occupied by the i-th cluster, of area A i and containing n i SPCs. Correspondingly, the colloid number density reads Since the colloid density fluctuations can be expressed as where |D| is the measure of the whole plane. The number of particles in a cluster n is known to scale with the cluster gyration radius R g as n ∼ R d f g , d f being the fractal (Hausdorff) dimension 57,58 ; n i and A i are, then, related by A i ∼ R 2 g ∼ n 2/d f . Plugging the latter relation and P(n) ∼ e −n/n c into (9), and approximating the sum with the integral, we finally get: where σ 2 ρ 0 stands for the fluctuations for a inactive particles and a is a phenomenological parameter, and where we have used the relation n c ∼ |µ| § (see inset of figure 3). Fig. 4 shows the quantitative agreement of the predicted power law, with the correct scaling exponent, with the numerical observations.

Role of hydrodynamic interactions.
Self-propelled colloids interact due to both the chemicals they produce and the flows they induce. Understanding the relative magnitude and competition between these two sources of dynamic interactions remains challenging. The model put forward allows us to switch off the hydrodynamic interactions (HI), setting the fluid velocity to zero at each time step, yet keeping the self-phoretic mechanism and the correct translational and rotational hydrodynamic friction (see Appendix for further details). Interestingly, our study reveals that, although the different dynamic scenarios at changing the sign of the phoretic mobility are preserved even without HI (being mainly determined by the chemical interaction), HI have a profound effect in the kinetics of formation and morphology of the observed aggregates. In the absence of particle induced flows in the solvent, attractive SPCs (µ < 0) show an enhanced tendency to form clusters, as it appears in figure Fig. 5, where we compare the time evolution of the mean cluster size n(t), with and without HI (no-HI). In the no-HI case, clusters coarsen, with n(t) growing in time as t 1/2 (top right inset), as for domains in a spinodal decomposition. The same behaviour (n(t) ∼ t 1/2 ) has been observed, indeed, in § In principle there can be a dependence on µ also of the fractal dimension d f ; we assume here, however, that the change in µ affects only the characteristic cluster size and not its "compactness" (or fractality). simulations of self-propelled Brownian particles interacting via a shifted-truncated Lennard-Jones potential 49,52 . With HI, instead, coarsening is arrested, as observed in experiments 26 . Simulations have suggested that in suspensions of attractive squirmers the emergence of continuous or arrested coarsening is selected depending on the form and intensity of the active stress (the coefficient B 2 in the squirmer terminology) 52 ; self-phoretic Janus colloids behave, in this respect, as squirmers with B 2 = 0 59 , for which, indeed, arrested coarsening was observed 52 . Due to the non stationarity of the coarsening process, steady state PDFs of cluster sizes cannot be computed in the no-HI simulations. Nevertheless, we observe that istantaneous cluster size distributions F(n,t) (i.e. the number of clusters of size n at time t) tend to assume a self-preserving scaling form F(n,t) ∼ n −2 f (n/n(t)), as it happens in classical colloidal aggregation phenomena for massconserving systems 60 . This is shown in the top left inset of Fig. 5, where we plot n 2 F(n,t) vs n/n(t) and see that all sets of points for different t's in the coarsening regime where n ∼ t 1/2 , within error bars, collapse onto each other. The statement on the different dynamics, with and without HI, is corroborated by the inspection of the radial distribution functions (RDFs) 61 (indicated as g HI (r,t) and g no−HI (r,t), respectively), defined as the probability of finding a particle between the distances r and r + dr from a reference particle (and averaged over all particles), i.e.
g(r,t) = 1 when hydrodynamics is active. Besides, clusters appear substantially more compact, as appreciated in the snapshots (insets) and quantified by the measurement of a larger fractal dimension (d  4). Hydrodynamics hinders, then, the colloidal aggregation process. Several complex mechanisms can be conjectured to cause this phenomenon: dynamically induced effective repulsion among particles, fluid flow generated disturbances in the chemical field distribution, etc. An effect, that we could clearly identify, is an enhanced tendency of SPCs to be oriented off-plane, when HI are present, not only when they hit a cluster (as, e.g., in the mechanism proposed in 62 ) but also for isolated particles. This may be attributed to fluid motion close to the wall, giving rise to hydrodynamic torques that rotate the particles. Actually the roto-translational dynamics of selfdiffusiophoretic colloids, at and close to interfaces, is an intricate problem 63-65 : in highly confined situations, one might indeed expect even the opposite trend (clusterisation enhancement) 62 ; the argument maintains, therefore, a qualitative character. Nevertheless, Fig. 6c provides a quantitative insight to the picture. There, we show the PDF of the degree of alignment of the particle orientation with the bounding solid wall, m 2 || = m 2 x + m 2 y . When HI are present, indeed, the peak of the PDF around m 2 || ∼ 0 is more pronounced, i.e. there is a larger fraction of colloids pointing out of the plane. Accordingly, the self-propulsion speed is effectively reduced, thus limiting the in-plane mobility and diminisihing the capability of particles to gather and clusterise. Before concluding, we would like to remind that it is still an open issue whether Janus particles can really have a homogeneous phoretic mobility; if the opposite is true, inhomogeneous µ(r) gives rise, in response to gradients of the concentration field, to chemical torques that can lead, themselves, to clustering inhibition [66][67][68] . With hydrodynamic interactions, the dynamics is, of course, even more complicated, due to the competition of these effects, and it is subject of ongoing research.

Conclusions
To conclude, we have used a mesoscopic numerical model of fully resolved spherical active colloids, propelled by self-generated gradients of a scalar field (e.g. a chemical product) where the selfinduced hydrodynamic flows can be accounted for. We have identified the role of the phoretic mobility as the key controlling parameter that determines two distinct dynamic regimes and the onset of a cluster phase. By means of a Voronoj tessellation we have characterised the cluster state finding that the probability distribution of sizes decays exponentially with a mean size growing linearly with the particle activity, in agreement with experimental results 9,26 . We have quantified the profound effect that hydrodynamics plays inhibiting clustering for negative phoretic mobilities. We have identified the interplay between induced flows and particle reorientation as a possible explanation to the strong slowing down of cluster coarsening, although it remains an open question, which needs a deeper analysis, whether fluid-wall interactions dominate over particle-particle hydrodynamic correlations. This study shows that our novel numerical method is powerful and enjoys some unique features, namely the explicit description of chemical signalling, through the production and diffusion of a solute concentration field and the solvent hydrodynamics, to simulate realistic systems. Moreover, it opens the way to address the dynamics of self-propelled colloids in general geometries and also for stronger activity (larger Pe), both in isotropic and unforced situations, where aggregation can lead to the formation of active colloidal gels, or under gravity as in the experimental sedimentation setup.

Conflicts of interest
There are no conflicts to declare. FP/2014-2020)/ERC Grant Agreement No. 739964 (COPMAT). This work was possible thanks to the access to MareNostrum Supercomputer at Barcelona Supercomputing Center (BSC) and also through the Partnership for Advanced Computing in Europe (PRACE).

The lattice Boltzmann method (LBM)
The fundamental quantities, in the LBM, are the discrete single particle probability density functions (pdfs), f l (x,t), defined in such a way that f l (x,t)∆x 3 is the probability of finding a fluid particle (of unit mass) in a small volume ∆x 3 , centred in x, at time t, moving at a velocity c l ; the index l runs over the set of discrete lattice velocities (nineteen in our case) 31,69 . The evolution equation for the f l 's is the so called lattice BGK-Boltzmann equation 70 : where ∆x and ∆t are the lattice spacing and time step, respectively, and τ LB is the relaxation time. Algorithmically, Eq. (12) consists of two steps: i) streaming, where the pdfs move ("stream") over the lattice and ii) collisions, which induce relaxation on the local equilibrium distribution functions, f (eq) l (x,t). The latters depend on space and time only through their prescribed dependence on the hydrodynamic fields, density ρ(x,t) and velocity u(x,t). We adopt, for the f (eq) l 's, a standard second order polynomial form 30 , that is: where w l and c s are constants (the lattice weights and speed of sound) characteristic of the particular lattice geometry and velocity scheme. The density and velocity fields are expressed, in terms of the f l 's, as: In the small Knudsen and Mach numbers limits, Eqs. (12)(13)(14) represent a numerical solver (second order accurate in space and time 69,71 ) for the incompressible Navier-Stokes equations.

Fluid-solid coupling
Coupling the solvent dynamics with that of the suspended particles requires tackling a boundary condition at the solid-fluid interface. A solid particle is discretised on the lattice, so that its surface crosses links between neighbouring lattice nodes. Let us assign a generic link the label j ( j = 1, . . . , N links ) and denote the associated outer and inner node as x o and x i = x o + c l j ∆t, where c l j is the corresponding lattice velocity (l j = 1, . . . , 18 * If a no-slip condition were to be enforced at the solid surface, the boundary probability density functions would be updated as follows: where the index l j is such that c l j = −c l j and f l j is the postcollisional distribution. However, the boundary velocity at a given location on the particle surface will be in general non-zero (due to slip and/or particle motion). The rule (15) must be, then, modified accordingly, as 41,42,44 : where ρ 0 is the mean fluid density and the local velocity u ( j) s at the particle surface location x ( j) s = x o + 1 2 c l j ∆t (the position of the boundary node is assumed to be in the middle of the link) depends on the particle linear V and angular Ω velocities and on the (local) slip velocity v where X is the position of the particle centre of mass. The boundary condition (16-17) induces a local fluid-particle momentum exchange (the total momentum being conserved) which entails a force at the boundary node x The total force and torque acting on the particle are, then, given by: where the sum runs over all links constituting the particle surface. Let us notice that, because of (17-18), such force and torque can be both split into four terms 44 : The expressions of F h and T h are as follows: The latter equations provide the force and torque exerted by the fluid on the particle and as such they mediate also the hydrodynamic interactions. The friction matrices and are diagonal, whereas are null (apart from discretisation effects, coming from surface irregularities, which can be reduced by increasing the particle radius). In what we call the 'no-HI' case (without hydrodynamic interactions) the effect of the fluid enters, then, in the drag terms −ζ F V and −Γ T Ω. The phoretic force F p and torque T p are given by: In order to get a better insight on the effect of these terms on the particle dynamics, we consider first the case of an axisymmetric configuration of the chemical field C(x) around the particle. This is the case, for instance, of the active Janus self-diffusiophoretic colloid we consider. The sum in (28), which runs over all links (or, equivalently, over all boundary nodes), can be split into two sums: a sum over the subset I x o of links connecting a certain fluid node x o , proximal to the particle surface, to the corresponding solid nodes inside the particle, and a sum over all such 'outer' fluid nodes, i.e.
The particle surface is endowed with a spherical coordinate system, such that the generic point x o can be associated to the usual polar and azimuthal angles, (θ , φ ), with respect to the symmetry axis identified by the the particle directorm (as in the main text). Bê 32) the triple of orthonormal vectors corresponding to the normal and the two tangential directions to the surface in x o † , the following decomposition c b = (c b ·n)n + (c b ·t 1 )t 1 + (c b ·t 2 )t 2 obviously holds for the lattice velocities, with b ∈ I (θ ,φ ) (we will use for x o its spherical coordinates from now on). We recall that the effective slip velocity, v s ∝ (1 −n ⊗n) ∇C, is, by definition, tangent to the surface, therefore: where v s,1 (θ ) and v s,2 (θ ) are the components of v s along the two tangent directions (polar and azimuthal), which depend only on θ owing to the axisymmetric character of the field C(θ ). Using (31) and (33), the force (28) can be rewritten as: s ∆t . Let us note, now, a property of the set of links: given x o = (θ , φ ) on the sphere, and its associated I (θ ,φ ) , there will be a x o = (θ , φ ) (with θ = θ and φ = φ + π) such that for every b ∈ I (θ ,φ ) there exists a b ∈ I (θ ,φ ) fulfilling where the superscripts "||" and "⊥" stand for the directions, respectively, parallel and orthogonal to the symmetry axism. Such property follows from the definition of I (θ ,φ ) , equation (30): since, by symmetry, upon the shift φ → φ + π the normal to the surfacen changes as n || = n || and n ⊥ = −n ⊥ (see (32)), in order to preserve the condition c b ·n < 0, the velocities have to transform as in (35). (34), upon summing on φ (i.e. on nodes all around the particle at a fixed latitude θ ), only the component parallel tô m survives, which reduces to: if the particle is large enough (R/∆x 1), so to minimise the discretisation errors, the coefficients are a θ ∝ 2πR sin θ , such that the sum in (36) approximates the integral over the polar angle. So, basically, for a solute distribution around the particle which is axisymmetric, the bounce-back-on-links algorithm yields a phoretic force which is proportional to the surface integral of the effective slip velocity, i.e.: It is easy to verify, by similar symmetry arguments, that the torque in (29) is T p = 0. In the general case it is not so straightforward to estimate such force, however we can grasp qualitatively what happens. We can always decompose the solute field C(x) around a (Janus) particle as the sum of the self-generated field (which is axisymmetric) and a background field (which will not be, in general, axisymmetric), C(x) = C axisymm (x) +C background (x). If we now assume that the background field does not vary significantly on the scale of the particle size (that is C 2 1/2 / ∇C(X) R), we can approximate it as C(x) ≈ C(X) + ∇C| X · (x − X). To leading order, then, the phoretic force F p will also split into a self-propulsion contribution and in a termF p ∝ ∇C| X ; the latter represents a sort of chemotactic drift towards regions of high or low solute, depending on the sign of the phoretic mobility.

Switching off hydrodynamics
When assessing the dynamics where no hydrodynamic interactions are present (see section 3.3), we set the fluid velocity field u(x,t) to zero everywhere. Consequently, the equation for the solute field C(x,t) reduces to a diffusion equation and the force F h and torque T h are null. To understand this latter point, we have first to notice that what enters in the expressions (22) and (23) are the post-collisional distribution functions. For τ LB = ∆t (as it is in all our simulations), these coincide exactly with the equilibrium distributions, equation (13), with u ≡ 0, that is f l (x,t) = w l ρ, ∀l. This implies: which, by symmetry, vanish (again, apart from small errors due to the surface discretisation) when summed over all boundary link velocities around the sphere. Let us remark that setting the fluid velocity to zero does not affect directly ‡ the expression for the phoretic force, equation (28), therefore the self-propulsion and phoretic interactions survive the operation of switching-off hydrodynamics. In this way we can decouple, in our mesoscopic framework, physical phenomena occurring at well separated scales: the fluid dynamic processes localised (on molecular scales) close to the particle (and effectively accounted for in the solute-field-dependent slip velocity) from macroscopic flows in the solvent (and associated hydrodynamic interactions). ‡ It affects it, of course, indirectly, since the dynamics of C(x,t) differs.

Particles close to contact: lubrication corrections and short-range repulsion
For particles close to contact, lubrication corrections are introduced: the forces and torques acting on two particles approaching each other are calculated, in terms of particle velocities and angular velocities, according to a grand-resistance-matrix formulation 44,45,72 . In particular, the lubrication correction takes the form of the difference between the lubrication force at a surface separation h and the force at a given cut-off separation h c ; for two particles of radii R 1 and R 2 (the particle-wall interaction corresponds to the limit R 2 → ∞) this reads 44 : where r 12 = X 1 − X 2 ≡ r 12r12 is the particle center-center distance vector and h = r 12 − R 1 − R 2 ; the cutoff distance is chosen to be h c = 0.67 lattice units, which is an optimal value to get good agreement with lubrication theory calculations, as shown in 44 . Lubrication forces may not be enough, though, to prevent particle overlap (as recognised also in 73,74 ), especially when the particle density is large (even just locally, as, for instance, inside clusters). Therefore, we add also a short-range soft-sphere repulsion modelled by the force F ss ∝ (h ss c /h) 3 − 1, with cutoff (coinciding with the soft-sphere radius) h ss c = 2h c .

Numerical tests
In this section we present results form simulations of simple test cases in order to check to which extent our method agree, qualitatively and/or quantitatively, with the theoretical predictions. Firstly, we consider a single SPC with constant phoretic mobility, µ(x s ) ≡ µ 0 > 0, and a Janus activity profile α(x s ) = α 0 ifm ·n ≤ 0 0 ifm ·n > 0, moving in a fully periodic tridimensional box. Such a particle is expected to perform a uniform rectilinear motion with speed 23 In figure 7 we report the coordinates of the position and the three velocity components (in lattice Boltzmann units) for an SPC of radius R = 2.5∆x in a cubic box of side L = 32∆x, initially located in X(0) = (L/2, L/2, L/2) with orientationm(0) = (0, 0, 1): the plot suggests that indeed the particle moves on the straight line, with director (0, 0, 1), with constant speed (the oscillations are due to recomputation of the set of links, hence of the discretised surface, as the particle moves across the lattice), albeit about half of the theoretical value. Such a quantitative disagreement is not too surprising, though, given the small size of the particle. The quality of the numerical result improves, in fact, significantly with the resolution. This is shown in figure 8, where we plot the steady state z-component of the velocity, V z , normalised by the expected V p , equation (41), as a function of the radius (the ratio R/L ≈ 0.08 is kept fixed over the various runs). We see that the numerics (red  (41)), for a single Janus SPC in a periodic box as a function of the particle radius (a dashed line at the value of 1 is plotted as a guide for the eye).
dots) approximate better and better the theoretical prediction as the particle radius is increased, and eventually tend to converge to the value ∼ 1.1 (i.e. the measured V z is about 10% larger than V p ). This residual mismatch might be due to the fact that, as R grows, the Péclet number, Pe, approaches unity, whereas the prediction (41) is valid for Pe strictly zero. For a smaller Pe ≈ 0.1 we find, indeed, a better agreement for the largest radius (blue square in figure 8), for which V z /V p ≈ 1.04. In the multi-particle case we decided to pay the price of a worse agreement with the single particle theoretical speed and adopted a smaller radius such to be able to simulate statistically significant systems with an affordable computational cost.