Self-assembly of active bifunctional Brownian particles

Abstract

In this work, with the intent of exploring the out-of-equilibrium polymerization of active patchy particles in linear chains, we study a suspension of active bifunctional Brownian particles (ABBPs). At all studied temperatures and densities, ABBPs self-assemble in aggregating chains, as opposed to the uniformly space-distributed chains observed in the corresponding passive systems. The main effect of activity, other than inducing chain aggregation, is to reduce the chain length and favour the alignment of the propulsion vectors in the bonding process. At low activities, attraction dominates over activity in the bonding process, causing self-assembly to occur randomly regardless of the particle orientations. Interestingly, we find that at the lowest temperature, as density increases, chains aggregate forming a novel state: MISP, i.e., motility-induced spirals, where spirals are characterised by a finite angular velocity. In contrast, at the highest temperature, density and activity, chains aggregate forming a different novel state (a spinning crystalline cluster) characterised by a compact and hexagonally ordered structure, both translating and rotating. The rotation arises from an effective torque generated by the presence of competing domains where particles self-propel in the same direction.

---

1 Introduction

In the last decades, significant progress in the comprehension of the structural and dynamical properties of liquids has been made through the investigation of colloidal particles interacting via spherically symmetric or (more realistic) anisotropic forces.¹ A practical model proposed to study anisotropic interactions between colloidal particles is the so-called “patchy particle” model, consisting of hard-spheres whose surface is decorated with a finite number of short-range attractive sites.^2,3 Patchy particles have allowed elucidating the behavior of network-forming materials,^4,5 such as water⁶ or silica,⁷ finite aggregates, such as surfactant micelles,⁸ or more complex structures, such as proteins.^9,10 Patchy particles have also represented a novel class of building blocks for constructing precise structures, where the arrangements and the selectivity of the sites dictate the overall structure of the assemblies.^11–14 Indeed, the bottom-up approach of patchy colloidal self-assembly has proven to be pivotal for technological advancements across diverse fields, including materials science,¹⁵ pharmaceutical industry,¹⁶ electronics,¹⁷ nanotechnology,¹⁸ and even food technology.¹⁹

Optimization of colloidal self-assembly, inspired by biological matter, has been further achieved by introducing activity²⁰ on simple colloids, with possible applications ranging from targeted drug delivery²¹ to autonomous depollution of contaminated water and soils.²²

The majority of the published work on active colloidal matter has focused on suspensions of active particles interacting via an isotropic potential (attractive²³ or repulsive²⁴). Only more recently, the field of active matter has branched out to explore the interplay between activity and anisotropic interactions, with the goal of developing a systematic understanding of how active forces can be exploited together with anisotropic forces to design assemblies with desired structural and functional features.^25–27

Active colloids have been shown to aggregate into functional structures not detected in equilibrium systems. When active particles are spherical and repulsive, they undergo a Motility Induced Phase Separation,^24,28 whereas when spherical and attractive they form living clusters.^23,29 Interestingly, when active particles are elongated, they aggregate into functional transient clusters capable of rotating, such as those reported in ref. 30 and 31. Thus, in order to detect a spinning state in a self-assembled suspension of active particles, one needs two main ingredients: particle activity and shape anisotropy. On the other side, one could consider spinning of an already formed structure. Spinning has also been observed in a passive gear embedded in an active bath of elongated particles, such as a bacterial suspension.^32,33 When dealing with active polymers, a spinning spiral state appears whenever the propulsion force along the polymer backbone is tightly parallel to the local tangent.^34,35 A recent work has revealed a spinning state in a suspension of active particles whose shape is more complex than spherical and is characterized by an attractive patch on their surface.³⁶

Full control of the complex dynamics of active patchy colloids remains yet challenging. So far, research has focused on tuning the shape, size, and composition of the patches in order to control autonomous locomotion and spontaneous assembly.^27,36,37 Specific interactions can be obtained by implementing lock and key groups on the particle surfaces, such as DNA oligonucleotides, protein cross-linkers or antibody–antigen binding pairs.³⁸ Due to their ability to self-assemble into chains, sheets, rings, icosahedra, tetrahedra, etc., patchy colloids provide access to a broad range of active colloidal materials.³⁹

In this article, we explore the effects of activity on a system of active patchy particles forming linear chains.⁴⁰ In Section 2, we report the numerical details of the system under study: a two-dimensional suspension of active Brownian repulsive particles whose surface is decorated with two diametrically opposed attractive sites that interact via a short-range attractive potential. In Section 3, we report the results, focusing on the structural features and on two observed novel states: active spirals and spinning crystalline clusters.⁴¹

2 Simulation details and analysis tools

We simulate a two-dimensional system of active bifunctional Brownian particles (ABBPs) in a square box with periodic boundary conditions. Particles are modeled as hard-disks with diameter σ, featuring two identical and diametrically opposed attractive sites, and self-propelling in the direction of the vector connecting the two sites (see Fig. 1a).

