Orientational order and morphology of clusters of self-assembled Janus swimmers

Due to the combined effect of anisotropic interactions and activity, Janus swimmers are capable to self-assemble in a wide variety of structures, many more than their equilibrium counterpart. This might lead to the development of novel active materials capable of performing tasks without any central control. Their potential application in designing such materials endows trying to understand the fundamental mechanism in which these swimmers self-assemble. In the present work, we study a quasi-two dimensional semi-dilute suspensions of two classes of amphiphilic spherical swimmers whose direction of motion can be tuned: either swimmers propelling in the direction of the hydrophobic patch (WP), or swimmers propelling in the opposite direction (towards the hydrophilic side) (AP). In both systems we have systematically tuned swimmers' hydrophobic strength and signature, and observed that the anisotropic interactions, characterized by the angular attractive potential and its interaction range, in competition with the active stress, pointing towards or against the attractive patch gives rise to a rich aggregation phenomenology.


I. INTRODUCTION
In the last decade, the pioneer experiments by Dombrowski and Cisneros [1,2] inspired the numerical work aimed at understanding the influence of hydrodynamics in the formation of coherent structures of microorganisms, simulated as spherical squirmers [3] in 3D [4,5] and quasi-2D [6].The features of a 3D squirmer suspension have been numerically characterized using the smoothed profile method [7][8][9] and Lattice Boltzmann (LB), either in bulk [10] or in response to a steady Couette flow [11].While Evans et al. [12] used Stokesian Dynamics to detect the instability of the isotropic state in a 3D suspension of weak pullers, Delmotte and coworkers [13] applied a force-coupling method to investigate, by large-scale simulations, the emergence of a polar order in a 3D squirmer suspension.
More recently, experiments by Ref. [14][15][16][17][18] have inspired numerical simulations of 2D and quasi-2D suspensions of squirmers interacting via a repulsive potential.* cvaleriani@ucm.esWhen a suspension (simulated by means of multi-particle collision dynamics) was confined between parallel walls (quasi-2D geometry) [19], active particles were forming clusters whose size depended on the system's concentration with the largest being made of pullers.The suspension prepared in such confined geometry would undergo phase separation into a dilute and a dense phase [20], differently to the 2D suspension of disk-like squirmers [21], where hydrodynamics suppressed phase separation.
Inspired by the self-assembly observed in experiments with either bacteria [16] or active colloids [17] and by numerical works on 2D dilute suspensions of active brownian particles (ABP) interacting via isotropic short-range attractions [22][23][24][25], we have studied a semi-diluted quasi 2D suspensions of attractive squirmers [26], to unravel the role played by hydrodynamics in suspensions of attractive active particles.
Most of the nowadays synthetic active colloids are colloids partially coated with Pt.These colloids self-propell by a catalytic reaction of H 2 O 2 and O 2 taking place at the Pt surface.This mechanism of colloid's propulsion is due to phoresis.On the one side, properly modeling the phoretic interactions can be quite complicated [27][28][29][30][31][32], eventhough there are experiments, simulations and theo-ries that consider them relevant for particles' propulsion.On the other side, which interactions dominate in these systems are still a matter of debate [33].
From the experimental point of view, several studies have focused on understanding the effect of hydrophobicity in catalytic Janus colloids.Modifying the surface of the particles, the authors of Ref. [34] found that particles with an hydrophobic patch would move faster than the ones with an hydrophilic patch.Similarly, Gao et al [35], studied self-assembly of amphiphilic particles made with hydrophobic octadecyltrichlorosilane (OTS)-modified silica microspheres capped with a catalytic Pt hemisphere patch.Additionally, Yan et al [36] have recently reported different collective states of active dipolar Janus colloids, whose motion is originated by an induced-charge electrophoresis.These colloids self-assemble in different morphologies, ranging from: gas, swarms, chains and clusters, by modifying both the charge imbalance of the colloids and the electric field frequency.Moreover, in [36], the authors were able to numerically reproduce all observed states by means of an overdamped molecular dynamics simulation whose particles were interacting via directed imbalanced interactions.
Motivated by these two experimental works [35,36], simulations of a dilute quasi-2D suspension of ABP interacting via anisotropic amphiphilic Janus interactions have been reported in Ref. [37].To the best of our knowledge, this is the first numerical attempt to unravel the collective behaviour of amphiphilic active Janus particles, taking into account their anisotropic nature.In Ref. [37] both the effect of interaction directionality and the propulsion speed have been used to control the physical properties of the assembled active aggregates.One might suggest that these two parameters could be easily tuned also in experiments, given that the concentration of a cationic surfactant could be used to reverse the propulsion's direction and adding pH neutral salts could be used to control the propulsion speed [38].
In order to establish the relevance of hydrodynamics in this suspension, we study the collective behaviour of quasi 2D suspensions of spherical squirmers interacting via an anisotropic Janus potential, considering different interaction ranges, interaction strengths and hydrodynamic signatures.When clusters appear, we characterise their morphology as a function of the above mentioned features.The objective of this manuscript is the analysis of minimal models that contain the essential ingredients of activity and steric interactions for heterogeneous particles.The scenarios we identify can be clearly correlated with the competition between self-propulsion, stress generation, and attraction among particles.It is true that phoresis is more complex and can also involve long range interactions among the active elements.Our results serve as a benchmarking scenario if phoresis is controlled by near field effects.In this sense, the systematic analysis we have carried out provides a valuable contribution in this general perspective.In a more specific perspective, this study also helps to understand the effect of the hydrodynamic interactions in the collective dynamics of active Janus colloids, inspired by experiments showing flow fields around a single catalytic Janus sphere [39].
The manuscript is organised as follows.In section II we present the anisotropic interaction model, the squirmer model, the lattice Boltzmann methodology and simulations details.In section III we describe the tools used to characterise the Janus squirmers.In section IV we present all results and discuss our conclusions in section V.

