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

TL;DR

This study investigates how anisotropic interactions and activity influence the self-assembly and morphology of Janus swimmer clusters, revealing diverse aggregation behaviors in semi-dilute suspensions with tunable propulsion directions.

Contribution

It introduces a systematic analysis of how anisotropic interactions and activity directionality affect cluster formation in Janus swimmers, advancing understanding of active material self-assembly.

Findings

Anisotropic interactions and activity directionality lead to diverse aggregation patterns.

Tuning hydrophobic strength influences cluster morphology.

Active stress competes with attractive interactions to determine structures.

Abstract

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.

Figures

Equations

$$V(r_{ij},{\theta }_i,{\theta }_j)=V_{rep}(r_{ij})+V_{att}(r_{ij})\varphi ({\theta }_i,{\theta }_j)$$

$$\begin{array}{rl}& V_{rep}\left(r_s\right)=\\ & \{\begin{array}{ll}{ϵ}_s\left[{\left(\frac{{\sigma }_s}{r_s}\right)}^{\nu }-{\left(\frac{{\sigma }_s}{h_0}\right)}^{\nu }\left(1-\frac{(r_s-h_0)\nu }{h_0}\right)\right]& 0<r_s\le h_0\\ 0& otherwise,\end{array}\end{array}$$

$$\begin{array}{rl}& V_{att}(r_{ij})=\\ & \{\begin{array}{ll}4ϵ\left[{\left(\frac{\sigma }{r_{ij}}\right)}^{12}-{\left(\frac{\sigma }{r_{ij}}\right)}^6\right]& (2^{1/6}\sigma )<r_{ij}<r_c,\\ -ϵ& r_{ij}<(2^{1/6}\sigma ),\\ 0& otherwise,\end{array}\end{array}$$

$$V_{att}(r_{ij})=4ϵ\left[{\left(\frac{\sigma }{r_{ij}}\right)}^{12}-{\left(\frac{\sigma }{r_{ij}}\right)}^6\right]+4ϵ\left[6{\left(\frac{\sigma }{r_c}\right)}^{12}-3{\left(\frac{\sigma }{r_c}\right)}^6\right]{\left(\frac{r_{ij}}{r_c}\right)}^2-4ϵ\left[7{\left(\frac{\sigma }{r_c}\right)}^{12}+4{\left(\frac{\sigma }{r_c}\right)}^6\right],$$

$$\displaystyle\begin{split}&\phi(\theta_{i},\theta_{j})=&\\ &\begin{cases}\cos\theta_{i}\cos\theta_{j}&\left(\hat{p}_{i}\cdot\hat{\mu}_{ji}\right)>0\cap\left(\hat{p}_{j}\cdot\hat{\mu}_{ij}\right)>0\,,\\ 0&otherwise\,.\\ \end{cases}\end{split}$$

$$\text{u}{|}_{R_p}=\left[B_1\sin \theta +B_2\sin \theta \cos \theta \right]\mathbf{\tau },$$

$$\xi =\frac{F_d}{F_{att}},$$

$$R_g(s)=\sqrt{\sum _{i,j=1}^s\frac{{\left({\vec{r}}_i-{\vec{r}}_j\right)}^2}{2s}},$$

$$P(s)=\left|\frac{{\sum }_{i=1}^s{\hat{e}}_i}{s}\right|$$

$$Q_{hk}(t)=\frac{1}{N}\sum _{i=1}^N\left(2e_{ih}(t)e_{ik}(t)-{\delta }_{hk}\right),$$

$$\frac{f(s)}{f(1)}=A\frac{\exp (-(s-1)/s_0)}{s^{{\gamma }_0}}+B\frac{\exp (-(s-1)/z_0)}{s^{-{\gamma }_0}},$$

$$\displaystyle\begin{split}&\phi_{a}\left(\theta_{i}\right)=\\ &\left\{\begin{array}[]{lc}1&\mbox{$\theta_{i}\leq\theta_{max}$}\\ \cos^{2}\left(\frac{\pi\left(\theta_{i}-\theta_{max}\right)}{2\theta_{tail}}\right)&\mbox{$\theta_{max}\leq\theta_{i}\leq\theta_{max}+\theta_{tail}$}\\ 0&\mbox{otherwise.}\end{array}\right.\end{split}$$

Full text

I Introduction

In the last decade, the pioneer experiments by Dombrowski and Cisneros Dombrowski et al. (2004); Cisneros et al. (2007) inspired the numerical work aimed at understanding the influence of hydrodynamics in the formation of coherent structures of micro-organisms, simulated as spherical squirmers Ishikawa et al. (2006) in 3D ISHIKAWA and PEDLEY (2007); ISHIKAWA et al. (2008) and quasi-2D Ishikawa and Pedley (2008). The features of a 3D squirmer suspension have been numerically characterized using the smoothed profile method Molina et al. (2013); Oyama et al. (2016, 2017) and Lattice Boltzmann (LB), either in bulk Alarcón and Pagonabarraga (2013) or in response to a steady Couette flow Pagonabarraga and Llopis (2013). While Evans et al. Evans et al. (2011) used Stokesian Dynamics to detect the instability of the isotropic state in a 3D suspension of weak pullers, Delmotte and coworkers Delmotte et al. (2015) 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. Palacci et al. (2013); Schwarz-Linek et al. (2012); Petroff et al. (2015); Theurkauff et al. (2012); Thutupalli et al. (2011) have inspired numerical simulations of 2D and quasi-2D suspensions of squirmers interacting via a repulsive potential. When a suspension (simulated by means of multi-particle collision dynamics) was confined between parallel walls (quasi-2D geometry) Zöttl and Stark (2014), 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 Blaschke et al. (2016), differently to the 2D suspension of disk-like squirmers Matas-Navarro et al. (2014), where hydrodynamics suppressed phase separation.

Inspired by the self-assembly observed in experiments with either bacteria Petroff et al. (2015) or active colloids Theurkauff et al. (2012) and by numerical works on 2D dilute suspensions of active brownian particles (ABP) interacting via isotropic short-range attractions Redner et al. (2013); Mognetti et al. (2013); Tung et al. (2016); Navarro and Fielding (2015), we have studied a semi-diluted quasi 2D suspensions of attractive squirmers Alarcón et al. (2017), 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 H2O2 and O2 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 Saha et al. (2014); Pohl and Stark (2014); Meyer et al. (2014); Liebchen et al. (2015, 2016); Nejad and Najafi (2017), eventhough there are experiments, simulations and theories that consider them relevant for particles’ propulsion. On the other side, which interactions dominate in these systems are still a matter of debate Liebchen and Löwen (2019).

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. Manjare et al. (2014) found that particles with an hydrophobic patch would move faster than the ones with an hydrophilic patch. Similarly, Gao et al Gao et al. (2013), 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 Yan et al. (2016) 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 Yan et al. (2016), 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 Gao et al. (2013); Yan et al. (2016), simulations of a dilute quasi-2D suspension of ABP interacting via anisotropic amphiphilic Janus interactions have been reported in Ref. Mallory et al. (2017). 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. Mallory et al. (2017) 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 Brown and Poon (2014).

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 Campbell et al. (2018).

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 Lighthill (1952); Blake (1971); Pedley (2016) interacting via a Janus anisotropic potential similar to the one used in Ref. Mallory et al. (2017) for Active Brownian Particles (ABP). In order to model the Janus swimmers, we have introduced a pair potential $V(r_{ij},\theta_{i},\theta_{j})$ which depends on the distance between two squirmers ($r_{ij}$) and their attractive patches’ orientations

[TABLE]

where $\theta_{i}$ is the angle between the patch unit vector $\vec{p_{i}}$ and the inter-particle vector $\vec{r}_{ji}=\vec{r}_{j}-\vec{r}_{i}$; and $\theta_{j}$ the angle between $\vec{p_{j}}$ and $\vec{r}_{ij}=-\vec{r}_{ji}$.

$V_{rep}(r_{ij})$ represents a short range soft repulsion Allen and Tildesley (2017) (to avoid overlapping) given by

[TABLE]

where $\epsilon_{s}=1.5$ is the energy scale, $\nu=2$, $r_{s}\equiv r_{ij}-\sigma$ is the inter-particle distance (with $r_{s}>0$ to ensure that particles never overlap), $\sigma_{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}$. 111Note that separation $h_{0}$ is relatively short, $h_{0}<\sigma_{s}$, and its value depends on the location of the minimum of the attractive potential.

The attractive potential consists of a radial $V_{att}(r_{ij})$ and an angular $\phi(\theta_{i},\theta_{j}$) contribution. Concerning the radial part of the potential, we have considered two different truncated Lennard-Jones: one long ($r_{c}=2.5\sigma$) and another short ($r_{c}=1.5\sigma$) range. While the following equation:

[TABLE]

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}$ Stoddard and Ford (1973) (see appendix A). Therefore, whenever $(2.2^{1/6}\sigma)<r_{ij}<r_{c}$

[TABLE]

whereas for $r_{ij}<(2.2^{1/6}\sigma)$ $V_{att}=-0.305\epsilon$ and $V_{att}=0$ otherwise.

Inspired by the work of Hong et al. Hong et al. (2008) for amphiphilic colloidal spheres, the anisotropic interaction is controlled by the angular contribution of the attractive potential

[TABLE]

where $\hat{\mu}_{ij}=\vec{r}_{ij}/r_{ij}$ is the unit vector joining the squirmers’ center of mass. This anisotropic potential ensures a non-vanishing torque when $\left(\hat{p}_{i}\cdot\hat{\mu}_{ji}\right)>0$ and $\left(\hat{p}_{j}\cdot\hat{\mu}_{ij}\right)>0$ (a more detailed analysis is discussed in appendix B).

Within our model, each Janus particle self-propels along a direction $\hat{e}_{i}$ rigidly bounded to the particle. As in Ref. Mallory et al. (2017), we study two types of active Janus particles depending of the orientation of $\hat{e}_{i}$ with respect to $\hat{p}_{i}$. When the propulsion direction ($\hat{e}_{i}$, red arrow in figure 1-left panel) is pointing towards the attractive patch $\hat{p}_{i}$ (red arrow in figure 1-left panel), swimmers are propelling ”with the patch” (WP). When the propulsion direction ($\hat{e}_{i}$, red arrow in figure 1-right panel) is pointing against the attractive patch $-\hat{p}_{i}$ (green arrow in figure 1-right panel), swimmers are propelling ”against the patch” (AP).

A spherical Janus squirmer of radius $R_{p}$ is characterized by the tangential surface slip velocity,

[TABLE]

where $\boldsymbol{\tau}$ is a unit vector tangential to the particle surface, $B_{1}$ quantifies the asymptotic self-propelling speed of an isolated squirmer ($v_{s}=\frac{2}{3}B_{1}$) Anderson and Prieve (1991); Stone and Samuel (1996) and $B_{2}$ is proportional to the active stress generated on the surrounding fluid. The ratio between active stress and self-propelling velocity, $\beta=B_{2}/B_{1}$ Ishikawa et al. (2006), quantifies the active state of the squirmer and its interaction with the fluid. A positive value of $\beta$ corresponds to the case when thrust is generated in front of the squirmer’s body (puller); a negative value of $\beta$ to the case when thrust is generated at the back (pusher)  ISHIKAWA and PEDLEY (2007) and $\beta=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 Cates et al. (2005); Succi (2001), implemented in a highly parallelized code Stratford and Pagonabarraga (2008) that exploit the excellent scalability of LB on supercomputing facilities. Spherical squirmers (with a diameter $\sigma$ of 4.6 lattice nodes Nguyen and Ladd (2002)) are individually resolved, imposing a modified bounce-back rule on the one-particle velocity distribution functions to nodes crossing the solid particle’s boundary (including the slip velocity to impose the squirming motion Llopis and Pagonabarraga (2010); Matas Navarro and Pagonabarraga (2010)). 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 $\phi=0.10$ (where $\phi=\frac{\pi}{4}\rho\sigma^{2}$ and $\rho=N/L^{2}$) in a quasi two-dimension geometry with periodic boundary conditions ($L\times L\times k\sigma$, with $k=5$ and $L\approx 83\sigma$) 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\approx 167\sigma$, in order to proof that results are not due to finite size effects, specially cases with large cluster sizes.

As in Ref. Alarcón et al. (2017), we quantify the competition between attractive and self-propelling forces via the dimensionless parameter

[TABLE]

where $F_{d}=6\pi\eta 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}\sigma=1.245\sigma$ for the long-range potential in Eq. 3 and to $r=1.239\sigma$ for the short-range potential in Eq. 4. On the one hand $\xi$ 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 $\tau=\frac{\sigma}{v_{s}}$, we perform simulations from $1450$ up to $3000$ $\tau$. Once the system reaches steady state, in a time window between $1000$ and $2000$ $\tau$, we carry out a systematic analysis of the dynamics of the squirmer suspension, considering $\xi$ from 0.1 to 10 and $\beta$ 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. Alarcón et al. (2017) and study the degree of aggregation as a function of $\beta$ and $\xi$. 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\sigma$ (first minimum of the $g(r)$, see appendix C).

To calculate the cluster-size distribution $f(s)$ we apply a criterion based on Valeriani et al. (2012): (i) We arbitrarily subdivide the range of $s$-values into intervals $\Delta s_{i}=(s_{i,max}-s_{i,min})$, where $n_{i}^{t}$ is the total number of clusters within each interval $\Delta s$; (ii) we assign the value $n_{i}=n^{t}_{i}/\Delta s_{i}$ to every $s$ within $\Delta s_{i}$, and compute the fraction of clusters of size $s$ as $f(s)=n_{i}/N_{c}$, where $N_{c}=\sum_{i}n_{i}\Delta s_{i}$ is the total number of clusters.

To characterise the clusters’ morphology, we compute the radius of gyration

[TABLE]

the average polar order parameter for each cluster size $s$

[TABLE]

the tensor order parameter (for a $2D$ system Barci and Stariolo (2009))

[TABLE]

(being $h$ and $k$ equal to ${x,y}$ and $N$ the total number of squirmers) and the nematic order parameter $\lambda(t)$ being the largest eigenvalue of $Q_{hk}$ Eppenga and Frenkel (1984).

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_{\infty}=P(t>>0)$ while the nematic order in steady-state is $\lambda_{\infty}=\lambda(t>>0)$.

IV Results

IV.1 Mean cluster size

We start with computing the mean cluster size $\langle s\rangle$ as a function of time for different values of $\beta$ and $\xi$. 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 ($\xi=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 (panel-c) 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 $\beta$. Pullers coarsening faster than pushers, a part from the system with $\beta=-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 $\beta$ 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 $\langle s\rangle<10$ even when $t>1000\tau$ (it is worth noting that we have never observed a steady state despite having run up to $6\times 10^{6}$ LB-time steps, corresponding to $\sim 8700\tau$). 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 $\langle s\rangle$ reaches a steady state at long $t$ for all values of $\beta$. A more detailed analysis is given in appendix D for the case of the WP system with $\xi=0.1$.

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

AP squirmers interacting via long-range attractions (panel a) quickly reach a steady state characterized by $\langle s\rangle\approx 2$ except for $\beta=3$ where $\langle s\rangle\approx 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 $\langle s\rangle\approx 2$ as well as pullers with $\beta=3$ (with a slightly larger mean cluster size of $\langle s\rangle\approx 2.5$). Interestingly, the mean cluster size for weak pullers (black curve in panel-b) first develops a metastable state at short times ($t\in(10,500)\tau$) and then another one for longer times ($t>1000\tau$): the latter due to the formation of dynamic clusters with polar order. Similarly, when $\beta=0$ the decay in $\langle s\rangle$ at long time corresponds to monomers aligning and moving in the same direction.

WP squirmers with $\beta<3$ relatively quickly reach steady-state with clusters with an average of 4 particles (panel-c). Whereas pullers with $\beta=3$ reach a slightly higher mean cluster size $\langle s\rangle\sim 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 $\sim 4$ particles, independently on their hydrodynamic signature.

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

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 $\beta$, fluctuating around a value of 3 for long times. In the case of neutral squirmers and large $\xi$ and after a long time ($t>1000\tau$), 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 $\xi$, 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 the mean cluster size as a function of $\beta$ for $r_{c}=2.5\sigma$ (panel a) and $r_{c}=1.5\sigma$ (panel b). 222The WP system with $\xi=0.1$, where clustering is observed (panel d of Fig. 3), is studied in more detail in appendix D.

As shown in figure 6, when $\xi=1$ WP swimmers mostly form tetramers and trimers, while WP swimmers with $\xi=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 $\xi=\{1$,$10$} and the two interaction ranges studied here.

For lower values of $\xi$, the dependence on $\beta$ becomes more evident for WP squirmers with $\xi=0.1$ (gray squares on panel-b): while for extreme values of $\beta$ (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\sigma$ (circular symbols panel-a) show that the mean cluster size increases as $\beta$ increases when $\xi=1$, whereas it does not change significantly when $\xi=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 $\xi=1$ and $10$, the mean cluster size $\langle s\rangle$ for WP squirmers with interaction range of $r_{c}=1.5\sigma$ behaves in the same way as the one with a larger cutoff ($r_{c}=2.5\sigma$). 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 $\langle s\rangle$ for AP squirmers with $\xi=1$ and $r_{c}=1.5\sigma$ (red circles in figure 6-(b)) does not behave in the same way than when $r_{c}=2.5\sigma$ (red cirlces in figure 6-(a)), the former showing a systematic lower value of $\langle s\rangle$ for all $\beta$ while the latter a monotonic increase of the cluster size for pullers. Finally, AP squirmers with $\xi=10$ (black circles in figure 6) show the same behavior varying $\beta$, independently on the attraction range. With a small peak for weak pullers, and a minimum at $\beta=0$.

IV.2 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 $\beta$ and $\xi$ 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 $\beta$ and the interaction range $r_{c}$ and strength $\xi$. 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 $\xi\approx 0.1$, (see Fig. 7-(a)) and self-organise into finite size clusters for $\xi\approx 1$ (see Fig. 7-(b)) (larger values of $\beta$ favoring the formation of larger finite-size clusters, accelerating coarsening). When $\xi=10$ , AP squirmers either form a polar ordered gas (see Fig. 7-(c), for weak positive stresslets, $0<\beta<1$) or an isotropic gas (for AP pushers and pullers with $\beta>1$). WP Janus squirmers have a strong tendency to form small micelles of three or four particles pointing at each other. For $\xi\approx 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 $\xi=1$ the competition with the active stress ($\beta$) gives rise to trimers and tetramers avoiding the formation of larger chains of particles (see Fig. 7-(e)). If activity dominates ($\xi=10$) and $|\beta|<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 $\beta$, we observe an isotropic system. In the Supplemental Material See Supplemental Material at [URL will be inserted by publisher] for additional information and videos 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.

IV.3 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 $\xi=1$ and $\xi=10$ and as a function of $\beta$, the direction of the attractive patch and the interaction range. The case when attraction dominates with respect to self-propulsion ($\xi=0.1$) for AP squirmers is analyzed in appendix D. Results for $\xi=1$ are presented in Fig. 8 and for $\xi=10$ in Fig. 9. Following previous studies Alarcón et al. (2017), where CSD was calculated, we use the same analytical function as a reference.

[TABLE]

with $\gamma_{0}$, $s_{0}$, $z_{0}$ and $B$ constants such that $A=1-B$.

Fig. 8-(a) represents the CSDs of AP squirmers for $\xi=1$ and $r_{c}=2.5\sigma$. The characteristic distribution width grows with $\beta$ until $\beta<3$. For all the explored parameters, monomers always represent the largest contribution. For large enough values of $\beta$, $\beta=3$, a second peak emerges for cluster sizes involving a few tens of particles (around 25 for the system parameters explored). Accordingly, CSDs with $\beta<3$ can be described by Eq. 11, with $B=0$ while the CSD for $\beta=3$ requires a non-vanishing 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 $\beta=3$. At larger values of $\beta$ the behavior changes due to the appearance of a second peak in the CSDs.

Fig. 8-(b) represents the CSDs of AP squirmers for $\xi=1$ and $r_{c}=1.5\sigma$. 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 $\beta=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\sigma$.

Fig. 8-(c) represents the CSDs of WP squirmers for $\xi=1$ and $r_{c}=2.5\sigma$. A peak appears between trimers and tetramers, while the characteristic CSD width increases with the active stress $\beta$. Interestingly, when $\beta<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 $\xi=1$ and $r_{c}=1.5\sigma$, 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 $|\beta|$, the wider the distribution (with distributions wider than the ones with longer interactions range $r_{c}=2.5\sigma$, and with a higher number of monomers especially when $\beta=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 $\beta$, 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 $\xi=10$ and $r_{c}=2.5\sigma$. The corresponding distributions are monotonically decreasing, particularly for $\beta=0.5$ a power law distribution emerges while in the other cases, the CSD can be approximated by Eq. 11 with $B\neq 0$.

Fig. 9-(b) displays AP squirmers’ CSDs for $\xi=10$ and $r_{c}=1.5\sigma$. The CSDs are essentially the same as in the previous case with longer interaction range for pushers and $\beta=0$, while for pullers the CSDs are lightly wider.

Fig. 9-(c) displays the CSDs for WP squirmers for $\xi=10$ and $r_{c}=2.5\sigma$. In this case, hydrodynamics does not significantly affect the CSD, since the CSDs are all monotonically decreasing. Although $\beta=3$ is slightly wider and $\beta=-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 $\beta$ 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 $\xi=10$ and $r_{c}=1.5\sigma$, the CSDs have indeed a different width depending on the hydrodynamic signature, corresponding $\beta=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.

IV.4 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 $\xi=1$ and $r_{c}=2.5\sigma$. We observe that $R_{g}\sim s^{0.64}$ for $\beta<3$, while for $\beta=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 $\beta<3$ have monomodal CSDs, with the width of the distribution depending on the value of $\beta$. Nonetheless, Fig. 10-(a) shows that clusters morphology is independent of the hydrodynamic signature when $\beta<3$, but once $\beta=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 $\beta<3$ while less dynamic ones are formed when $\beta=3$: therefore, attractive interaction is enhanced when particles are strong pullers.

Similarly, $R_{g}\sim 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 $\beta$; 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 $\beta$, 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\sigma$, Fig. 11-(b). In Fig. 11-(c) we report the radius of gyration for the WP squirmers when $r_{c}=2.5\sigma$. 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}\sim s^{0.64}$ for all $s$. Moreover, the clusters’ morphology is the same despite of the active stress $\beta$. However, when the interaction range is reduced to $r_{c}=1.5\sigma$, as in Fig. 11-(d), clusters of pullers are more compact when larger than 20 particles.

IV.5 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 $\xi=1$ using Eq. 9. In general, the polar order inside a cluster usually decreases with cluster size as a power law Alarcón et al. (2017).

As reported in Fig. 12-(a), increasing $|\beta|$ leads to a faster decay of the local polar order for AP squirmers with $\xi=1$ and $r_{c}=2.5\sigma$. 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 practically 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 $\beta$ (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 ($\xi=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 $\beta$ 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 $\beta$ and on the interaction range (except for $\beta=0.5$).

In Fig. 14 we report results for the nematic order $\lambda(s)$.

The local nematic order for AP squirmers with $\xi=1$ decays with the cluster size $s$ independently on $\beta$ for $\beta<3$ while presents an even faster decay when $\beta=3$ (see Fig. 14-(a)). When $r_{c}=1.5\sigma$, $\lambda(s)$ decays with $s$ in the same way as $P(s)$ does for all $\beta$. 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 $\xi=1$ is quite striking. The reason why $\lambda(1)=\lambda(2)=1$ is that when the interaction range is $2.5\sigma$, 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\sigma$, the local nematic order behaves in the same way as for $2.5\sigma$, 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 $\lambda(s)$ in between.

In Fig. 15 we report the local nematic order $\lambda(s)$ for $\xi=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: $\lambda(s)\approx 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 $\xi=10$ and long range interactions develop a high nematic order for neutral and weak pushers ($\beta=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 $\lambda(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.

IV.6 Global polar and nematic order

To study the orientational order of the system, we compute as in Alarcón and Pagonabarraga (2013) the global polar order parameter and the nematic order tensor Eq. 10 Barci and Stariolo (2009), at steady-state at long times. In Fig. 16-(a) we represent the global polar order $P_{\infty}$ for AP squirmers in steady state for $\xi=1$ and $\xi=10$ and both interaction ranges.

When $\xi=10$, activity dominates and the polar order parameter has a non-zero value for weak pullers ($0\geq\beta\geq 1$), similarly to the 3D case without attractive interactions Alarcón and Pagonabarraga (2013). When $\xi=1$, the range of interaction affects to the formation of aligned states: for $2.5\sigma$ the system is completely isotropic despite of the value of $\beta$ (red circles in Fig. 16-(a)), whereas for $r_{c}=1.5\sigma$ the system still has polar order for weak pullers (black circles), as with $\xi=10$. The Fig. 16-(b), shows that the nematic order $\lambda_{\infty}$ 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_{\infty}$ for WP squirmers in steady state for $\xi=1$ and $\xi=10$ and both interaction ranges. For strongly anisotropic interactions ($\xi=1$) the suspension does not develop any significant degree of polar ordering. For weakly interacting WP squirmers ($\xi=10$) and $r_{c}=1.5\sigma$, the polar order emerges in the same range of $\beta$’s (black squares) as the AP squirmers. In Fig. 16-(d) we report the nematic order parameter $\lambda_{\infty}$. When $\xi=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 $\xi=10$, when $r_{c}=2.5\sigma$ and $r_{c}=1.5\sigma$: in the latter, the nematic order is due to the polar order (black squares), while in the former mid-range 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 self-assembly 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 Liebchen and Löwen (2019). As in the case of active dipolar Janus particles Yan et al. (2016) 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 Gao et al. (2013) 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 Campbell et al. (2018) or in experiments where the self-propulsion direction can be tuned Brown and Poon (2014).

In a previous work, we have studied aggregation in a semi-dilute quasi two dimensional suspension of isotropically attractive squirmers Alarcón et al. (2017), 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 Alarcón et al. (2017); Navarro and Fielding (2015); Redner et al. (2013); Mognetti et al. (2013); Prymidis et al. (2015) where $r_{c}=2.5\sigma$ has been used for isotropic potentials, while in Mallory et al. (2017); Miller and Cacciuto (2009); Pu et al. (2017) an interaction range of $r_{c}=1.5\sigma$ 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 $\langle s\rangle$, but also the time scale at which the steady state (if any) is reached.

When the attractive interaction dominates over propulsion, e.g. at $\xi=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 ($\xi=0.7$) Alarcón et al. (2017), 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 $\xi=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 Alarcón et al. (2017).

The competition between interaction and activity always in terms of $\xi$.

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 $\xi=1$ and $\xi=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 $\xi=1$, CSDs are characterized by a cutoff-algebraic shape except for strong pullers, which develop a bimodal distribution not present for short interaction range ($1.5\sigma$). In this case, the CSD’s width follows a behaviour reminding that of isotropic squirmers with $\xi>1$, i.e. wider CSDs for weak pullers. When $\xi=10$, the effect of the interaction range on the CSDs is not observed for any value of $\beta$. On the other side, the CSDs for WP squirmers has a peak for trimers when $\xi=1$. The shape of the CSDs is approximately the same independently on the attraction range. When $\xi=10$, even though the CSD is described by a power law with an exponential tail, when $r_{c}=2.5\sigma$ the larger the value of $\beta$ the wider the distribution; whereas when $r_{c}=1.5\sigma$ 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 $\xi=1$ and $r_{c}=2.5\sigma$, the morphology of AP squirmers follows the same behaviour as the one detected in the isotropic case Alarcón et al. (2017), with more compact clusters for strong pullers; differently from the case when $r_{c}=1.5\sigma$, where clusters are more compact for weak pullers. When $\xi=10$ the clusters’ $R_{g}$ does not change, independently on the interaction range or the hydrodynamic signature. When $\xi=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 $\xi=10$ and $r_{c}=2.5\sigma$, the $R_{g}$ for WP squirmers is the same despite $\beta$; however, when $r_{c}=1.5\sigma$ pullers form more compact clusters.

It is worth mentioning that these chains resemble the chains observed with amphiphilic active brownian particles Mallory et al. (2017), 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 $\xi=1$, the polar order of the clusters for AP squirmers strongly depends on $r_{c}$: when $r_{c}=2.5\sigma$, $P(s)$ quickly decays with the cluster size (quicker for pushers than for pullers); whereas when $r_{c}=1.5\sigma$, $P(s)$ corresponds to high polar order for large clusters with $\beta=0$ and weak pullers. The same features can be observed when $\xi=10$, independently on $r_{c}$. $P(s)$ for WP squirmers vanishes for all cluster sizes, independently of $\beta$ and $r_{c}$ when $\xi=1$, since the formation of trimers and chains assemble the particles head to head, thus it avoids the polar order. However, when $\xi=10$, the polar order decays with $s$: for $r_{c}=1.5\sigma$ weak pullers form clusters with high polar order, while for $r_{c}=2.5\sigma$ $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 $\beta=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 $\xi\approx 0.1$. When the interaction strength competes with self-propulsion ($\xi=1$), dynamic clusters emerge whose mean size depend on hydrodynamic stresses. When activity dominates ($\xi=10$), particles form a gas and depending on the value of $\beta$ 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 ($\xi\approx 0.1$), or a suspension of trimers and tetramers if the interaction competes with the active propulsion ($\xi=1$). Whenever activity dominates ($\xi=10$), particles form a gas and depending on the value of $\beta$ 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 $\xi=10$, AP squirmers are more sensitive to the interaction range ($r_{c}$) when $\xi=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.

Acknowledgments

We thank to S.A. Mallory and A. Cacciuto for helpful discussions. This work was possible thanks to the access to MareNostrum Supercomputer at Barcelona Supercomputing Center (BSC). This article is based upon work from COST Action MP1305, supported by COST (European Cooperation in Science and Technology). This work was supported by FIS2016-78847-P of the MINECO and the UCM/ Santander PR26/16-10B-2. IP acknowledges MINECO and DURSI for financial support under projects FIS2015-67837- P and 2017 SGR-884, respectively. FA acknowledges funding from Juan de la Cierva-formación program.

Appendix A Attractive part of the potential

Fig. 17 compares the radial part of the total short range potential given by Eq. 4 (red curves) and the classical LJ potential Eq. 3 and cut-off distance of $2.5\sigma$ (blue curves), $\sigma$ is the particles’ diameter. Both curves in Fig 17 are the sum of the repulsive soft-sphere potential (Eq. 2) and the respective attractive potential.

Appendix B Different types of Janus interactions

There are different models to simulate the amphiphilic behavior of Janus colloids, in this appendix we have compared two different models for the angular potential $\phi$ in equation 1. The angular potential used in this manuscript (equation 5) has been used to simulate spherical particles with one hydrophobic hemisphere and charged on the other Hong et al. (2008) where both hemispheres have the same size, thus lets call it symmetric potential $\phi_{s}(\theta_{i},\theta_{j})$.

On the other hand, in reference Miller and Cacciuto (2009), they developed an asymmetric potential, where the size of the hydrophobic region is tuned, in order to study 3D Janus colloid suspensions in equilibrium and more recently, some of us have used this potential to study the self-assembly of Active amphiphilic Janus particles with attractive patches of different sizes Mallory et al. (2017). Analogously to the symmetric potential above, the the orientation depending term is given by

[TABLE]

The orientational part of the potential $\phi_{a}$, is a a smooth step function that modulates the angular dependence of the potential, the smoothness is modulated by the parameter $\theta_{tail}$, which it is set to $\theta_{tail}=0.2618$ in Miller and Cacciuto (2009) for 3D Janus colloids while in Mallory et al. (2017) was tuned in $\theta_{tail}=0.436$ to generate a sufficiently smooth potential at the Janus interface and $\theta_{max}$ tunes the size of the hydrophobic region, for a patch coverage of 50$\%$ then $\theta_{max}=\pi/2$.

By the analysis of the curves shown in Fig. 18, it is clear that interactions among the Janus particles will depend on the angular function we choose, despite of the fact that $\theta_{max}=\pi/2$ in the asymmetric potential.

In terms of the interaction strength, it is clear that $\phi_{s}$ will give us a lower attraction than $\phi_{a}$ for any direction of $\theta_{i}$. An important remark is that $\phi_{s}=0$ when two particles are parallel while $\phi_{a}\neq 0$ in the same configuration, therefore $\phi_{a}$ enhance the emergence of chain-like aggregates Mallory et al. (2017).

To understand the relevance of the orientational part in the potential and the hydrodynamic interactions, we are working in a model of amphiphilic squirmers using $\phi_{a}$ in order to compare with the ’dry’ case of Ref. Mallory et al. (2017), this work is in progress.

Appendix C Identifying clusters in the suspension

In order to determine the criterion to characterize whether a particle belongs or not to a cluster, we computed first the radial distribution function of our swimmer suspensions.

We observe a local peak at $r=1.6\sigma$ for the $g(r)$ corresponding to $\beta=3$ with $r_{c}=2.5\sigma$ with AP interaction. We took $r_{cl}=1.75\sigma$ for all cases, although for this case it is evident the first minimum cannot be located there. We analyze how the mean cluster size and the cluster size distribution vary for this troublesome case with the choice of $r_{cl}$, comparing $r_{cl}=1.75\sigma$ to a cutoff in the peak $r_{cl}=1.6\sigma$ and a cutoff lying in the first valley $r_{cl}=1.4\sigma$.

We observe that as we reduce $r_{cl}$, the mean cluster size is reduced, being the difference between $r_{cl}=1.75\sigma$ and $r_{cl}=1.4\sigma$ of 0.9 units ($\approx 10\%$). As for the cluster size distribution, we find that all three of them are very similar. Even if the frequencies for the largest observed clusters are reduced with $r_{cl}$, all distributions are cut at the same bin (clusters sized between 240 and 318 particles), with a frequency $7.7\cdot 10^{-6}$ for $s=279$ and $r_{cl}=1.75\sigma$, and $4.6\cdot 10^{-6}$ for the same size with $r_{cl}=1.4\sigma$. The monomer frequency goes from 0.680 for $r_{cl}=1.4\sigma$ to 0.656 for $r_{cl}=1.75\sigma$, and all three curves have a very similar shape, sharing even local maxima and minima along the whole profile.

Appendix D Clustering analysis for WP Janus with $\xi=0.1$.

Janus squirmers with interaction strength $\xi=0.1$ develop very interesting dynamics, on one hand, for strong pullers with $\beta=3$, particles form mainly trimers and tetramers that eventually aggregate due to hydrodynamic interactions (panel a of Fig. 21 and Fig. 22). The other case is the chains system, where particles self-assembly in chains of different sizes, such chains are form by two lines of squirmers with the attractive patch pointing towards the other line of squirmers that form the chain (panel b of Fig. 21 and Fig. 22). The chain-length depends on the value of $\beta$. This chains system has been observed for pushers and weak pullers as well.

When interaction range is $r_{c}=2.5\sigma$, the clusters either trimers or chains, keep growing slowly, as far as our calculation capacity could reach, we have observed coarsening for all $\beta$. In contrast with WP squirmers with interaction range of $r_{c}=1.5\sigma$. In this case WP squirmer suspensions reach a clustering state for all hydrodynamic signature $\beta$. We have found two general types of clustering:

We have characterized the clusters observed for these WP squirmers. In (panel a in Fig. 23) we show the cluster size distribution $f(s)$. In general, $f(s)$ has two main peaks, one for trimers and another one in the bin of cluster-size around 10 particles for pushers and weak pullers, while strong pullers (red squares) have just one peak for trimers. Moreover, none of these distributions show monomers and just pullers show dimers in their distributions. Another interesting feature of the cluster size distributions are the fact that both strong pushers and pullers have bigger clusters than the rest of squirmers (gray diamonds and red squares respectively) and both distributions follow the analytical shape of equation (11) (dotted line) for large clusters.

Even though that both cluster size distributions are similar for strong pushers and pullers, the morphology of such clusters are not the same. As we can observe in Fig. 23-panel b, where despite that the radius of gyration for small clusters ($s\leq 10$) follows the same behavior for all $\beta$, the large clusters for strong pullers (red squares) are more compact (small slope) than clusters of the rest of squirmers, in fact, clusters of pushers and weak pullers are chains, thus we can observe a greater slope for large clusters (gray diamonds and pink triangles).

We have also calculated both polar and nematic order parameters at each cluster, we have observed that clusters have no polar order at any value of $\beta$ (see panel c of Fig. 23), but an alignment, given by a non-zero nematic order is observed (see panel d of Fig. 23). For pullers, all dimers are aligned by a head-to-head configuration and then the nematic order is 1, particles in trimers and tetramers in all cases are pointing out to the center of the clusters therefore, the nematic order is zero. But, in the case of pushers to weak pullers all clusters have a large non-zero nematic order, since particles are forming chains of two lines of particles aligned to the center of the chain giving rise to a nematic order. In contrast, for strong pullers (red squares) large clusters are formed by trimers gather together by hydrodynamic interactions, thus particles are pointing in different directions and nematic order parameter is reduced. In fact, the magnitude of nematic order parameter for large clusters is given by the $\beta$ value. The lower the $\beta$ the greater the nematic order.