Fig. 1 (a) Pictorial representation of ABBPs: hard-disks featuring two identical and diametrically-opposed attractive sites (small golden disks) and self-propulsion (red arrow) along the segment connecting the two sites. (b) Interaction potential between two ABBPs in the most favorable bonding configuration, i.e., when two sites are facing each other (see inset). The dashed red line represents the hard-core plus square-well potential used as a reference for the choice of the interaction potential.

The two-body interaction potential between particles i and j is given by:

V(i,j) = V_CM(i,j) + V_S(i,j) (1)

where V_CM(i,j) represents the hard-core interaction between the centers of mass and V_S(i,j) represents the directional attractive interaction between the sites. Specifically, as in ref. 42, we choose:where r_ij is the distance between the centers of mass of the two particles and r^ab_ij is the distance between sites a and b located on particles i and j, respectively. The selected interaction potential incorporates the following assumptions: (1) particles are hard (m = 200); (2) the site–site V_S interaction resembles a square-well (n = 10); (3) the single bond per site condition is fulfilled (α = 0.12); and (4) the potential depth u₀ is set to 1 (ε = 1.001). The choice of α = 0.12 rises from purely geometric considerations. Indeed, geometric considerations for a three touching sphere configuration show that the choice of a well-width 0.119σ guarantees that each site is engaged at most in one bond.⁴⁰Fig. 1b depicts the shape of the interaction potential (black line) when two particles are in the most favorable bonding configuration.

Each particle is characterized by the position vector of its center of mass r = (x, y, 0) and the orientation angle θ representing the direction of the vector connecting the two sites with respect to the x-axis. The orientation vector η = (cos θ, sin θ, 0) is applied at each particle's center of mass and is restricted to rotate in the two-dimensional plane of the system.

While self-propelling in the same direction of the orientation vector with a constant speed v, each particle undergoes Brownian motion, in both position and orientation, at a constant temperature T. Thus, for a particle i, the translational and rotational equations of motion read as:

The diffusion coefficients D_T and D_R relate to each other via D_R = 3D_T/σ². In both the translational and rotational equations, the Gaussian white-noise terms are characterized by 〈ξ(t)〉 = 0 and 〈ξ(t)ξ(t′)〉 = δ(t − t′). The total force F_i and torque T_i acting on each particle are, respectively, given by:

where F_ij and T_ij are, respectively, the force and the torque between particles i and j interacting via the potential V_ij = V(i,j) described in eqn (1). The potential only depends on the distance between centers of mass r_ij and the orientations of both particles θ_i and θ_j.