II. SIMULATION DETAILS
In this work, we have studied a quasi two dimensional dilute suspension of squirmers [40][41][42] interacting via a Janus anisotropic potential similar to the one used in Ref. [37] for Active Brownian Particles (ABP).In order to model the Janus swimmers, we have introduced a pair potential V (r ij , θ i , θ j ) which depends on the distance between two squirmers (r ij ) and their attractive patches' orientations where θ i is the angle between the patch unit vector p i and the inter-particle vector r ji = r j − r i ; and θ j the angle between p j and r ij = − r ji .
V rep (r ij ) represents a short range soft repulsion [43] (to avoid overlapping) given by where s = 1.5 is the energy scale, ν = 2, r s ≡ r ij − σ is the inter-particle distance (with r s > 0 to ensure that particles never overlap), σ s = 0.5 is the characteristic separation length between particles, and h 0 the distance at which the potential is shifted to ensure its continuity at r s = h 0 .
[? ] The attractive potential consists of a radial V att (r ij ) and an angular φ(θ i , θ j ) contribution.Concerning the radial part of the potential, we have considered two different truncated Lennard-Jones: one long (r c = 2.5σ) and another short (r c = 1.5σ) range.While the following equation: represents the former, it has been modified for the latter by adding a smooth transition function, so that the potential V att and its derivative vanish continuously as r ij approaches r c [44] (see appendix A).Therefore, when-ever (2. whereas for r ij < (2.2 1/6 σ) V att = −0.305and V att = 0 otherwise.Inspired by the work of Hong et al. [45] for amphiphilic colloidal spheres, the anisotropic interaction is controlled by the angular contribution of the attractive potential where μij = r ij /r ij is the unit vector joining the squirmers' center of mass.This anisotropic potential ensures a non-vanishing torque when (p i • μji ) > 0 and (p j • μij ) > 0 (a more detailed analysis is discussed in appendix B).Within our model, each Janus particle self-propels along a direction êi rigidly bounded to the particle.As in Ref. [37], we study two types of active Janus particles depending of the orientation of êi with respect to pi .When the propulsion direction (ê i , red arrow in figure 1-left panel) is pointing towards the attractive patch pi (red arrow in figure 1-left panel), swimmers are propelling "with the patch" (WP).When the propulsion direction (ê i , red arrow in figure 1-right panel) is pointing against the attractive patch −p i (green arrow in figure 1-right panel), swimmers are propelling "against the patch" (AP).The swimming direction ê is the red arrow and direction of the attractive patch p is the green one.The center-to-center distance between two swimmers rij is the black vector.
A spherical Janus squirmer of radius R p is characterized by the tangential surface slip velocity, where τ is a unit vector tangential to the particle surface, B 1 quantifies the asymptotic self-propelling speed of an isolated squirmer (v s = 2 3 B 1 ) [46,47] and B 2 is proportional to the active stress generated on the surrounding fluid.The ratio between active stress and self-propelling velocity, β = B 2 /B 1 [3], quantifies the active state of the squirmer and its interaction with the fluid.A positive value of β corresponds to the case when thrust is generated in front of the squirmer's body (puller); a negative value of β to the case when thrust is generated at the back (pusher) [4] and β = 0 corresponds to a neutral swimmer (see Fig. 2).In our simulations we disregard thermal fluctuations; thus velocity fluctuations are simply induced by the particles' activity, where B 2 can be understood as an effective temperature.
In order to characterize how hydrodynamics affects the assembly of a dilute suspension of Janus swimmers, we model the embedding solvent using a Lattice Boltzmann (LB) method with a D3Q19 three dimensional lattice [48,49], implemented in a highly parallelized code [50] that exploit the excellent scalability of LB on supercomputing facilities.Spherical squirmers (with a diameter σ of 4.6 lattice nodes [51]) are individually resolved, imposing a modified bounce-back rule on the oneparticle velocity distribution functions to nodes crossing the solid particle's boundary (including the slip velocity to impose the squirming motion [52,53]).The total force and torque the fluid exerts on the particle is obtained by imposing that the total momentum exchange between the particle and the fluid nodes vanishes (as a result of the modified bounce-back).While accounting for all forces acting on each squirmer allows to update them dynamically, the torque exerted by the fluid determines the change of the squirmer 's self-propulsion direction.
The simulated system consists of a suspension of N = 888 spherical squirmers at φ = 0.10 (where φ = π 4 ρσ 2 and ρ = N/L 2 ) in a quasi two-dimension geometry with periodic boundary conditions (L×L×kσ, with k = 5 and L ≈ 83σ) needed to capture three dimensional hydrodynamics effects, whereas both colloids' position and orientation are confined to move in 2D.Additionally, we have carried out bigger simulations with N = 3552 squirmers at the same packing fraction and L ≈ 167σ, in order to proof that results are not due to finite size effects, specially cases with large cluster sizes.
As in Ref. [26], we quantify the competition between attractive and self-propelling forces via the dimensionless parameter where F d = 6πηR p v s is the friction force associated to the squirmer self-propulsion and F att is the absolute value of the attractive force at its minimum, corresponding to r = (26/7) 1/6 σ = 1.245σ for the long-range potential in Eq. 3 and to r = 1.239σ for the short-range potential in Eq. 4. On the one hand ξ controls the competition between two mechanisms: the self-propulsion and the interaction strength.On the other hand the mechanisms that control the re-orientation of the particles are the active stress given by the B 2 squirmer parameter and the orientational-dependence of the potential which is modulated by the range and the strength of the interaction.Given that a squirmer travels its own size in a time τ = σ vs , we perform simulations from 1450 up to 3000 τ .Once the system reaches steady state, in a time window between 1000 and 2000 τ , we carry out a systematic analysis of the dynamics of the squirmer suspension, considering ξ from 0.1 to 10 and β from -3 up to 3.

III. ANALYSIS TOOLS
In order to be able to establish the effect played by the anisotropy of the interaction between swimmers, we follow the procedure in Ref. [26] and study the degree of aggregation as a function of β and ξ.Whenever the system forms a steady-state cluster distribution, we explore the clusters morphology by computing their mean size, the cluster-size distribution, their radius of gyration, and both the polar and nematic order parameter.
When clusters are present, we identify them following a distance criterion: at each time step, two particles belong to the same cluster of size s whenever their distance is smaller than r cl = 1.75σ (first minimum of the g(r), see appendix C).
To calculate the cluster-size distribution f (s) we apply a criterion based on [54]: (i) We arbitrarily subdivide the range of s-values into intervals ∆s i = (s i,max − s i,min ), where n t i is the total number of clusters within each in-terval ∆s; (ii) we assign the value n i = n t i /∆s i to every s within ∆s i , and compute the fraction of clusters of size s as f (s) = n i /N c , where N c = i n i ∆s i is the total number of clusters.
To characterise the clusters' morphology, we compute the radius of gyration the average polar order parameter for each cluster size s the tensor order parameter (for a 2D system [55]) (being h and k equal to x, y and N the total number of squirmers) and the nematic order parameter λ(t) being the largest eigenvalue of Q hk [56].
We have also computed both global nematic and global polar order parameters, substituting s with N in the global polar order parameter (Eq.9).Therefore, the polar order in steady-state is P ∞ = P (t >> 0) while the nematic order in steady-state is λ ∞ = λ(t >> 0).

A. Mean cluster size
We start with computing the mean cluster size s as a function of time for different values of β and ξ.This allows us to distinguish coarsening from clustering, establishing, in the latter case, when the system reaches a steady state.
Fig. 3 represents the mean cluster size as a function of time for suspensions where attraction dominates over propulsion (ξ = 0.1) of AP squirmers (interacting via long-range (panel-a) or short-range (panel-b) attraction) and of WP squirmers (interacting via long-range (panelc) or short-range (panel-d) attraction).In all AP suspensions the system coarsen, within the simulated time window as shown by the continuously growing mean cluster size, with a speed dependent on β.Pullers coarsening faster than pushers, a part from the system with β = −3 where particles form clusters independently on the interaction range.WP suspensions only coarsen when interactions are long-range (panel-c) with a speed dependent on β and pullers coarsening faster than pushers, even though coarsening is much slower than in the AP case.In panel-c, the coarsening dynamics is so slow that s < 10 even when t > 1000τ (it is worth noting that we have never observed a steady state despite having run up to 6 × 10 6 LB-time steps, corresponding to ∼ 8700τ ).Comparing panel-c and panel-d, we conclude that the attraction range affects aggregation, given that when the interaction is short-range (panel-d), coarsening is suppressed and s reaches a steady state at long t for all values of β.A more detailed analysis is given in appendix D for the case of the WP system with ξ = 0.1.AP squirmers interacting via long-range attractions (panel a) quickly reach a steady state characterized by s ≈ 2 except for β = 3 where s ≈ 10 (red curve in panel-a).Whereas when dealing with AP squirmers interacting via short-range attractions (panel b), pushers reach a steady state quite quickly with s ≈ 2 as well as pullers with β = 3 (with a slightly larger mean cluster size of s ≈ 2.5).Interestingly, the mean cluster size for weak pullers (black curve in panel-b) first develops a metastable state at short times (t ∈ (10, 500)τ ) and then another one for longer times (t > 1000τ ): the latter due to the formation of dynamic clusters with polar order.Similarly, when β = 0 the decay in s at long time corresponds to monomers aligning and moving in the same direction.
WP squirmers with β < 3 relatively quickly reach steady-state with clusters with an average of 4 particles (panel-c).Whereas pullers with β = 3 reach a slightly higher mean cluster size s ∼ 5: such increase in the cluster size is observed only for this interaction range, since WP squirmers with short-range attractions (panel d) develop a mean cluster-size of ∼ 4 particles, independently on their hydrodynamic signature.When self-propulsion dominates over attraction, squirmers in general form smaller clusters (never larger than 4 particles on average).In either AP suspensions, the mean cluster size for weak pullers (black circles in panels-a and b) is larger than the one for any other β, fluctuating around a value of 3 for long times.In the case of neutral squirmers and large ξ and after a long time (t > 1000τ ), particles in the suspension become strongly aligned and most of the clusters are monomers, for both interaction ranges.
When dealing with WP squirmers and high ξ, we have observed a small value of the mean cluster size (around 2), which indicates that particles form dimers on average.A similar behavior is observed for WP squirmers with short-range interactions, with slightly larger fluctuations in the mean cluster size for weak pullers at long times.
Having established the conditions needed for clustering to appear, we represent in Fig. 6  As shown in figure 6, when ξ = 1 WP swimmers mostly form tetramers and trimers, while WP swimmers with ξ = 10 form dimers.The clusters observed for WP squirmers have an average size rather insensitive to the squirmer hydrodynamic signature for these two interaction strength ξ = {1,10} and the two interaction ranges studied here.
For lower values of ξ, the dependence on β becomes more evident for WP squirmers with ξ = 0.1 (gray squares on panel-b): while for extreme values of β (either pushers or pullers) the clusters on average form pentamers, whereas neutral and weak pullers or weak pushers are mostly form tetramers.
AP squirmers with r c = 2.5σ (circular symbols panela) show that the mean cluster size increases as β increases when ξ = 1, whereas it does not change significantly when ξ = 10 or when the interaction range is shorter (circular symbols panel b).
As summary, on the one side, comparing the two panels in Fig. 6, when ξ = 1 and 10, the mean cluster size s for WP squirmers with interaction range of r c = 1.5σ behaves in the same way as the one with a larger cutoff (r c = 2.5σ).Therefore, when attraction does not dominate, the mean cluster size for WP squirmers is independent on the interaction range.On the other side, the mean cluster size s for AP squirmers with ξ = 1 and r c = 1.5σ (red circles in figure 6-(b)) does not behave in the same way than when r c = 2.5σ (red cirlces in figure 6-(a)), the former showing a systematic lower value of s for all β while the latter a monotonic increase of the cluster size for pullers.Finally, AP squirmers with ξ = 10 (black circles in figure 6) show the same behavior varying β, independently on the attraction range.With a small peak for weak pullers, and a minimum at β = 0.

B. Coarsening, clustering and aligned states
Fig. 7 displays the different states, types of collective motion or self-assembly that WP and AP Janus squirmers form as β and ξ vary.This way to classified our simulations, allow us to observe in a better way all the results.
We have identified seven different cases: 3 gas systems, 3 clustering systems and 1 coarsening state.All the different cases, show different types of cluster morphologies and/or particles' alignment.The emerging structures result to be sensitive to the patch direction (WP or AP), the hydrodynamic signature β and the interaction range r c and strength ξ.The three gas systems are characterized by the different particles' alignment: isotropic (with particles randomly swimming), polar (with most of the particles swimming in the same direction) and nematic (with approximately half of the total amount of particles swimming in one direction and the other half in the opposite direction).The three clustering systems can be classified as: finite-size dynamic clusters, chains and trimers.In the coarsening case all particles form a unique macroscopic aggregate.
AP Janus squirmers coarsen into a single cluster for ξ ≈ 0.1, (see Fig. 7-( favoring the formation of larger finite-size clusters, accelerating coarsening).When ξ = 10 , AP squirmers either form a polar ordered gas (see Fig. 7-(c), for weak positive stresslets, 0 < β < 1) or an isotropic gas (for AP pushers and pullers with β > 1).WP Janus squirmers have a strong tendency to form small micelles of three or four particles pointing at each other.For ξ ≈ 0.1 these micelles coalesce and end up forming chains, which in most of the cases coexist with smaller micelles (see Fig. 7-(d)).At ξ = 1 the competition with the active stress (β) gives rise to trimers and tetramers avoiding the formation of larger chains of particles (see Fig. 7-(e)).If activity dominates (ξ = 10) and |β| < 1, even though a polar order is not favored due to the anisotropic interactions, we still observe a global nematic order: a significant portion of particles move in one direction and the rest in the opposite one (see Fig. 7-(f)).For other values of β, we observe an isotropic system.In the Supplemental Material [57] we have included the movies of the six simulations represented in Fig. 7, to better capture both the dynamical and morphological features reported in this paper.

C. Cluster size distribution
Additionally to the mean cluster size, we calculate the cluster size distribution (CSD) to study the statistics of the different cluster sizes found it for Janus swimmers, computing the CSD for ξ = 1 and ξ = 10 and as a function of β, the direction of the attractive patch and the interaction range.The case when attraction dominates with respect to self-propulsion (ξ = 0.1) for AP squirmers is analyzed in appendix D. Results for ξ = 1 are presented in Fig. 8 and for ξ = 10 in Fig. 9. Following previous studies [26], where CSD was calculated, we use the same analytical function as a reference.
with γ 0 , s 0 , z 0 and B constants such that A = 1 − B.  Accordingly, CSDs with β < 3 can be described by Eq. 11, with B = 0 while the CSD for β = 3 requires a nonvanishing value of B in Eq. 11 to capture the observed maximum.The width of the CSDs is measured by the value of the exponential tail parameter s 0 , thus pushers CSDs have smaller s 0 , whereas the corresponding s 0 for pullers is larger until β = 3.At larger values of β the behavior changes due to the appearance of a second peak in the CSDs.
Fig. 8-(b) represents the CSDs of AP squirmers for ξ = 1 and r c = 1.5σ.CSDs of AP pushers and neutral squirmers follow the analytical behavior of a power law with an exponential tail (Eq.11 with B = 0: gray dashed line) with the same width.The squirmers with β = 3 present a CSD wider than the pushers ones, whereas weak pullers CSD follows a power law behavior instead of the power law with an exponential tail (Eq.11) which represents clusters of all sizes, such effect is not observed when interaction range is 2.5σ.
Fig. 8-(c) represents the CSDs of WP squirmers for ξ = 1 and r c = 2.5σ.A peak appears between trimers and tetramers, while the characteristic CSD width increases with the active stress β.Interestingly, when β < 3 monomers are not observed.Strong WP pushers are characterized by a suspension dominated by trimers, while strong WP pullers form a polydisperse suspension of monomers, dimers, trimers and even chains of tens of particles.
Fig. 8-(d) represents the CSDs of WP squirmers for ξ = 1 and r c = 1.5σ, with CSDs presenting a peak between trimers and tetramers.In this case, the CSD width depends on the active stress magnitude rather than in its sign: the higher the value of |β|, the wider the distribution (with distributions wider than the ones with longer interactions range r c = 2.5σ, and with a higher number of monomers especially when β = 3 ).
Therefore, CSDs present relevant differences depending on whether squirmers are AP or WP-like and on the interaction range of the potential.The CSDs for AP squirmers decay monotonously as a function of cluster size, with a large tail corresponding to cases in which there is the formation of finite-size aggregates.On the contrary, CSDs for WP usually display a peak for trimers.Moreover, the interaction range affects how the CSDs width depends on the hydrodynamic signature: while for large interaction range, strong AP pullers develop a wider distribution than any other β, for shorter interaction range weak pullers, develop a power law distribution.
In contrast, neither the range of interaction among WP squirmers nor the hydrodynamic signature affect the shape of the CSDs, since CSDs with a peak on trimers is present in all cases, with the widest distribution appearing for the stronger WP or AP pullers.
Fig. 9-(a) represents the CSDs of AP squirmers with ξ = 10 and r c = 2.5σ.The corresponding distributions are monotonically decreasing, particularly for β = 0.5 a power law distribution emerges while in the other cases, the CSD can be approximated by Eq. 11 with B = 0.  Fig. 9-(c) displays the CSDs for WP squirmers for ξ = 10 and r c = 2.5σ.In this case, hydrodynamics does not significantly affect the CSD, since the CSDs are all monotonically decreasing.Although β = 3 is slightly wider and β = −3 is slightly narrower than the others, the distributions are similar among them, with approximately the same width.Since activity dominates over the interaction strength, we can infer that this independence on β is due to two main effects: on the one hand, the anisotropy of the potential and on the other hand, the range of the interaction.Since in Fig. 9-(d) for WP squirmers and ξ = 10 and r c = 1.5σ, the CSDs have indeed a different width depending on the hydrodynamic signature, corresponding β = 0.5 to the widest distribution.When activity dominates over the interaction strength, we can infer that CSDs depend mainly on hydrodynamics, but the interaction range plays an important role too, mainly for AP and WP pullers.

D. Radius of gyration
The next feature we compute to study the morphology of the clusters is their radius of gyration Rg(s) (Eq.8), as a function of the cluster size s.Fig. 10-(a) displays the radius of gyration R g (s) of AP squirmers with ξ = 1 and r c = 2.5σ.We observe that R g ∼ s 0.64 for β < 3, while for β = 3 (red squares in Fig. 10-(a)), R g starts following the same power law before a crossover to a lower value of R g .It is insightful to compare these results with the CSDs, since AP squirmers with β < 3 have monomodal CSDs, with the width of the distribution depending on the value of β.Nonetheless, Fig. 10-(a) shows that clusters morphology is independent of the hydrodynamic signature when β < 3, but once β = 3 the CSD is bimodal and it turns out that the change in R g coincides with the changes in the CSD.Moreover, dynamic clusters are observed when β < 3 while less dynamic ones are formed when β = 3: therefore, attractive interaction is enhanced when particles are strong pullers.
Similarly, R g ∼ s 0.64 in Fig. 10-(b), when the interaction range is short and for small clusters (s < 50); while for larger clusters of pullers (s > 50) the R g power law is smaller than 0.64.The slower growth of R g depends on the fact that the clusters' structures are mostly controlled by hydrodynamics.For WP squirmers, R g shows a richer behavior: particles first aggregate in trimers and then organize in chain-like structures.As in Fig. 10-(c), R g presents a crossover to a different asymptotic regime.
Chains are observed only for WP independently of β; hence, chain formation is driven by the anisotropic interaction and not by hydrodynamics, and is rather insensitive to the range of the potential, as in Fig. 10-(d).
When self-propulsion dominates over attraction, R g = ks 0.64 for small clusters of AP squirmers independently on β, while for large clusters of pullers R g has a smaller exponent Fig. 11-(a).The same happens when the attraction range is 1.5σ, Fig. 11-(b).In Fig. 11-(c) we report the radius of gyration for the WP squirmers when r c = 2.5σ.Given that activity dominates over attraction, particles do not assembly as chains: R g does not present the cross-over as in Fig. 10-(c), but rather R g ∼ s 0.64 for all s.Moreover, the clusters' morphology is the same despite of the active stress β.However, when the interaction range is reduced to r c = 1.5σ, as in Fig. 11-(d), clusters of pullers are more compact when larger than 20 particles.

E. Local polar and nematic order
We now analyze the degree of alignment as a function of the cluster size, computing in Fig. 12 the polar order P (s) when ξ = 1 using Eq. 9.In general, the polar order inside a cluster usually decreases with cluster size as a power law [26].As reported in Fig. 12-(a), increasing |β| leads to a faster decay of the local polar order for AP squirmers with ξ = 1 and r c = 2.5σ.When decreasing the interaction range, we observe a stronger dependence of the polar order with the hydrodynamic signature: clusters of neutral squirmers (blue triangles in Fig. 12-(b)) are prac-tically swimming in the same direction, while clusters of weak-pullers (black circles in Fig. 12-(b)) show different behavior depending on their size: for small and medium size (up to 100 particles) the polar order decays with size monotonically; the polar order has a crossover for larger clusters, decreasing faster with the cluster size.In general, weak pullers present higher degree of polar order than pushers and strong pullers behave in the same way as the corresponding pushers, (with P (s) decaying with the cluster-size monotonically).
However, the results of the local polar order for WP clusters are quite different.Given that squirmers form elongated clusters from the assembly of trimers, the local order rapidly vanishes as the cluster size grows (Fig. 12-(c)).This tendency is less marked for short range anisotropic interactions and large β (strong pullers in Fig. 12-(d) (red squares)).
In Fig. 13 we report the polar order as a function of the cluster size for both AP and WP squirmers when self-propulsion dominates over attraction (ξ = 10).For AP squirmers, P (s) decays faster for pushers than for pullers, independently on the interaction range (panels a and b of Fig. 13).The degree of polar order is more sensitive to β for pullers than for pushers; weak puller's clusters are more aligned than clusters of stronger pullers, while neutral squirmer's clusters are the most aligned ones for all cluster sizes.Panels c and d of Fig. 13 report the polar order of clusters of WP squirmers, with a monotonous decay of the polar order with cluster size, weakly dependent on β and on the interaction range (except for β = 0.5).
In Fig. 14 we report results for the nematic order λ(s).The local nematic order for AP squirmers with ξ = 1 decays with the cluster size s independently on β for β < 3 while presents an even faster decay when β = 3 (see Fig. 14-(a)).When r c = 1.5σ, λ(s) decays with s in the same way as P (s) does for all β.Neutral squirmers are completely polarized (as also shown in the nematic order, blue triangles in Fig. 14-(b)).The high polar order observed in clusters of puller and neutral squirmers is also detected by the nematic order.
The behaviour of the local nematic order for WP clusters with ξ = 1 is quite striking.The reason why λ(1) = λ(2) = 1 is that when the interaction range is 2.5σ, once two particles collide they are kept together by the attractive patch.Particles in trimers and tetramers point towards the center of the cluster and are kept together by the attractive patch, reason why they have a small nematic order (all curves in Fig. 14-(c)).For larger clusters, active stresses become more important: hydrodynamics attract trimers so that particles' attractive patches reorient in the clusters forming aligned chains.This does not happen for strong pullers, where the nematic order remains low and constant (red squares in Fig. 14-(c)), due to the fact that trimers form disordered chains (made of trimers formed by hydrodynamics).When the interaction range is 1.5σ, the local nematic order behaves in the same way as for 2.5σ, but in this case oriented chains are observed only for strong pushers, whereas the rest of squirmers form less oriented chains with a value of λ(s) in between.
In Fig. 15 we report the local nematic order λ(s) for ξ = 10 and AP squirmers panels a and b.In clusters made of neutral squirmers (blue triangles), particles swim in the same direction independently on the interaction range: λ(s) ≈ 1 for all cluster sizes.The nematic order for all other squirmers decays in the same way as the polar order.To conclude, in the case of clusters of AP squirmers, alignment is not affected by the interaction range.
WP clusters with ξ = 10 and long range interactions develop a high nematic order for neutral and weak pushers (β = 0 and −0.5) as shown in panel c of Fig. 15.This nematic order is not due to a global polar order, (being very low) but to the competition between the anisotropic potential and the hydrodynamic signature: even though activity dominates over attraction, the interaction range is long enough to allow for the development of a nematic order in the system.Undoubtedly, the interaction range is important to develop clusters with nematic order, since λ(s) for WP squirmers with short range has a different behavior, where the nematic order only shows a non-zero value due to the polar order shown before.

F. Global polar and nematic order
To study the orientational order of the system, we compute as in [10] the global polar order parameter and the nematic order tensor Eq. 10 [55], at steady-state at long times.In Fig. 16-(a) we represent the global polar order P ∞ for AP squirmers in steady state for ξ = 1 and ξ = 10 and both interaction ranges.
When ξ = 10, activity dominates and the polar order parameter has a non-zero value for weak pullers (0 ≥ β ≥ 1), similarly to the 3D case without attractive interactions [10].When ξ = 1, the range of interaction affects to the formation of aligned states: for 2.5σ the system is completely isotropic despite of the value of β (red circles in Fig. 16-(a)), whereas for r c = 1.5σ the system still has polar order for weak pullers (black circles), as with ξ = 10.The Fig. 16-(b), shows that the nematic order λ ∞ follows the polar order: for AP squirmers there is not an additional orientation besides the polar order.
In Fig. 16-(c) we represent the global polar order P ∞ for WP squirmers in steady state for ξ = 1 and ξ = 10 and both interaction ranges.For strongly anisotropic interactions (ξ = 1) the suspension does not develop any significant degree of polar ordering.For weakly interacting WP squirmers (ξ = 10) and r c = 1.5σ, the polar order emerges in the same range of β's (black squares) as the AP squirmers.In Fig. 16-(d) we report the nematic order parameter λ ∞ .When ξ = 1 WP squirmers do not develop any nematic order (red and black circles): either polar or nematic order are absent for WP squirmers with strongly anisotropic interactions.A non-zero nematic order is observed for ξ = 10, when r c = 2.5σ and r c = 1.5σ: in the latter, the nematic order is due to the polar order (black squares), while in the former midrange interactions create a new orientational order (red squares).Interestingly, the nematic order parameter for WP squirmers with weak interactions shows a non-zero nematic order when the polar order parameter is zero.

V. DISCUSSION AND CONCLUSIONS
In this paper, we report a numerical study of the selfassembly of Janus squirmers in a quasi two dimensional semi-dilute suspension.
Even though we are aware of the fact that experimental catalytic active colloids might present other types of interactions, like phoretic ones, the nature of such interactions is going to depend on the particular system under study [33].As in the case of active dipolar Janus particles [36] where electrophoresis is effectivelly modelled by the interaction between two imbalanced charges on each particle, in our work we suggest that amphiphilic active Janus colloids [35] can be modelled by means of an effective anisotropic interaction.For this 2D system, we study the effect of the hydrodynamic interactions in the self-assembly of such particles, hopefully shedding some light on experiments with active Janus colloids where the fluid flow field can be measured [39] or in experiments where the self-propulsion direction can be tuned [38].
In a previous work, we have studied aggregation in a semi-dilute quasi two dimensional suspension of isotropically attractive squirmers [26], demonstrating that both clusters' morphology and alignment within clusters depended on particles' interaction strength and hydrodynamic signature.In the current work, by changing the symmetry of the inter-particle interaction (from isotropic to anisotropic) we have observed that the interaction anisotropy, modulated by its range and strength, compete with the active stress to re-orient the particles and therefore, giving rise to a rich aggregation phenomenology.
The study of two different interaction ranges was originated from previous results [22,23,25,26,58] where r c = 2.5σ has been used for isotropic potentials, while in [37,59,60] an interaction range of r c = 1.5σ is used for anisotropic potentials.Both cutoff values are typical in the context of attractive Brownian colloids either passive or active.
As shown in all studied cases, the amount of activity versus attraction has a strong effect on the aggregation dynamics, affecting not only the characteristic cluster size s , but also the time scale at which the steady state (if any) is reached.
When the attractive interaction dominates over propulsion, e.g. at ξ = 0.1, AP squirmers coarsen and the velocity of coarsening will depend on the hydrodynamic stresses for both interaction ranges, similarly to the isotropic case with lower interaction strength (ξ = 0.7) [26], where the system coarsens for pullers and weak pullers and form clusters for strong pushers.As opposed to AP, WP squirmers develop more strinking effects, for long range interactions squirmers coarsen very slowly as far as our calculation time is concerned.While for short range, WP squirmers have reached a clustering steady state, where clusters are self-assembly from small chains or trimers.WP cases at ξ = 0.1 are explained deeply in appendix D.
AP squirmers display a behavior closer to that of squirmers interacting with an isotropic potential when interaction competes with activity and when activity dominates over attraction [26].The competition between interaction and activity always in terms of ξ.For example, the fluctuation of the mean cluster size for weak pullers observed in panel b of Fig. 4 and panels a and b of Fig. 5 or the power law behavior of the CSD for weak pullers in panel b of Fig. 8 and panels a and b of Fig. 9.In contrast, WP squirmers have a stronger tendency to develop novel collective behavior, from polydisperse suspension of monomers, dimers, trimers and even chains of tens of particles.
Our interpretation here is that, since two AP squirmers interact when their attractive patches are within their range, this effect is reduced due to their propulsion in opposite directions.Thus the contribution to the torque is mostly given by hydrodynamic stresses rather than by the anisotropic interactions, as in the isotropic case.
In order to perform a statistical study of the clusters, we focus on ξ = 1 and ξ = 10 when the system is in steady state.We first compute the cluster size distribution (CSD).On the one side, CSD for AP squirmers can be described by a power law with an exponential tail.When ξ = 1, CSDs are characterized by a cutoffalgebraic shape except for strong pullers, which develop a bimodal distribution not present for short interaction range (1.5σ).In this case, the CSD's width follows a behaviour reminding that of isotropic squirmers with ξ > 1, i.e. wider CSDs for weak pullers.When ξ = 10, the effect of the interaction range on the CSDs is not observed for any value of β.On the other side, the CSDs for WP squirmers has a peak for trimers when ξ = 1.The shape of the CSDs is approximately the same independently on the attraction range.When ξ = 10, even though the CSD is described by a power law with an exponential tail, when r c = 2.5σ the larger the value of β the wider the distribution; whereas when r c = 1.5σ the shape of the CSD is mainly dictated by hydrodynamics, with weak pullers showing a wider distribution.
Next, we compute the radius of gyration as a function of the cluster size.When ξ = 1 and r c = 2.5σ, the morphology of AP squirmers follows the same behaviour as the one detected in the isotropic case [26], with more compact clusters for strong pullers; differently from the case when r c = 1.5σ,where clusters are more compact for weak pullers.When ξ = 10 the clusters' R g does not change, independently on the interaction range or the hydrodynamic signature.When ξ = 1, the R g for WP squirmers varies with the clusters sizes since, for both interaction ranges, trimers and tetramers gives rise to the formation of chains, which lead to different exponents.In contrast, when ξ = 10 and r c = 2.5σ, the R g for WP squirmers is the same despite β; however, when r c = 1.5σ pullers form more compact clusters.
It is worth mentioning that these chains resemble the chains observed with amphiphilic active brownian particles [37], but clearly their nature is different since WP squirmers form chains due to their stresslet, while the Brownian case depends on the activity and the size of the hydrophobic patch (see appendix B).Moreover, it seems that hydrodynamic interactions cancel the formation of chains for AP squirmers, since we have observed chains just for WP squirmers.
To establish a potential alignment within clusters, we compute both polar and nematic order.When ξ = 1, the polar order of the clusters for AP squirmers strongly depends on r c : when r c = 2.5σ, P (s) quickly decays with the cluster size (quicker for pushers than for pullers); whereas when r c = 1.5σ,P (s) corresponds to high polar order for large clusters with β = 0 and weak pullers.The same features can be observed when ξ = 10, independently on r c .P (s) for WP squirmers vanishes for all cluster sizes, independently of β and r c when ξ = 1, since the formation of trimers and chains assemble the particles head to head, thus it avoids the polar order.However, when ξ = 10, the polar order decays with s: for r c = 1.5σ weak pullers form clusters with high polar order, while for r c = 2.5σ P (s) decays faster for pushers than for pullers.
In the latter case, the WP squirmers show another type of alignment, with particles with a high nematic order in the range of β = 0 and weak pushers.
During aggregation, we have identified seven different cases depending on whether the attractive patch is oriented against the propulsion direction (AP squirmers) or directed towards it (WP squirmers): 3 gas states, 3 clustering cases and 1 coarsening.On the one side, AP squirmers coarsen isotropically if the interaction is strong enough ξ ≈ 0.1.When the interaction strength competes with self-propulsion (ξ = 1), dynamic clusters emerge whose mean size depend on hydrodynamic stresses.When activity dominates (ξ = 10), particles form a gas and depending on the value of β this gas can be isotropically oriented or polarly oriented.On the other side, WP squirmers form chains of particles if the interaction is high enough (ξ ≈ 0.1), or a suspension of trimers and tetramers if the interaction competes with the active propulsion (ξ = 1).Whenever activity dominates (ξ = 10), particles form a gas and depending on the value of β this gas can be isotropic, polar or even nematically oriented.
Therefore, on one hand anisotropy drives the formation of structures of trimers and tetramers for WP squirmers; on the other hand anisotropy modulates the sensitivity of the hydrodynamic signature: while WP squirmers are more sensitive to r c when ξ = 10, AP squirmers are more sensitive to the interaction range (r c ) when ξ = 1.To conclude, the rich morphology of the detected squirmers' aggregates is the result of the anisotropic interactions, characterized by the angular attractive potential and its interaction range, in competition with the active stress, that can be pointing towards or against the attractive patch.

FIG. 1 .
FIG. 1. Panel a) WP Janus swimmers.Panel b) AP Janus swimmers.The attractive patch is represented in blue (dark gray), and the non-attractive patch is in green (light gray).The swimming direction ê is the red arrow and direction of the attractive patch p is the green one.The center-to-center distance between two swimmers rij is the black vector.