In this article, all results are reported in reduced units. The unit length is σ (one particle's diameter, which is set to 1) and the energy unit is u₀ (the potential depth, which is also set to 1). With k_B = 1, temperature is measured in units of energy. Time is in units of σ²/D_T. All simulations are run for at least 10⁸ steps with an integration time step of 10⁻⁶ units.

We set the number of particles to N = 5000 and simulated the system at four different number densities ρ = N/A (with A the total area): ρ = 0.1, 0.2, 0.3, and 0.4. Activity is quantified by means of the Péclet number, defined as in ref. 43:

τ_R = 1/D_R being the reorientation time. We fix the value of the rotational diffusion (D_R = 3) and vary the value of the propulsion speed. Specifically, the Péclet number varies among the following values: Pe = 0, 1.66, 3.33, 5, 10, and 20 (corresponding to speeds v = 0, 1.66, 3.33, 5, 10, and 20, respectively). Simulations in the passive regime (Pe = 0) are performed as a reference. We choose to study the behaviour of the system at two temperatures: a lower one T = 0.07 and a higher one T = 0.1.

In order to study the assembly features of the suspension, we evaluate the chain length distribution ρ_ch and the cluster size distribution ρ_cl according to:

where N_l is the number of chains of length l, N_s is the number of clusters of size s, run, respectively, over all chain lengths and all cluster sizes and 〈⋯〉 averages over steady state configurations. The maximum values of l and s are fixed by the largest chain and largest cluster found in the entire simulation. On the one hand, the chain length distribution relies on an energetic criterion: two particles form a chain bond when their interaction energy is lower than −0.3 units. On the other hand, the cluster size distribution relies on a geometric criterion: two particles belong to the same cluster when the distance between their centers of mass is smaller than 1.2 units.

To describe the structural properties of the suspension, we compute the system structure factor:

where q is the exchanged wave vector, r_m is the coordinate of particle m and 〈⋯〉 averages over steady state configurations.

In order to be able to extract more detailed conclusions, we have decided to compare the S(q) calculated from the simulations (eqn (11)) with a theoretical S(q) representative of an ideal gas of polydisperse straight chains (eqn (12)). In the ideal gas limit, correlations between different chains can be neglected, and the structure factor of the system should be provided by the structure of a single chain, weighted by the appropriate chain length distribution:

where S_l(q) is the structure factor of a chain of length l:which, under the approximation that chains are straight, and averaging over all possible orientations of a chain, becomes:where represents the Bessel function of order zero. Any deviation we observe with respect to the theoretical predictions is due to a correlation between chains or the bending of the chains.

To analyze the bonding dynamics, we compare the orientation of the first and the second particle of each chain and assign +1 if the orientations are the same and −1 if they are not (i.e., we measure the scalar product between the two propulsion vectors and assign +1 if the product is greater than 0 or −1 if is less than 0). We repeat the same procedure for each pair of particles in the chain and then sum all values. Hence, for each chain i of length l, we obtain a value B_i ranging between +(l − 1) and −(l − 1). The upper limit +(l − 1) represents the case in which all particles are assembled with the same orientation and the lower limit −(l − 1) represents the case with alternating ones. We evaluate the average over all N_l chains of length l (E_B(l)) and the variance (Var_B(l)) as:

This method is intended to determine whether two particles prefer to assemble with propulsion vectors aligned in the same direction, in opposite directions, or randomly, rather than evaluating the chain propulsion. Given the flexibility of the chains, simply summing all propulsion vectors in one direction and subtracting the ones in the opposite direction would not adequately evaluate the chain propulsion.

When the formed structures are more compact, we evaluate the hexagonal order parameter,⁴⁴ whose expression, for each particle m, is given by:

where the sum runs over the k = 6 nearest neighbors and θ_mn is the angle formed by the vector r_mn and the x-axis. In particular, we are interested in: , where 〈⋯〉 average over steady-state configurations.

In the context of the crystalline structure, predominantly characterized by straight chains, we evaluate the chain propulsion as follows. For each chain, we compare the orientation of the first particle with the orientation of each other particle and assign +1 if they are the same (scalar product greater than 0) or −1 if they are not (scalar product less than 0). Then, chain propulsion is obtained by summing 1 (value for the first particle) to all values, once computed its absolute value. We normalize the chain propulsion by dividing it by the chain length, with 0 indicating no propulsion and 1 indicating maximum possible propulsion. Intermediate values provide insight into the chain's propulsion relative to its maximum possible value.

To demonstrate that a phase separation takes place, we assess the local density distribution by applying a Voronoi tessellation to the system.⁴⁵ Each cell in the Voronoi tessellation corresponds to the area of a particle identified by all points that are closer to that particle than to any other. The reciprocal areas of these Voronoi cells can be interpreted as local densities.

As far as we are aware, one can characterise a spiral-like structure quantifying the number of turns of the chain, by computing either the turning number³⁴ or the spiral number.⁴⁶ The two quantities, although defined in a slightly different way, identically give the same information.

The turning number is computed,³⁴ for each chain i with length l, as:

where β_j is defined by t̂_j = (cos β_j, sin β_j), which represents the bond unit vector t̂_j = (r_j+1 − r_j)/|r_j+1 − r_j|. Thus, (β_j+1 − β_j) gives the angle increment between two consecutive bonds. The turning number defines the transition from an elongated to a spiral state, by quantifying the number of turns of the chain between its two ends: χ_i = 0 (no turns), χ_i = ±1 (one turn), χ_i = ±2 (two turns), and so on. In particular, we are interested in the average turning number, defined as: , where N̄ is the total number of chains and 〈⋯〉 average over steady state configurations.

The spiral number,⁴⁶ for each chain i with length l, is defined as:

where α_l is the bond orientation of the last monomer and α₁ of the first one. To note that, in this case, the bond orientation α takes into account all full rotations. Therefore, α can be larger than 2π. The spiral number defines the transition from an elongated to a spiral state, by quantifying the number of turns of the chain between its two ends: s_i = 0 (no turns), s_i = ±1 (one turn), s_i = ±2 (two turns), and so on. Thus, the definition of the turning number is equivalent to that of the spiral number (see the ESI†).

To conclude, we underline that one could also use the end-to-end distance of the chain to characterise the spiral state. Our choice is to focus on the turning number in the main text as a way to characterise the system in a spiral state. The results on the spiral number are reported in the ESI.†

3 Results

We define the steady state as the state where the total number of bonds is stationary in time (the same applies for the total potential energy, see the ESI†). We cannot exclude that at longer time intervals a coarsening (or phase separation) could occur, but our evidence suggests that the system enters a stationary state where all static quantities (such as chain length and cluster size distributions) do not change over time. We start by analysing the phase behaviour of the suspension when varying activity and density. Fig. 2 reports two panels, each showing snapshots taken once the system was in steady state: the top one represents the system at a lower temperature (T = 0.07) and the bottom one at a higher temperature (T = 0.1).

Fig. 2 Steady state configurations, as a function of activity and density, at temperature T = 0.07 (top panel) and T = 0.1 (bottom panel). Activity increases horizontally (from left to right) and density increases vertically (from bottom to top).

Passive particles (left-most column) self-assemble into uniformly space-distributed linear chains (independently of temperature and density). With increasing activity, active particles self-assemble into linear chains that aggregate with each other. These dense and compact structures are more clearly visible at the highest density (top row in both panels of Fig. 2). At the lowest temperature, aggregates appear even at lower densities (bottom rows of the top panel of Fig. 2). For the largest simulated activity Pe = 20 (at the highest temperature, bottom panel), the system forms a compact and ordered structure.

Counterintuitively, this compact and ordered structure forms at the highest temperature: this is coherent with the fact that shorter chains are present at this temperature. This is further supported by the observation of this structure at the highest activity, when chains are the shortest.

We underline that the spinning crystalline phase differs from the well known MIPS phase for many reasons. First, this compact structure is rotating and translating, while the dense phase in MIPS does not rotate. Second, the spinning crystalline phase is composed of chains of different lengths instead of single particles as in the MIPS phase. Third, the spinning crystalline phase is characterized by a monocrystalline structure instead of a polycrystalline one, as in the MIPS case. To better characterize the region of the state diagram where the crystal is observed, we report a zoom of the state diagram in the ESI.†

At the lowest temperature and highest densities (rows ρ = 0.4 or ρ = 0.3 of top panel of Fig. 2), especially when the activity assumes low or mid-range values (such as Pe = 1.66, 3.33, or 5), chains aggregate forming spirals which are rotating at a finite angular velocity, reminiscent of the recently experimentally observed spirals in driven actin filaments on a motility assay.⁴⁷ We note that these spirals differ on the basis of their structures from the vortices detected in experiments of active filaments in ref. 48. These aggregates are very different from the density fluctuations observed in the suspensions of purely repulsive active Brownian particles (MIPS).²⁴ For this reason, we define this novel state with the acronym MISP, i.e., Motility-Induced SPirals.

The most common aggregated states are temporary chains of different lengths or compact spinning crystals. In the former case, i.e. when particles form a chain, they are able to adjust their propulsion direction while maintaining the bond with neighboring particles. However, if a particle's propulsion direction changes such that the particle can no longer stay connected to its neighbor, the chain breaks. This constrains how much particles in a chain can change their propulsion direction. On the other hand, a crystal is a more compact structure, even though it consists of several chains merged together. In this case, if particles change their propulsion direction enough to break the bonds within the chain they belong to, they are trapped by their neighbors and the compact structure does not break (unless they are located at the outer surface of the crystal and are free to move away).

3.1 Chain and active spiral phase

In either passive or active systems, particles self-assemble into linear chains whose length distribution decays exponentially (see Fig. 3a and c). In active systems, chains tend to be shorter than in the corresponding passive system. In particular, the higher the activity, the shorter the chains. This is observed at all temperatures and densities within the studied ranges (see the ESI† for the chain length distributions and cluster size distributions at different temperatures and densities). Therefore, activity clearly affects chain formation.

Fig. 3 Chain length distributions (left) and cluster size distributions (right) at temperature T = 0.1, density ρ = 0.3 (top) and density ρ = 0.4 (bottom), and all Péclet numbers studied (as indicated in the legend). Note that the bottom panels do not include the case Pe = 20 (red line), where the system is in a crystalline phase.

Moreover, as activity increases, a consistently more pronounced peak is observed at small chain lengths, likely due to the higher diffusivity of small chains as compared to long ones. This faster diffusion of smaller chains likely results in their quicker assembly with other particles, leading to a smaller number of small chains than what would be expected from the exponential trend.

On top of that, as density increases, chains aggregate forming clusters. At the highest density (see Fig. 3d), both passive and active clusters exhibit percolation. In fact, all cluster size distributions follow the same power law: ρ_ch(s) ∼ s^−τ with τ = 2.05, which is consistent within the random percolation universality class.⁴⁹

At the highest density (Fig. 3d), the cluster size distributions for all active systems present a peak for large cluster sizes (of the order of the total number of particles). A similar behavior has already been observed by the authors of ref. 50, which have demonstrated that the cluster size distribution of a system of self-propelled soft disks exhibits a peak when the system phase separates. The peak corresponds to a cluster size equal to the average number of particles in the dense phase. In contrast, the dilute phase contributes to the same cut-off power law observed in the homogeneous state.

Using a kinetic model in a finite-size system of active particles, the authors of ref. 51 and 52 have quantitatively demonstrated that active systems can exhibit not only an individual phase (characterized by a cluster-size distribution dominated by an exponential form) but also a clustering phase, characterized by a non-monotonic cluster-size distribution. In the latter case, a peak appears towards the tail of the distribution, which is an indication of particles aggregating in one large cluster. We observe the same features in our dense active system.

Instead, at a lower density (see Fig. 3b), only the clusters in the more active systems percolate. Indeed, the passive cluster size distribution follows an exponential power law, indicating that clusters simply coincide with chains.

As clearly reported in Fig. 3, neither the chain length nor the cluster size distribution show any relevant density-dependent behavior. For this reason, from now onward, we will mostly present our results at density ρ = 0.4 (indicating when not otherwise).

Even though clusters percolate in both passive and active systems, their structures differ significantly due to activity. Fig. 4 illustrates the structure factor S(q) at different Péclet numbers. As expected, due to the excluded volume effects, S(q) oscillates with a periodicity set determined by the diameter (first neighbors peak in Fig. 4).

Fig. 4 Structure factors S(q) at temperature T = 0.1, density ρ = 0.4, and all Péclet numbers studied, except Pe = 20, where the system is in a crystalline phase (see legend). The peak at qσ ∼ 3, indicating chain alignment, disappears as activity increases. This observation holds also at the other studied temperature and densities.

An important finding in S(q) of passive systems is the presence of a peak at qσ ≈ 3, which indicates the alignment within the chains. As activity increases, the peak shifts towards the first neighbors and, so, disappears. Its disappearance implies that chains are aggregating among each other instead of being uniformly distributed, as a characteristic distance between chains is no longer evident.

The peak at qσ ≈ 3 also vanishes when decreasing the density, implying that the chain alignment cannot generate a sufficiently strong signal when the system is too diluted. Fig. 5 shows the comparison of S(q) of passive systems at different densities with S(q) of an ideal gas of polydisperse straight chains (dotted black line). In the latter case, where chains are straight and non-interacting, the peak at qσ ≈ 3 is not expected.

Fig. 5 Structure factors S(q) for the passive system at temperature T = 0.1 and varying density as indicated in the legend. The dotted black line represents the structure factor of an ideal gas of polydisperse straight chains.

Chaining manifests in the non-negligible values of S(q) at small q. In particular, higher values are observed in the active case compared to the passive one (see Fig. 4). In the passive case, higher values are observed at lower densities (see Fig. 5), with particularly elevated values in the case of the ideal gas of polydisperse straight chains.

Non-negligible values of S(q) at small q are present. In particular, higher values are observed in the active cases (see Fig. 4). In the passive case, higher values are observed at lower densities (see Fig. 5), with particularly elevated values in the case of the ideal gas of polydisperse straight chains.

When analyzing the structure factor at different Péclet numbers, we do not observe any significant signal indicating the presence of spiral structures (see the ESI† for the analysis of the structure factor of spiral configurations).

To better understand chaining, we study the dynamics of the bonding process, to unravel whether particles have a tendency to self-assemble into chains with similar or opposite orientations and whether such tendency is related to activity. We will present our results at density ρ = 0.3 at which the system (for the chosen temperature and activity range) is never in a crystalline state. This allows us to present the full range of simulated activities. At low activities, E_B(l) ∼ 0 and Var_B(l) ∼ l − 1 (see Fig. 6). This means that B_i follows a Bernoulli distribution with equal probability of success (particles placed with the same orientation of the preceding one) and failure (particle placed with the opposite orientation of the preceding one). Hence, attraction dominates over activity in the bonding process, leading self-assembly to occur randomly regardless of the particle orientations.

Fig. 6 E _B(l) (left) and Var_B(l) (right) at temperature T = 0.1, density ρ = 0.3, and all Péclet numbers studied (see legend). As activity increases, E_B(l) and Var_B(l) take larger value consistently. This observation holds also for the other studied temperature and densities.

As activity increases, Fig. 6 shows that both E_B(l) and Var_B(l) consistently take larger values. This indicates, for every bonding event, an increase of the probability of two particles to self-assemble with the same orientation and a decrease of the probability to self-assemble with opposite ones. In this instance, self-assembly occurs favouring an alignment of the propulsion vectors. Hence, activity affects bonding, as we observe the propulsion vectors of two particles aligning when forming a bond.

At the lowest temperature and highest values of density and activity, the chains aggregate forming rotating spirals. A movie representing spiral formation is shown in the ESI.† The pathway for spiral formation is characterised by a few steps: (1) ABBPs self-assemble to form long chains; (2) chains aggregate due to their persistent velocity; (3) due to combined density and temperature effects, chains merge forming spirals, that spin due to the alignment within the chains.

We characterise the MISP state via the average turning number χ (being χ small when spirals are not present and large when spirals are present in the system). Setting the temperature at the lowest value, where we know the system can be in a spiral state, we report the average turning number χ as a function of density ρ for different Péclet numbers (Fig. 7a).

Fig. 7 (a) Average turning number χ as a function of density ρ at temperature T = 0.07 and all Péclet numbers studied at this temperature. (b) Typical spiral configuration (T = 0.07, Pe = 1.66, and ρ = 0.4) with particles color-coded according to the absolute value of the spiral number |χ_i|.

As expected, spiral formation occurs at higher density and smaller activity values. This is because when the system is too diluted (low density), the chains do not need to compete for space and move freely, and when activity is too high, chains are too short to coil into a spiral. In the ESI† we also show the probability density distribution of the average turning number.

To understand the reason why the values of the turning number reported in Fig. 7a are so small, we plot a typical snapshot of the system in a spiral state as shown in Fig. 7b: particles color-coded according to the absolute value of the spiral number. Spirals are usually formed by more than one chain and only a few of them are able to form strongly wound-up conformations. The reason for this is that these structures are quite unstable (unfolding and folding all the time) while also varying in size.

As in ref. 46, we have also computed the absolute value of the spiral number, averaging over all chains and over all configurations (Fig. 8 of the ESI†). Comparing the spiral number to the turning number computed for the same system, we observe they coincide, being both very small for all the chosen parameters (as shown in the snapshot of Fig. 7b).

Fig. 8 (a) Local density distributions at temperature T = 0.1 and density ρ = 0.4. Péclet numbers vary according to the legend. At the highest Péclet number (red color), the distribution is characterized by two distinguished peaks, which indicate phase separation. (b) Crystal steady-state configuration with particles color-coded according to the reciprocal volumes of the associated Voronoi cells.

3.2 Spinning crystalline cluster phase

At the highest temperature, density and activity (top-right configuration of bottom panel of Fig. 2), the system forms a crystalline cluster. Crystalline clustering is a two-step self-assembly process (first particles self-assemble into chains, next chains self-assemble into a cluster). For a more comprehensive overview of this process, in the ESI,† we show local density distribution computed for the location of the centers of mass of the chains in the crystalline configuration, together with a typical snapshot showing their location in the crystalline cluster.

Fig. 8a shows the local density distributions (evaluated performing a Voronoi tessellation of the system) at the highest temperature and density. A phase separation takes place at the highest activity (red line), as shown by the rise of two distinguished peaks. The bimodal behavior of the distribution is visible even though the distribution assumes non-zero values at mid-range densities, which is due to the contribution of particles located at the boundary of the crystalline cluster. Fig. 8b shows all particles of the crystalline configuration color-coded according to the reciprocal volumes of the associated Voronoi cells, clearly illustrating their influence on the local density distribution.

Interestingly, the pathway towards crystal formation follows the steps reported in Fig. 9 (also a movie of the entire process is shown in the ESI†).

Fig. 9a shows an initial state where chains start to form and aggregate but not in a stable way. Fig. 9b shows a stable cluster of chains with a head (bluish chains) and a tail (yellowish chains). Yellowish chains are chains where particles are pointing all in one direction and so have a non zero chain propulsion. Fig. 9c shows that the pushing chains in the tail allow the cluster to explore the system, leading to its growth due to the aggregation of other rather slow chains. All chains aggregate in a compact way and activity helps to anneal defects present in the cluster increasing its crystalline order (Fig. 9d).

Fig. 9 Snapshots taken along the crystallization process. Particles belonging to the same chain are depicted with the same color. A different color indicates a different value of the chain propulsion (as introduced in Section 2) divided by the chain length. (a) Chains are forming and aggregating but there is not a stable nucleus. (b) A stable nucleus is formed with a head (bluish chains) and a tail (yellowish chains). (c) The nucleus is moving and aggregating chains in the head. (d) The growing nucleus becomes a stable crystalline structure.

Once the crystal is formed, it is interesting to notice that it translates and rotates. To characterise its structure, we compute several properties. Fig. 10 depicts this steady state configuration in three different panels.

Fig. 10 Crystalline structure formed at the highest temperature (T = 0.1), density (ρ = 0.4), and Péclet number (Pe = 20). (a) Particles are colored according to the value of the crystalline order parameter ψ_i. While ψ_i = 1 represents perfect order, ψ_i = 0 no order at all. (b) Particles belonging to the same chain are depicted with the same color. A different color indicates a different value of the chain propulsion (as introduced in Section 2) divided by the chain length. While 1 indicates that all particles are arranged in the chain with same orientation, 0 with orientations alternated. (c) Particles are colored according to the value of the orientation vector along the x-axis. Note that cos θ_i = 1 indicates self-propulsion toward the right side of the box, cos θ_i = −1 toward the left side, and cos θ_i = 0 toward the top or bottom side. The big black arrows indicate the main direction of self-propulsion for each domain.

In Fig. 10a, particles are colored according to the value of the crystalline order parameter ψ_i (whose averaged value over all particles is ψ ≈ 0.87). In the core of the dense structure, particles are arranged as in a hexagonal lattice with disclinations. All particles not belonging to the crystal are monomers, dimers or trimers. In Fig. 10b, particles are colored with the same color when belonging to the same chain and according to the value of the chain propulsion (as introduced in Section 2) divided by the chain length. In the core of the dense structure, particles are arranged in straight chains, with alternated orientations in the innermost region and similar orientation in the outermost regions. This is due to the fact that, once the crystal core is formed, particles aggregate to it at the interface (see crystallization process in Fig. 9). In Fig. 10c, particles are colored according to the value of the orientation vector η along the x-axis. We note the presence of domains where particles are self-propelling in the same direction. The direction of such domains are indicated with black arrows. This results in an applied torque to the crystalline cluster. Thus, the crystalline cluster has a finite angular velocity other than a translating motion of its center of mass dictated by the evaporating front.

4 Conclusions

We investigate the phase behaviour of a model system made of active Brownian particles with two opposite-located short-range attractive sites. Our work explores the role of activity, temperature and density in the process of polymerization of active patchy particles in linear chains.

If particles are active, they self-assemble into chains which then aggregate (as opposed to uniformly space-distributed chains observed in the passive corresponding systems), forming from motility-induced spirals (lowest temperature and higher densities) to crystalline clusters (highest temperature, density and activity).

To characterise the structural features of the aggregated chains, we evaluate the cluster size distributions on an energetic and geometric basis. In particular, the first method (energetic bonds) allows us to characterize the length of the chains in function of the density, temperature and, most importantly, activity. The presence of activity reduces the average chain length at every temperature and density combination. On the other hand, the second method (geometric bonds) is evaluated with the purpose of characterizing the spatial chain aggregates. As a result, we observe the onset of a percolation phenomenon in a wide range of densities.

Then, we keep trying characterization the chain aggregates. This time we exploit a well-known quantity in the description of a system's structural properties, which is the structure factor. In all passive systems at the highest density, we observe the presence of an anomalous peak in the function that we attributed to the alignment of the chains. Furthermore, as activity increases, we observe that such peak shifts towards the first neighbors peak. This can be explained by the fact that when activity is introduced in the system, chains do not uniformly distribute but aggregate, and thus we cannot identify a characteristic length anymore.

Next, the analysis proceeds by investigating the arrangement of the particles within the chains based on the direction of the propulsion vectors. Specifically, our interest focuses on understanding whether the particles were bonding with propulsion vectors in the same or opposite direction. In systems with low activity, we find the probability of bonding in the same direction to be equal to the probability of bonding in the opposite direction, i.e., the attraction being dominant over activity in the bonding process. Interestingly, as activity increases, we note the probability of bonding in the same direction increases. This result is in agreement with our predictions of a bonding process being mainly determined by activity.To summarize, in passive systems, clusters are made of long, aligned and slow chains, while in active systems, clusters are made of short, aggregating and fast chains.

Finally, we focus on the formation and features of the crystalline structure observed at the highest values of temperature, density and Péclet number. In particular, we discover a significant rotation of the crystal cluster. The rotation arises from an effective torque generated by the presence of domains where particles are self-propelling in the same direction. What controls the nature of the fluid-to-solid transition in this active system surely deserves further investigation.

In this study, we can tune the flexibility of the chains by changing the angular aperture of the patch interaction. But in order to compare these results with those reported in our manuscript we would need to keep the bonding volume constant (since the dimerization constant depends on it). Unfortunately this cannot be done with the current functional form of the potential (which is square-well patchy) and would require us to switch to a different model (Kern–Frenkel⁵³). This will be the subject of a future research.

In conclusion, the presented results demonstrate the rich dynamics and emergent phenomena in active bifunctional Brownian particles highlighting the potential for a deeper understanding of out-of-equilibrium systems, and for novel applications in colloidal science.

Data availability

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

C. V. acknowledges fundings IHRC22/00002 and PID2022-140407NB-C21 from MINECO. J. R. and F. S. acknowledge support by ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by European Union – NextGenerationEU. The authors thank José Martín-Roca, Horacio Serna, Giulia Janzen and Daniel A. Matoz-Fernandez for insightful suggestions.

References

1. J.-P. Hansen and I. R. McDonald, Theory of simple liquids: with applications to soft matter, Academic press, 2013.
2. E. Bianchi, R. Blaak and C. N. Likos, Phys. Chem. Chem. Phys., 2011, 13, 6397.
3. L. Rovigatti, J. Russo and F. Romano, Eur. Phys. J. E, 2018, 41(59) DOI:10.1140/epje/i2018-11667-x.
4. F. Sciortino, Eur. Phys. J. B, 2008, 64, 505–509.
5. J. Russo, F. Leoni, F. Martelli and F. Sciortino, Rep. Prog. Phys., 2022, 85, 016601.
6. I. Nezbeda, J. Kolafa and Y. Kalyuzhnyi, Mol. Phys., 1989, 68, 143–160.
7. M. Ford, S. Auerbach and P. Monson, J. Chem. Phys., 2004, 121, 8415.
8. D. J. Kraft, R. Ni, F. Smallenburg, M. Hermes, K. Yoon, D. A. Weitz, A. V. Blaaderen, J. Groenewold, M. Dijkstra and W. K. Kegel, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 10787–10792.
9. R. P. Sear, J. Chem. Phys., 1999, 111, 4800–4806.
10. H. Liu, S. K. Kumar and F. Sciortino, J. Chem. Phys., 2007, 127, 084902.
11. E. G. Noya, I. Kolovos, G. Doppelbauer, G. Kahl and E. Bianchi, Soft Matter, 2014, 10, 8464–8474.
12. D. Tracey, E. Noya and J. Doye, J. Chem. Phys., 2021, 154, 194505.
13. E. Noya, C. Wong, P. Llombart and J. Doye, Nature, 2021, 596, 367–371.
14. H. Liu, M. Matthies, J. Russo, L. Rovigatti, R. P. Narayanan, T. Diep, D. McKeen, O. Gang, N. Stephanopoulos, F. Sciortino, H. Yan, F. Romano and P. Šulc, Science, 2024, 384, 776–781.
15. Y. Mai and A. Eisenberg, Chem. Soc. Rev., 2012, 41, 5969–5985.
16. Y. Zhang, J. Guo, C. Zhang and K. Yang, Adv. Mater., 2018, 30, 18721–18728.
17. S. Casalini, C. A. Bortolotti, F. Leonardi and F. Biscarini, Chem. Soc. Rev., 2017, 46, 40–71.
18. A. K. Boal, F. Ilhan, J. E. DeRouchey, T. Thurn-Albrecht, T. P. Russell and V. M. Rotello, Nature, 2000, 404, 746.
19. R. Mezzenga, Curr. Opin. Food Sci., 2015, 4, 294–308.
20. C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 045006.
21. W. Francis and D. Jörn, Nature, 2016, 536, 81–85.
22. W. Francis and D. Jörn, Nat. Commun., 2017, 8, 15169.
23. B. M. Mognetti, A. Šaric, S. Angioletti-Uberti, A. Cacciuto, C. Valeriani and D. Frenkel, Phys. Rev. Lett., 2013, 111, 245702.
24. J. Stenhammar, R. Wittkowski, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2015, 114, 018301.
25. F. Alarcon, E. Navarro-Argemí, C. Valeriani and I. Pagonabarraga, Phys. Rev. E, 2019, 99, 062602.
26. S. A. Mallory and A. Cacciuto, J. Am. Chem. Soc., 2019, 141, 2500–2507.
27. S. A. Mallory, F. Alarcon, A. Cacciuto and C. Valeriani, New J. Phys., 2017, 19, 125014.
28. M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys., 2015, 6, 219–244.
29. G. S. Redner, A. Baskaran and M. F. Hagan, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 012305.
30. A. Suma, G. Gonnella, D. Marenduzzo and E. Orlandini, Europhys. Lett., 2014, 108, 56004.
31. J. Schwarz-Linek, C. Valeriani, A. Cacciuto, M. E. Cates, D. Marenduzzo, A. N. Morozov and W. C. K. Poon, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 4052–4057.
32. A. Sokolov, M. M. Apodaca, B. A. Grzybowski and I. S. Aranson, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 969–974.
33. R. Di Leonardo, D. Dell'Arciprete, L. Angelani and V. Iebba, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 9541–9545.
34. G. Janzen and D. A. Matoz-Fernandez, Soft Matter, 2024, 20, 6618–6626.
35. G. Zhu, L. Gao, Y. Sun, W. Wei and L.-T. Yan, Rep. Prog. Phys., 2024, 87, 054601.
36. Z. Wang, Z. Wang, J. Li, C. Tian and Y. Wang, Nat. Commun., 2020, 11, 2670.
37. Z. Wang, Z. Wang, J. Li, S. T. H. Cheung, C. Tian, S.-H. Kim, G.-R. Yi, E. Ducrot and Y. Wang, J. Am. Chem. Soc., 2019, 141, 14853–14863.
38. Y. Wang, Y. Wang, D. R. Breed, V. N. Manoharan, L. Feng, A. D. Hollingsworth, M. Weck and D. J. Pine, Nature, 2012, 491, 51–55.
39. Z. Zhang and S. C. Glotzer, Nano Lett., 2004, 4, 1407–1413.
40. F. Sciortino, E. Bianchi, J. F. Douglas and P. Tartaglia, J. Chem. Phys., 2007, 126, 194903.
41. M. P. Holl, A. B. Steinberg and U. Thiele, arXiv, 2024, preprint, arXiv:2408.06114 DOI:10.48550/arXiv.2408.06114.
42. J. Russo, P. Tartaglia and F. Sciortino, J. Chem. Phys., 2009, 131, 014504.
43. J. Martin-Roca, R. Martinez, L. Alexander, A. Diez, D. Aarts, F. Alarcón, J. Ramirez and C. Valeriani, J. Chem. Phys., 2021, 154, 164901.
44. D. R. Nelson and B. I. Halperin, Phys. Rev. B: Condens. Matter Mater. Phys., 1979, 19, 2457–2484.
45. G. Voronoi, J. Reine Angew. Math., 1908, 134, 198–287.
46. R. E. Isele-Holder, J. Elgeti and G. Gompper, Soft Matter, 2015, 11, 7181–7190.
47. V. Schaller, C. Weber, C. Semmrich, E. Frey and A. R. Bausch, Nature, 2010, 467, 73–77.
48. A. Sciortino and A. R. Bausch, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2017047118.
49. D. Stauffer and A. Aharony, Introduction to Percolation Theory, Taylor & Francis, London, 2nd edn, 1992.
50. Y. Fily, S. Henkes and M. C. Marchetti, Soft Matter, 2014, 10, 2132–2140.
51. F. Peruani, A. Deutsch and M. Bär, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 030904.
52. F. Peruani and M. Bär, New J. Phys., 2013, 15, 065009.
53. N. Kern and D. Frenkel, J. Chem. Phys., 2003, 118, 9882–9889.

---

Footnote

† Electronic supplementary information (ESI) available: We report two movies to describe the novel active states. Movie 1 represents the spiral formation for the system at Pe = 1.66, ρ = 0.4 and T = 0.07. The spiral formation is characterized by a few steps: self-assembly of ABBPs in chains, chain aggregation and formation and breakage of the spirals. Movie 2 represents the formation of the spinning cluster for the system at Pe = 20, ρ = 0.4 and T = 0.1. The crystal formation is characterized by a few steps. After an initial chain formation and aggregation in an unstable way, it is possible to observe the formation of a stable cluster of chains characterized by a head (bluish chains) and a tail (yellowish chains). This comet-like cluster grows in size and order through the incorporation of new chains, while moving in the space available in the system. See DOI: https://doi.org/10.1039/d4sm00805g

---