FIG. 2 .
FIG. 2. XY projection of a squirmer at β = −1 (pusher, left panel), β = 0 (neutral, middle panel) and β = 1 (puller, right panel).The red arrow represents the swimming direction, the black arrows the velocity directions of the fluid, the color code the speed of the fluid normalized by the Stokes velocity (vs = 2/3 B1).Velocities are computed in the squirmer frame of reference.

Fig. 4
represents the mean cluster size as a function of time for suspensions where attraction competes with propulsion (ξ = 1) of AP squirmers (interacting via longrange (panel-a) or short-range (panel-b) attraction) and of WP squirmers (interacting via long-range (panel-c) or short-range (panel-d) attraction).

Fig. 5
represents the mean cluster size as a function of time for suspensions where propulsion dominates over attraction (ξ = 10) of AP squirmers (interacting via longrange (panel-a) or short-range (panel-b) attraction) and of WP squirmers (interacting via long-range (panel-c) or short-range (panel-d) attraction).

FIG. 7 .
Fig.7displays the different states, types of collective motion or self-assembly that WP and AP Janus squirmers form as β and ξ vary.This way to classified our simulations, allow us to observe in a better way all the results.We have identified seven different cases: 3 gas systems, 3 clustering systems and 1 coarsening state.All the different cases, show different types of cluster morphologies and/or particles' alignment.The emerging structures result to be sensitive to the patch direction (WP or AP), the hydrodynamic signature β and the interaction range r c and strength ξ.The three gas systems are characterized by the different particles' alignment: isotropic (with particles randomly swimming), polar (with most of the particles swimming in the same direction) and nematic (with approximately half of the total amount of particles swimming in one direction and the other half in the opposite direction).The three clustering systems can be classified as: finite-size dynamic clusters, chains and trimers.In the coarsening case all particles form a unique macroscopic aggregate.AP Janus squirmers coarsen into a single cluster for ξ ≈ 0.1, (see Fig. 7-(a)) and self-organise into finite size clusters for ξ ≈ 1 (see Fig. 7-(b)) (larger values of β

Fig. 8 -
Fig.8-(a) represents the CSDs of AP squirmers for ξ = 1 and r c = 2.5σ.The characteristic distribution width grows with β until β < 3.For all the explored parameters, monomers always represent the largest contribution.For large enough values of β, β = 3, a second peak emerges for cluster sizes involving a few tens of particles (around 25 for the system parameters explored).

Fig. 9 -
Fig.9-(b) displays AP squirmers' CSDs for ξ = 10 and r c = 1.5σ.The CSDs are essentially the same as in the previous case with longer interaction range for pushers and β = 0, while for pullers the CSDs are lightly wider.Fig.9-(c)displays the CSDs for WP squirmers for ξ = 10 and r c = 2.5σ.In this case, hydrodynamics does not significantly affect the CSD, since the CSDs are all monotonically decreasing.Although β = 3 is slightly wider and β = −3 is slightly narrower than the others, the distributions are similar among them, with approximately the same width.Since activity dominates over the interaction strength, we can infer that this independence on β is due to two main effects: on the one hand, the anisotropy of the potential and on the other hand, the range of the interaction.Since in Fig.9-(d)for WP squirmers and ξ = 10 and r c = 1.5σ, the CSDs have indeed a different width depending on the hydrodynamic signature, corresponding β = 0.5 to the widest distribution.When activity dominates over the interaction strength, we can infer that CSDs depend mainly on hydrodynamics, but the interaction range plays an important role too, mainly for AP and WP pullers.

FIG. 10 .
FIG. 10.Radius of gyration normalized by the particle radius as a function of the cluster size for suspensions with ξ = 1.Top panels are AP squirmers (a) with rc = 2.5σ and (b) rc = 1.5σ, whereas bottom panels are WP squirmers (c) with rc = 2.5σ and (d) rc = 1.5σ.Grey dashed lines in all panels are the function Rg = Ks 0.64 .

FIG. 11 .
FIG. 11.Radius of gyration normalized by the particle radius as a function of the cluster size for suspensions with ξ = 10.Top panels are AP squirmers (a) with rc = 2.5σ and (b) rc = 1.5σ, whereas bottom panels are WP squirmers (c) with rc = 2.5σ and (d) rc = 1.5σ.Grey dashed lines in all panels are the function Rg = Ks 0.64 and most of the curves fit with this power law.