Activation of colloidal patchy particle networks

Abstract

Active physical gels, exemplified by the cytoskeleton in muscle and plant tissues, are characterized by continuous energy injection, leading to rich but poorly understood non-equilibrium physics. Activated self-assembled colloidal architectures consisting of patchy particles and self-propelled particles can provide a well-controlled (experimental) model system that allows exploring the non-equilibrium behavior of such active physical gels. We conduct a numerical investigation of the effect of introducing self-propelled colloids modeled as active Brownian particles into a network-forming colloidal dispersion of dipatch and tripatch particles. We find a rich response of the self-assembled networks upon increasing activity. At low active forces, the networks form inhomogeneous void-rich structures. At medium active force, the network fragments into clusters of chains, and develops broad local density distributions. Finally, at high active force, the system exhibits motility-induced phase separation. These structural and dynamical responses are intimately related to the system's bond probability that can increase or decrease as a function of active force magnitude and direction, as well as attraction strength, affecting both the rate of bond formation and breakage. We discuss how our predictions compare to experiments.

---

1 Introduction

Active matter consists of individual components consuming energy to produce motion, resulting in collective behavior.¹ Examples range from biological phenomena such as the motility and morphogenesis of cells² and bird swarms, to synthetic systems such as colloidal Janus particles and colloidal rollers.^3–7 Active particles with different shapes,^8,9 repulsive potentials,^10,11 including attractive potentials,^5,12 alignment,¹³ or with different velocities,¹⁴ and mixtures of passive and active particles^9,15–18 have demonstrated a diverse phase behavior and dynamics. Being distinct from equilibrium systems, where dynamics arise from interparticle forces and thermal fluctuations, active matter is intrinsically out of equilibrium, leading to changes in the onset of phase transitions, and the emergence of entirely new dynamical phases such as flocking motion,¹³ non-equilibrium clustering,¹⁹ and motility-induced phase separation (MIPS) where active particles segregate into high and low-density regions.^20,21

Active physical gels, an intriguing class of active matter exemplified by the cytoskeleton found in muscle and plant tissues, are characterized by their viscoelastic properties and inherent activity, resulting in distinct non-equilibrium behavior.^22,23 Central to the cytoskeleton's function is the actin protein, which can polymerize, forming a variety of structures such as the crosslinked and branched networks seen in the cell cortex and lamellipodium.² The mechanical properties of actin, combined with its dynamic assembly and disassembly in conjunction with the molecular motor myosin, enable actin to play a central role in processes such as cell motility, replication, growth, and tissue repair.^24,25

Therefore, such soft biological materials are of great fundamental and technological relevance. Indeed, one of the promising directions in material science is to mimic driven biological systems in the form of active architectures, where active particles provide continuous energy injection. While much research is being done in this area, both on biological materials^25–30 as well as simplified physical systems,^31–36 there is a need for well-controlled model systems that would allow us to investigate the fundamental physical properties of such driven materials in a systematic way.

Experimental breakthroughs in nanostructure assembly and active matter provide such prototypical systems. Recent progress in colloid chemistry enable the design of micrometer-sized colloidal particles surface-decorated with patches of a different material.³⁷ Suspending such patchy particles in a near-critical binary mixture (e.g. water and lutidine) induces attractive directed bonds (i.e., only one bond per patch) between the patches on the surface of neighboring particles via a solvent-mediated critical Casimir force.^38,39 These attractive bonds allow controlled self-assembly into complex structure such as chains, rings, and networks.^40–44

While patchy particles experience thermal motion and adhere to Boltzmann statistics, they can be directly observed using techniques such as confocal microscopy.⁴³ Thus, they can be regarded as mesoscopic analogs of atoms⁴⁵ and can act as an experimental model system to explore complex self-assembled colloidal architectures analogous to their molecular counterparts.^40–44

Colloidal networks can be constructed by mixing patchy particles of different valencies, with an average valency beyond two, and can be viewed as a physical gel, characterized by reversible bond formation and adaptive character. Patchy colloidal gels create open architectures that achieve equilibrium,^46,47 in contrast to non-equilibrium gels made of isotropically interacting particles.^12,48–50 This distinction is for example evident in their spatially uniform dynamics.⁵¹ Patchy particle gels therefore are likely to show a different structural and dynamic response to active forces, compared to gels formed by isotropic particles.

Here, we aim to explore the behavior of physical gels under controlled activity, by combining such colloidal patchy particle architectures with well-controlled self-propelling particle systems, which allow experimental control of microscopic energy injection.^52,53 Examples of such particle systems are platinum-coated Janus particles that are self-propelled by catalysis of e.g. hydrogen peroxide, and colloids that are externally driven via external electric, magnetic or optical fields. These processes induce an active force aligned along the particle orientation. As the particle is still free to rotate in the suspension, such active particles’ dynamics are often modeled as active Brownian particles (ABPs).

Previously, we combined colloidal patchy particles with active Brownian particles to explore the collective non-equilibrium response of small colloidal architectures to the introduction of activity.⁵⁴ Introducing activity alters the dynamics, deforms the architectures, and enhances or reduces the breaking of colloidal bonds. Here we shift our focus to entire networks, and investigate how activity leads to the rearrangement of large-scale colloidal architectures. While we limit ourselves to modeling efforts, we stay close to the experimental conditions, to make a connection to the experimental realization of these systems.⁵⁵ Therefore, we employ an accurate potential model for patchy particle systems, capable of quantitatively predicting the outcome of (passive) experiments,³⁹ and combine this potential with the active Brownian particle (ABP) model, which mimics the behavior of the active particles. This means we have to employ Brownian dynamics instead of a Monte Carlo algorithms as used in ref. 39, to explore the phase behavior of activated colloidal patchy particle networks.

We find that the response to activity yields a rich behavior, in which activity can enhance as well as reduce the stability of the network architecture, induce inhomogeneity and alter the mechanism of fragmentation. The connectivity (i.e. bond probability) of the architecture determining the phase behavior displays non-monotonic trends, correlating with the attractive strength of patchy particles and the magnitude and direction of active forces. These findings corroborate well with our recent experimental work.⁵⁵

The paper is organized as follows. In the next section we introduce the colloidal model systems, the interaction potential, the simulation specifics and the analysis techniques. In the Results section, we present long simulations that explore the global behavior of activated colloidal patchy particle networks. We discuss the effects of activity on bond formation and dissociation rates, which together dictate the bond probability. We end with conclusions and a perspective on future research.

2 Methods

2.1 Model

2.1.1 Patchy particle models. In this work, we investigate self-assembled quasi-2D networks composed of colloidal dipatch (DP) and trigonal planar patchy (TPP) particles. As depicted in Fig. 1, DP particles have two attractive patches at opposite sites, while TPP particles possess three patches arranged in a planar configuration with a 120° angle between them. The DP particles are allowed to form bonds with both DP and TPP particles, but the TPP particles are restricted to only bind with DP particles. This is to prevent the TPP particles from clustering and forming small rings or honeycomb lattices.⁵⁶ In this way, the (active) TPP particles are distributed throughout the architectures and activate the passive DP chains. We consider two types of active TPPs: type I has the active force directed towards a patch, type II has the active force in the opposite direction (see Fig. 1). Type I and II could thus be classified as pushers and pullers, respectively.

Fig. 1 Illustrations of particle architectures. The divalent patchy (DP) particle exhibits a 180° angle between its patches, whereas the trigonal planar patchy (TPP) particle displays a 120° angle. There are two varieties of active TPP particles: one where the active force (depicted by the orange arrow) points toward one patch (type I) and the other where it lies between two patches (type II).

2.1.2 Interaction potential. The effective anisotropic pair interaction between two patchy particles i and j with orientation Ω_i and Ω_j, respectively, and interparticle distance r_ij is given by³⁹

V_pair(r_ij, Ω_i, Ω_j) = V_Yukawa(r_ij) + V_C(r_ij)S(Ω_i, Ω_j). (1)

Here, V_Yukawa signifies an isotropic repulsion, while the second term represents the directed patch–patch attraction between the particles i and j. It is important to note that we assume each particle pair can form only a singular bond. Given the relatively limited range and width of the patchy critical Casimir interaction, this condition is readily met in our systems, ensuring that only one combination of patches results in an effective attraction.

The anisotropy of the patch interactions is captured by two switching functions S′(θ). While these functions, in principle, depend on the orientations Ω of both particles, they are simplified to depend solely on the angles θ associated with each particle.

where k and l run over the n_p patches of each particle. The position of each patch in the particle reference frame is given by n_p unit patch vectors p, which point from the particle's center to the center of the patch as shown in Fig. 2. The angle θ is the angle between the patch vector and the interparticle vector as depicted in Fig. 2.

Fig. 2 A schematic illustration of the inter-particle vector r (dotted arrow), patch vectors p on each particle (solid arrows), and the angles θ.

Fig. 3 presents the optimized patchy particle potential in units of k_BT, where k_B is the Boltzmann constant and T the temperature. This potential is capable of reproducing the experimental system of dipatch particles.³⁹ The specific functional forms of the Yukawa electrostatic repulsion, critical Casimir attraction V_C, and the switching functions S′ can be found in the Appendix A and ref. 39.

Fig. 3 The patchy particle radial potential for dipatch particles composed of Yukawa repulsion V_Yukawa (eqn (14)) and critical Casimir attraction V_C (eqn (16)). The inset shows the switching functions S′ that are additionally a function of dT.

The total potential energy of the system is given by

where the first term sums over all patchy particle pair interactions between the N colloidal particles. The second term denotes the gravitational potential (eqn (18), Section A.2), mimicking the experimental system in which the particles sediment to the bottom surface, where z is the vertical distance to the surface. In addition, for the active TPP particles we include an alignment potential V_align(ê_A), with ê_A the unit vector of the direction of the active force (see Section 2.2.3). This potential mimics the experimentally observed fact that active Janus particles, when subjected to an AC field, exhibit a propelling motion and hover parallel over the capillary surface.³ To this effect we define the alignment potentialwhere ê_A,z signifies the z-component of the active force direction. This component is set to zero for all passive particles. The prefactor ε_align = 500k_BT was chosen such that a lift-off of the particles against the gravitational field does not occur.

The active force alignment of the tripatch particle (TPP) will cause one patch to always direct along the surface plane which is the patch in the direction ê_A or −ê_A, for type I and II, respectively. In quasi-2D systems, patches oriented along the surface plane show an enhanced binding propensity compared to patches facing the surface as shown previously in ref. 57. Therefore, in active systems, the reactivity of active TPP particles will be influenced by both their enhanced diffusion and the alignment. While in a passive system there would in principle be no active alignment, this would lead to a discontinuous jump from completely freely rotating TPPs, to completely aligned TPPs, the moment one turns on the activity, irrespective of its magnitude. To remedy this we have included the alignment force also in the passive case, so that in the limit of F_A → 0, the dynamics of the active systems approaches the dynamics of the passive system and any observed distinctions between passive and active systems are attributable exclusively to active propulsion effects.

2.2 Simulation details

2.2.1 System. Networks are constructed by mixing 70% DP and 30% TPP particles of in total N = 1000 or 500 patchy particles in a cubic simulation box of length 51.17 or 36.18σ with periodic boundary conditions, respectively. The strong attractive systems at dT = 0.12 K contain 1000 particles, while the others comprise 500 particles. The resulting quasi-2D system has an area coverage of roughly η = 30% or density ρ = 0.382σ⁻². This combination of density and ratio of DP : TPP particles was chosen as it resulted in a highly connected network at the strongest attraction strength at passive conditions (dT = 0.12 K).

2.2.2 Monte Carlo. We employ Monte Carlo (MC) to create and structurally equilibrate the colloidal network, with single particle and cluster moves. In a single particle move, a randomly selected particle is randomly rotated (50% of the single move MC steps) around their center of mass with a randomly selected rotation magnitude dq ∈ [0, θ_max] or randomly translated with a magnitude dr ∈ [0, r_max] (in the other 50% of the single moves). If the gravitational potential is applied, the z-direction of the random translation is divided by 10 to prevent placing the particle inside the surface wall or against the gravitational field which leads to high energies.

In addition to single particle moves (comprising 95% of all MC steps), cluster moves are performed, in which a cluster – a colloidal structure that is connected by bonds, i.e. E_pair < 0 for each bond – is rotated (50%) or translated (50%) in the x,y-plane only, due to the presence of gravity. In the first move, the cluster is rotated in the x,y-plane with an angle dq ∈ [0, θ_c,max] around a randomly selected particle in the cluster. In the second move, the cluster is translated in a random direction in the x,y-plane with magnitude dr ∈ [0, r_c,max]. If a bond is formed during a cluster move, the new configuration is rejected as an easy solution to obey detailed balance.⁵⁸ See Fig. 4 for an example of a self-assembled network equilibrated with MC.

Fig. 4 An example of a top view of a quasi-2D colloidal network composed of trigonal planar and dipatch particles in a ratio of 30 : 70 percent with color coding indicating its cluster size N_c, ranging from 1 (dark blue) to >4 (bright pink). The density ρ = 0.382σ⁻² and the patchy pair potential is strong (dT = 0.12 K).

2.2.3 Brownian dynamics. The dynamics of colloidal particles is described by overdamped Langevin dynamics, i.e. Brownian molecular dynamics (BMD), that includes the stochastic thermal motion of the particles.⁵⁹ We introduce activity by implementing the active Brownian particle (ABP) model that is widely used to simulate self-propelling colloidal particles. The translational equation of motion of a Brownian particle iswhere r⃑ is the positional vector of the particle, Δt is the timestep, t is the time, F⃑ is the conservative force acting on the particle coming from the potential V (eqn (3)), is the translational mobility tensor with D_T = 0.0034σ² per s as measured in an experiment for dipatch particles with diameter σ = 3.2 μm⁴² and is the identity matrix, β = 1/k_BT the inverse temperature, and [small xi, Greek, vector] is a Markovian vector where each element is i.i.d., with a zero mean and unit variance over time. The self-propulsion force F⃑_A acting on the center of mass of the colloidal particle is given by

F⃑_A = F_Aê_A, (6)

where F_A is the magnitude of the active force, and ê_A is a particle fixed unit vector that expresses the direction of the force.

Besides translational motion, the particle changes orientation. The rotational equation of motion is given by

where [small tau, Greek, vector] is the conservative torque acting on the particle derived from V (eqn (3)), and is the rotational mobility tensor with D_R = 0.05 rad² s⁻¹.⁶⁰ Note that for the active particles the alignment (eqn (4)) is included in the torque.

The activity of the ABP is expressed in this work by the active force magnitude F_A. Alternatively, one can compare the timescale related to the active velocity and the rotational and translational time using the particle Péclet number:

where v₀ (= F_AD_T) is the self-propulsion velocity, D_T is the translational diffusion coefficient and D_R is the rotational diffusion coefficient.⁶¹ The active forces 0, 10, 50, and 250k_BT/σ employed in this work, would correspond to Pe = 0, 2.6, 13.1, and 65.3. Note that the Pe number is also dependent on the effective radius of the particle when using a soft repulsion.³²

2.3 Analysis

We perform long straightforward BMD simulations on three independent samples and after reaching a steady state, i.e. when the bond probability does not drift over time, measurements are taken from the self-assembled architectures.

The static observables include: bond probability P_b = N_p,bound/N_p,tot, which is the number of bound sites N_p,bound divided by the total number of sites N_p,tot. The cluster size distribution, , where n_j is the number of clusters of size j, and is the total number of clusters. The bond occupancy of TPP particles, which is the distribution of the number of bonds the TPP particles make. Finally, the local density ρ_L is measured by drawing a square grid (with a grid cell size of approximately 26σ²) and counting the number of particles in each grid cell. (The grid size choice is a compromise between a too large grid cell, which would not lead to a local average, and a too small one, which would not lead to meaningful density estimates.)

The bonding of two patches can be seen as a (chemical) reaction in which the patch transitions from an unbound to a bound state:

where K is the bonding reaction constant, while ρ_b and ρ_u are the density of bound and unbound (free) patches, respectively. The k_ub and k_bu are the binding and unbinding rate constants, respectively. The bound state b is defined when a patch has an attractive interaction with another patch, i.e. their attractive part of the pair potential V_C·S < 0 (see eqn (1)). A breakage event is detected if two particles have separated more than half a diameter 0.5σ, which indicates the transition to the unbound state u. Note that k_ub and k_bu do not refer to a single binding or breakage mechanism but, instead, contain many possible breakage rates depending on the colloidal architecture and breakage location as observed in ref. 54. Note that this simple approach is thus different from analysis of a regular association process, which is more complicated.

In the brute force simulations, the state of each patch is monitored and the survival probability P_S, i.e. the probability that the patch did not change state, as function of time t is recorded. Assuming that binding and unbinding times are exponentially distributed processes and are independent events, dominated by a single timescale, the survival probability then follows the relationship

P_S(t) = exp(−t/τ), (10)

where the binding and unbinding rates, k = 1/τ, are subsequently deduced through a fit. Fig. 5 shows an example of such a measurement and fit of the bound and unbound survival probability.

Fig. 5 The survival probability of the bound and unbound state of the patches as measured in brute force simulations at dT = 0.12 K, F_A = 50k_BT/σ, and type II active particles. A fit of eqn (10) measures the binding and unbinding rate k_ub = 1/τ_ub and k_bu = 1/τ_bu, respectively. Solid colored lines are the measurements from three independent samples, and the black dotted lines are the fits.

We can use the resulting binding and unbinding rate constants, k_ub and k_bu, respectively, to determine the bond probability P_b in an alterative way as follows

Thus, we can explain the changes observed in the measured P_b in terms of the effects of activity on k_ub and k_bu.

3 Results and discussion

3.1 Dynamical phase diagram

3.1.1 Passive systems form (percolated) clusters of chains. We performed long BD simulations of a mixture consisting of DP and TPP particles in a ratio of 70 : 30, respectively, for different attractive Casimir interactions dT = 0.12, 0.16, 0.20 K, as well as for different active forces F_A = 0, 3, 10, 25, 50, 100, 250kT/σ. Initiated as a random fluid, the system forms a network structure after reaching equilibration with 1 × 10⁷ to 1 × 10⁶ MC cycles for the strongest to weakest interaction strengths, respectively. While we could also have used BD to equilibrate the system instead of using MC, the latter is more efficient due to the presence of global unphysical cluster moves. Snapshots of the resulting structures can be seen in the left most column of Fig. 6. Upon increasing the bond attraction from (dT = 0.20 K to dT = 0.12 K), a well-connected equilibrium colloidal network, or physical gel, forms, where chains of DP particles are interconnected via nodes composed of TPP particles (top left snapshot in Fig. 6 and 7).⁶² Despite achieving thermodynamic equilibrium, we still observe an overabundance of monomers with respect to the expected exponential chain length distribution. Previous work has indicated that this ‘anomalously’ high monomer concentration can be attributed to the relatively high rotational entropy of the monomers in the gravitational confinement inherent in the experimental setup we aim to model.⁵⁷ To follow the effect of activity on the phase behavior of the colloidal network, we conduct extensive Brownian dynamics simulations (lasting more than 20 000 seconds of simulation time) for various active force magnitudes, i.e., F_A = 0, 3, 10, 25, 50, 100, 250kT/σ, that reached a steady state in the bond probability after approximately 10 000 seconds of simulation time. We observe structural changes of the architectures attributable to active propulsion effects, as depicted by the snapshots at the end of the runs in Fig. 6 and 7 for type I and II active particles, respectively. We remind the reader that a type I tripatch particle has the active force aligned with a patch, while in type II the active force is anti-aligned, i.e. points away from a patch (see Fig. 1). We have color coded the snapshots background with similar structural and dynamical features (orange, blue and green).

Fig. 6 “Phase diagram” composed of snapshots of the colloidal architecture for type I active particles (with the active force directed to a patch, see Fig. 1) at various attraction strengths (dT) and active force magnitudes (F_A). The color of the particles indicates its cluster size N_c, ranging from 1 (dark blue) to >4 (bright pink). The shaded background color represents the architectures’ different global structure: homogeneous densities (bottom-left corner, orange), void forming, inhomogeneous densities (center, blue), and separated densities that form MIPS regions (top right, green).

Fig. 7 “Phase diagram” composed of snapshots of the colloidal architecture for type II active particles (with the active force directed opposite to a patch, see Fig. 1) at various attraction strengths (dT) and active force magnitudes (F_A). The color of the particles indicates its cluster size N_c, ranging from 1 (dark blue) to >4 (bright pink). The shaded background colors represent architectures’ global structure: homogeneous densities (bottom-left corner, orange), void forming, inhomogeneous densities (center, blue), and separated densities that form MIPS or nematic regions (top right, green).

3.1.2 Low to moderate active forces cause inhomogeneity. In the low to moderate activity regime and with both directions of active forces, an inhomogeneous structure with more dense and more dilute regions (voids) compared to the passive case appears, which is observed in the blue colored systems in Fig. 6 for type I and Fig. 7 for type II. What differentiates the inhomogeneous (blue) from the homogeneous (red) structures is the distribution of the local density ρ_L as depicted in Fig. 8. In the inhomogeneous structures, we observe values of local density, specifically at ρ_L < 0.15 and ρ_L > 0.6, that are not present in the passive systems. The onset of this inhomogeneous structure is almost immediate, i.e. for low F_A. For example, for strong bonds (dT = 0.12 K), void formation starts at F_A ≥ 3k_BT/σ (see Fig. 6 and 7 for snapshots of such systems). However, as bond strength weakens (dT > 0.12 K), the minimum force magnitude required for void formation increases. Such immediate structural responses to activity are not observed for isotropically interacting particles, as they quickly form a non-equilibrium gel at low temperatures, arresting the dynamics of the system.^12,49,50

Fig. 8 The probability distribution of the local density ρ_L at various active force magnitudes (F_A = 0, 10, 50, and 100k_BT/σ) for type I (a)–(c) and type II (d)–(f), and attraction strengths dT = 0.12 K (a) and (d), dT = 0.16 K (b) and (e), and dT = 0.20 K (c) and (f). The density distributions broaden upon increasing activity and even become bimodal for the largest force.

If the system is well-connected and active TPPs are anchored to the architectures, low activity leads to increased void formation w.r.t. moderate activity. Compare, for example, the snapshots at dT = 0.12 K and F_A = 10 and 50k_BT/σ in Fig. 6 and 7, and their local density distribution in Fig. 8a and d. Note that the distribution for F_A = 10k_BT/σ is unexpectedly wider than the other ones, reflecting the void formation in the network. While relatively mild activity initiates separation by persistently pushing the colloidal chains and clusters together, it also leads to enhanced bond breakage. Smaller (active) clusters possess faster effective rotations and translations compared to large clusters,^63,64 counteracting the void formation at higher active forces.

When the passive DP particles exhibit limited connectivity with the active TPP particles, especially in combination with type II active particles, the system resembles a system of passive particles immersed in an active bath.^49,50,65 Indeed, for dT = 0.20 K activity initially has hardly any effect on the local density distribution (see Fig. 8f).

3.1.3 High activity induces (nematic) phase separation. Increasing activity leads eventually to a completely separated phase at F_A = 250k_BT/σ. The separated phase is identified by the appearance of two peaks in ρ_L (light blue curves in Fig. 8). At dT = 0.12 K, due to high activity from F_A ≥ 100k_BT/σ onward, the colloidal network completely falls apart, as type II TPP particles break their bonds, thereby losing their ability to function as nodes connecting chains (Fig. 9b), while type I particles still exhibit significant binding (Fig. 9a). The passive DP particles and chains are then pushed together to high density regions; prompting them to aggregate into a nematic phase at F_A = 100k_BT/σ as seen in the green systems of Fig. 7, enlarged in Fig. 10b. Interestingly, the chains are significantly longer compared to lower activity. These chains themselves do not contain active particles, but the active TPPs push the passive dipatch particles towards each other so they can assemble and grow. In contrast, the corresponding type I system in Fig. 6 and 10a has many chains that do contain active tripatch particles, so that the chains themselves are pushed around as well, preventing formation of large clusters of DP chains. When increasing the activity further to F_A = 250k_BT/σ, the long nematic chains fragment again into shorter chains. Nevertheless, the high-density phase remains. The strength of the interaction between DP particles influences when this phase separation occurs; stronger interactions lead to earlier separation as the patchy interaction prevents the particles from moving back to less dense areas. Simultaneously, for high forces the active particles gather at the boundary of the cluster, due to the slow rotational diffusion. Such mechanism is then similar to the MIPS transition. One striking difference from standard MIPS is that active particles can also float atop of the cluster by overcoming the gravitational potential.

Fig. 9 Bonding occupancy of the active TPP particle type I (a) and II (b) at dT = 0.12 K.

Fig. 10 Enlarged snapshots from Fig. 6 and 7 at dT = 0.12 K and F_A = 100k_BT/σ, for type I (a) and II (b).

Despite the pronounced structural changes witnessed in the high activity regime, our primary interest lies in the intermediate, experimentally accessible,⁵⁵ region (up to F_A = 25k_BT/σ). In this regime, a significant portion of the network structure remains intact, albeit being strongly influenced by activity.

3.2 Kinetics of bond breaking and formation

3.2.1 Kinetics of an association reaction. The structure of the colloidal network, which determines its phase behavior, is directly related to the bond probability P_b. We observe a nontrivial behavior, either effectively reducing or increasing bond probability depending on bond strength, active force magnitude, and direction, as shown in Fig. 11a.

Fig. 11 The bond probability measured at steady state (a), the binding (b) and bond breaking (c) rate constants, and their ratio (d) at various active force magnitude and direction. Note that dT = 0.14 K is not shown for clarity in (b)–(d), as it is very similar to dT = 0.12 K. The inset of (c) shows a zoom-in of k_ub at dT = 0.20 K and an effective binding rate (black dotted line) calculated by multiplying the measured passive binding rate k⁰_ub with an effective diffusion and the normal diffusion D_T for the active TPP and DP, respectively.⁷ The shaded area represents a 95% confidence interval of the error of the mean from three independent samples.

To discuss this behavior we can relate the bond probability P_b to the binding and unbinding rate. At F_A = 0k_BT/σ, the measured binding rate constant k_ub is lower for systems with stronger attraction potentials due to fewer available patches, as seen in Fig. 11b. While appearing counterintuitive at first, this can be explained as follows. Consider the reaction 2P ⇌ P₂, where two free patches P form a bond P₂ with its equilibrium constant

defined by k_on and k_off, the association and dissociation rate constants, respectively.⁶⁶ Both reaction equations in eqn (9) and (13) describe the same process, and setting
yields k_ub = [P]k_on, and k_bu = k_off. Thus, k_ub as defined in eqn (9) is directly proportional to the density of free patches [P] in the system. This confirms the observed trend in Fig. 11b since the number of free patches decreases with their interaction strength.

3.2.2 Enhanced diffusion partly explains the effect on association. Upon increasing activity, the binding rate constant k_ub rises with the active force magnitude for all studied systems. Besides the free patch effect, this increased binding rate is partially attributed to the activity-enhanced diffusion of the TPP particles. To shed light on this underlying mechanism, consider the system at dT = 0.20 K, which contains predominantly free (TPP) particles. In this context, we can approximate the effect of enhanced diffusion on k_ub by assuming that rates are proportional to the diffusion constant k ∝ D as is the case for the Smoluchowski diffusion-limited reaction rate constant.^66,67 Given that the system has a 70 : 30 ratio of DP to TPP particles, the measured passive binding rate k⁰_ub is multiplied by the factor (0.3D^eff_T + 0.7D_T) with and v₀ = βF_AD_T.⁷ In the system with type I activity, the approximation is qualitatively similar to the measured k_ub at low activity (F_A < 5k_BT/σ) as shown in the inset of Fig. 11b by the black dotted line. Note that the effect of enhanced diffusion is not visible for the type II system, possibly because the active force direction is not aligned with a patch, leading to fewer binding events.

3.2.3 Breakage induces non-monotonic kinetics. The breakage rate constant k_bu at F_A = 0k_BT/σ is directly related to the potential depth, set by dT; deeper attraction translates into slower breakage, as depicted in Fig. 11c. However, upon increasing activity, we observe non-monotonic behavior: at weak attraction strength, type I active particles show reduced breakage rate constants k_bu. In contrast, type II particles and those with stronger attraction strengths display an increase in these rate constants.

To delve deeper into this topic, we recall that in ref. 54 we examined the influence of activity on bond breakage rates. The main results are summarized here. When the active force induces buckling, it becomes the prime factor in amplifying the rate of bond breakage. This particularly impacts the first stage of bond breakage, i.e. the particles escaping the potential well. The activation of the initial stage of bond breakage can lead to rate constants orders of magnitudes different compared to passive cases. This is, for example, evident in applying sliding forces on the decamer.⁵⁴ For systems that do not buckle—such as compressed dimers, extended forces applied to a decamer, or forces directed outwardly in a ring structure—the orientation of the active force plays a decisive role. Depending on its direction, either toward or away from the bonding volume, it can either suppress or enhance the bond breakage rate, respectively.

3.2.4 Cluster size distributions. In passive systems with strong attraction, the assembled clusters are relatively large with long DP chains, visible in the cluster size distribution in Fig. 12. (Note that these distributions also show the ‘anomalous’ overabundance of monomers,⁵⁷ although the effect is suppressed by the alignment of the tripatch particles.) Most TPP particles are bound to the colloidal architectures as shown in Fig. 9a. Upon increasing activity in both type I and II systems to F_A = 10k_BT/σ, forces are exerted on the colloidal architectures leading to the fragmentation of the clusters and a noticeable decline in triple-bound TPP particles. Both type I and II TPP particles mainly become double and single-bound (Fig. 9a and b). Akin to what is observed in the lifetime of (inward facing) colloidal rings in ref. 54, a maximum in double-bound type II TPP particles around F_A = 10k_BT/σ is observed (Fig. 9b). Moreover, the cluster size distribution at dT = 0.12 K (Fig. 12a) shows an increased numbers of small clusters upon applying activity in both force directions because DP–DP bonds break faster under the influence of activity as for example observed in decamers with sliding and extending forces in ref. 54. The cluster size distribution at dT = 0.12 K (Fig. 12a) indeed shows a shift to smaller clusters upon applying activity in both force directions.

Fig. 12 Cluster size distributions at dT = 0.12 (a), and 0.20 K (b). The magnitude of activity is F_A = 0, 10, or 50k_BT/σ and has directions as type I or II as indicated by the colors and symbols, respectively.

At activity levels of F_A = 50k_BT/σ, a clear distinction between type I and II particles is observed. While half of the type II TPP particles are not bound to any colloidal architecture (Fig. 9b), their potential to catalyze bond breakage via collisions remains. Conversely, type I TPP particles predominantly exist in a singly-bound state (Fig. 9a), bonding with DP particles to form what we term an “active chain”. By inspecting the snapshots, we identify two variants of these active chains: those with one or with two active particles capping a DP chain. Both chain varieties exhibit enhanced stability due to activity, evident from the increased prevalence of small clusters (of size N_c < 20) under both weak and strong attraction potentials as a function of F_A, as depicted in Fig. 12a and b.

Can we understand why these active chains are stabilized? For chains featuring a single type I particle, the active force, aligned along the patch bond, reduces the likelihood of particle separation i.e. the second stage of bond breakage (see ref. 54). Coupled with enhanced diffusion, the free patch at the other side of the chain becomes significantly more reactive yielding more short chains (N_c < 20). See Fig. 13 for an illustration.

Fig. 13 An illustration of how chains with type I particles experience increased reactivity due to enhanced diffusion, while type II, shows decreased reactivity due the diffusion in the opposite direction.

For chains with two active particles exerting a compressive force, our measurements in ref. 54 indicated that a decamer breaks faster at dT = 0.12, yet displays lifetimes similar to passive systems at dT = 0.22 K. For non-buckling chains, such as the dimers covered in ref. 54, or those made up of one or two DP particle paired with two TPP particles, compressing forces will extend their lifetime, again leading to enhanced stability for short chains at higher active force as observed in Fig. 12.

3.3 Summary of the simulation results

Network architectures, initially self-assembled by passive divalent colloidal patchy particles and trigonal planar patchy particles in a 70 : 30 ratio, were activated by modeling the trigonal planar patchy particles as active Brownian particles. Independent of the direction of the active force (type I or II), we observe that upon increasing the magnitude of the active force the network responds by forming inhomogeneous structures characterized by low density voids and high density clusters, broadening the local density distribution toward densities that are not observed in the passive systems. These inhomogeneities occur at active forces far below that required for MIPS. We find that the stronger the attraction strength of the patchy particles, the lower the active force magnitude needed to induce the inhomogeneties, which is correlated with the longer chain length in the strong attraction systems.

For the strongly attractive system, where the majority of active particles are bound within the colloidal architectures, an increase in the active force magnitude results in effectively fewer bonds (see Fig. 9a and 11a). Even though both the rate of bond breakage, and the rate of bond formation increase, their ratio confirms the trend in observed reduced bond probability. Increasing activity leads eventually to fragmentation of the colloidal network, and suppresses again the formation of voids and inhomogeneous structures. For the system at dT = 0.12 K, there is an optimal activity for F_A ≈ 10k_BT/σ at which a balance is struck: the architectures develop an inhomogeneous structure akin to phase separation but still have a sufficient number of bonds to retain a coherent network structure. As a result, this specific system displayed a broader local density distribution opposed to those at directly lower or higher activity levels.

In systems where bonding is sparse either due to weak patchy attractions or due to high activity for the strong attractive systems, there is a noticeable difference between the two active force directions. Type I particles, which have the active force pointing in the direction of one patch, tend to bind with DP particles, forming small, stable chains through enhanced effective diffusion of the chains and a suppression of the second stage of bond breakage as discussed in ref. 54. Interestingly, the formation of small stable chains with type I active particles leads to an effective decrease of the breakage rate for the weakly attractive systems at dT = 0.16 and 0.20 K compared to passive systems. In these active systems, we can observe a growth of bond probability as a function of active force magnitude before it declines again, something which is not observed for type II active particles (see Fig. 11a).

At high active forces F_A ≥ 100k_BT/σ, type II active particles become monomeric and push the passive DP particles into high density regions reminiscent of MIPS or, specifically at dT = 0.12 K, a nematic phase with long chains. The type I active particles, forming stable shorter chains, suppress the formation of this separated phase which is only observed at F_A = 250k_BT/σ and dT = 0.12 K.

3.4 Comparison with experiments

In ref. 55 we report on experimental observation of chain formation and fragmentation in activity-induced colloidal architectures. While these experiments employed a Janus particle instead of a TPP as the self-propelled particle, the findings are in line with our simulation results. The experiments show that upon increasing the active force until F_A = 12k_BT/σ the average chain length grows, caused by an increasing binding rate as well as a decreasing dissociation rate. This cause, in turn, is mostly explained by the enhanced diffusion mechanism in which active chains pick up new particles, and at the same time a force-induced decrease in breakage due to straightening up of the chains.

It turns out the experimental rate constants and bond probability roughly correspond to the type I system around dT = 0.18 K. Our simulations indicate chain growth for type I particles as a function of the active force magnitude (see Fig. 12b). Also, the relative binding rate constant k_ub increases for dT = 0.20, but barely increases for dT = 0.16 K (see Fig. 11b). The relative breakage rate constant k_bu decreases with active force magnitude, although for the dT = 0.16 K case we seem to observe non-monotonic behavior. The increased binding rate is attributed to the activity-enhanced diffusion of the TPP particles. Note, however, that we only can compare up to F_A = 12.

Further increment of the active force beyond F_A = 12k_BT/σ resulted in breakage of the network and chains. This qualitative behavior is also clearly observed in our simulations, e.g. qualitatively in Fig. 6 and also in Fig. 11a, in which P_b for the type I system at dT = 0.16 K first increases and then decreases again for F_A ≈ 20k_BT/σ. For the experimental temperature dT = 0.18 K this turnover will most likely shift to somewhat lower values, making it agree with experiments even better.

4 Conclusions

Drawing inspiration from biology, colloidal models of active physical gels, which instead of molecular motors are activated by self-propelling particles, enable the investigation of the response of physical gels to external stimuli. In this study, we conducted extensive numerical Brownian dynamics simulations on colloidal networks representing such active physical gels. Using accurate models for patchy particles, we predict that introducing activity into self-assembled colloidal network architectures leads to a rich response. Depending on the attractive patchy potential, and magnitude and direction of the active force, the system exhibits an increase or decrease in bond probability and displays homogeneous, inhomogeneous, and separated structures. At finite but low active forces the network deforms, but can still sustain its topology. In fact, the network is in some respect reinforced since the activity is inducing more bonds in some cases. In this stage the system starts to develop some inhomogeneities in density, as voids appear in the network, due to active particles pushing on chains and pushing them together. This behavior is especially prevalent in the strong interaction systems with long chains. This holds only up to a certain point – once the active force reaches a threshold, it begins to disrupt the network, leading to fragmentation, structural inhomogeneities, and eventually a motility-induced phase transition.

Recent experiments tested some of these predictions, and found remarkable similarity in the behavior, even though the details of the self-propelled particles were different.⁵⁵ The use of accurate effective patch potentials is certainly beneficial in this respect.

The observed increase in bonding probability and the onset of inhomogeneities at relatively low activity in this study, suggests that colloidal active gel networks have the opportunity to restructure themselves upon activation. Since the network topology is not destroyed, such behavior might be interpreted as a kind of resilience against activity, and a way to rebuild the network when activity is lowered again. Such repair dynamics in a far-from-equilibrium state could offer valuable insights into the adaptability of e.g. the actin network within the cytoskeleton, which can enhance our understanding of cellular resilience and functionality.

For future research one suggestion is to simply consider larger systems with more particles. The phase behavior may be shifted or influenced by the constrained size of the current system. A logical next step could be to change the density, or adjust the fraction of active particles. For a higher fraction of TPPs, we expect a denser network in which average chain length is shortened. For a significantly lower fraction of TPPs only disconnected clusters will occur. Of course, for very strong attraction, percolated networks could still form, in which nodes are connected by long chains of dipatch particles. Void formation and inhomogeneity due to activity will still occur for different TPP fractions, although the onset will likely require a larger active force for the open network at a lower fraction of TPPs. Additionally, we have not discussed any (global) relaxation pathways of bond breakage and formation. Future investigations might focus on such aspects.

As a final comment, our study did not include hydrodynamic interactions as these are only minor for the passive system due to the low particle Reynolds number. Moreover, we investigated the influence of activity using the ABP model, which in its standard form does not include hydrodynamic interactions. Future work might examine the effect of the hydrodynamics in the large activity regime.

Author contributions

HJJ – equal: conceptualization, formal analysis, investigation, methodology, writing – original draft, writing – review and editing. Lead: data curation, software, validation, visualization. NO – support: conceptualization, investigation, writing – review and editing. PS – support: conceptualization, investigation, supervision, writing – review and editing. Equal: funding acquisition, project administration. PGB – equal: conceptualization, formal analysis, funding acquisition, investigation, methodology, project administration, writing – original draft, writing – review and editing. Support: software. Lead: supervision.

Conflicts of interest

There are no conflicts to declare.

Data availability

The code for the Brownian dynamics simulations, and for the functions to measure the survival probability, local density distribution, TPP bond occupancy, bond probability, and cluster size distribution from trajectory-files can be found on Github https://github.com/hjjonas/pp_codes. Python codes for plotting the observations are available at Figshare with https://doi.org/10.6084/m9.figshare.29333234.v1.

Appendices

A Potentials

A.1 The patchy particle pair potential. The pair potential that governs patchy particles is composed of an isotropic repulsive component and a directed attractive critical Casimir component mediated through the patches (eqn (1)). In a previous study, dipatch particles were immersed at a surface coverage of 30% area in a 30–70% vol lutidine–water mixture with 1 mM MgSO₄ to form chains of specific lengths and stiffness depending on the strength of the interaction.⁴⁴ To reproduce the observed chain length distributions and bending rigidities across a range of weak to strong interactions as a function of temperature, we optimized the patchy particle potential based on the physical dimensions of the dipatch particles and theoretical critical Casimir potentials in ref. 39.

The repulsive potential is given by the electrostatic Yukawa potential:

withwhich is an exponentially decaying repulsion, when the particles do not overlap, i.e. r > σ. It is dependent on the charge of the particles Z = πσ_c²ϒ, the diameter of the charged colloids σ_c, the surface charge density ϒ, the Bjerrum length of the solvent λ_B = βe²/4πε with the permittivity of the solvent ε, and the elementary charge e. The screening length, i.e. Debye length, is where ρ_i is the number density of monovalent ions in the solvent.⁶⁸ The values of σ_c = 2.0 μm, λ_B = 2.14 nm, κ⁻¹ = 2.78 nm, and the optimized ϒ = −0.09e nm⁻² have been used in the model.³⁹

The attractive potential is based on the (isotropic) critical Casimir potential V_C between two spherical particles:

where A and B are fit parameters depending on the wetting set to w = 0.462 and temperature dT ∈ [0.12, 0.22] K and can be found in Appendix B1 in ref. 39. The wetting w is determined by the interplay between the patch material and the surrounding binary liquid, and the temperature dT is defined as its distance from the phase separation temperature T_cx of the binary liquid dT = T_cx − T.

The switching function was obtained by explicit integration using V_Yukawa and V_C at an effective patch width of θ^eff_p = 19.5° and fitted to the functional form:

where the coefficients c_l are given in Table VI in ref. 39.

A.2 External gravitational potential. The total potential is a summation over the pair potential, an alignment and an external gravitational potential (eqn (3)). The external gravitational potential is composed of two terms: a hard wall represented by a steep Lennard-Jones potential and a gravitational potential that depends, among other things, on the mass of the particle. See Appendix B2 of ref. 39 for more details. The resulting total gravitational potential is:where ε_LJ is an arbitrary (high) value set to 500k_BT and the gravitational force F_g is equal to −7.70k_BT/σ. Parameters b and z_cut are chosen such that both the potential energy and the force are continuous at z_cut.

Acknowledgements

The authors acknowledge the funding (Grant No. 680.91.124) from the Netherlands Organization for Scientific Research (NWO).

Notes and references

1. T. Vicsek and A. Zafeiris, Phys. Rep., 2012, 517, 71–140.
2. L. Blanchoin, R. Boujemaa-Paterski, C. Sykes and J. Plastino, Physiol. Rev., 2014, 94, 235–263.
3. D. Nishiguchi, J. Iwasawa, H.-R. Jiang and M. Sano, New J. Phys., 2018, 20, 015002.
4. J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine and P. M. Chaikin, Science, 2013, 339, 936–939.
5. I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger and T. Speck, Phys. Rev. Lett., 2013, 110, 238301.
6. A. Bricard, J.-B. Caussin, N. Desreumaux, O. Dauchot and D. Bartolo, Nature, 2013, 503, 95–98.
7. C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 045006.
8. S. E. Moran, I. R. Bruss, P. W. A. Schönhöfer and S. C. Glotzer, Soft Matter, 2022, 18, 1044–1053.
9. J. Eglinton, M. I. Smith and M. R. Swift, Phys. Rev. E, 2022, 105, 044609.
10. Y. Fily and M. C. Marchetti, Phys. Rev. Lett., 2012, 108, 235702.
11. G. S. Redner, M. F. Hagan and A. Baskaran, Phys. Rev. Lett., 2013, 110, 055701.
12. G. S. Redner, A. Baskaran and M. F. Hagan, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2013, 88, 012305.
13. T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen and O. Shochet, Phys. Rev. Lett., 1995, 75, 1226–1229.
14. T. Kolb and D. Klotsa, Soft Matter, 2020, 16, 1967–1978.
15. J. Stenhammar, R. Wittkowski, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2015, 114, 018301.
16. A. Wysocki, R. G. Winkler and G. Gompper, New J. Phys., 2016, 18, 123030.
17. R. Wittkowski, J. Stenhammar and M. E. Cates, New J. Phys., 2017, 19, 105003.
18. J. Stürmer, M. Seyrich and H. Stark, J. Chem. Phys., 2019, 150, 214901.
19. F. Peruani, A. Deutsch and M. Bär, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2006, 74, 030904.
20. M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys., 2015, 6, 219–244.
21. J. O'Byrne, A. Solon, J. Tailleur and Y. Zhao, Out-of-equilibrium Soft Matter, 2023, pp. 107–150.
22. J. Prost, F. Jülicher and J.-F. Joanny, Nat. Phys., 2015, 11, 111–117.
23. B. Najma, M. Varghese, L. Tsidilkovski, L. Lemma, A. Baskaran and G. Duclos, Nat. Commun., 2022, 13, 6465.
24. E. H. Barriga, K. Franze, G. Charras and R. Mayor, Nature, 2018, 554, 523–527.
25. F. Burla, Y. Mulla, B. E. Vos, A. Aufderhorst-Roberts and G. H. Koenderink, Nat. Rev. Phys., 2019, 1, 249–263.
26. M. P. Murrell and M. L. Gardel, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 20820–20825.
27. J. Alvarado, M. Sheinman, A. Sharma, F. C. Mackintosh and G. H. Koenderink, Nat. Phys., 2013, 9, 591–597.
28. C. Alkemade, H. Wierenga, V. A. Volkov, M. Preciado López, A. Akhmanova, P. R. ten Wolde, M. Dogterom and G. H. Koenderink, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, 1–12.
29. Y. Mulla, M. J. Avellaneda, A. Roland, L. Baldauf, W. Jung, T. Kim, S. J. Tans and G. H. Koenderink, Nat. Mater., 2022, 21, 1019–1023.
30. H. N. Verwei, G. Lee, G. Leech, I. I. Petitjean, G. H. Koenderink, R. M. Robertson-Anderson and R. J. McGorty, J. Visualized Exp., 2022, 1–27.
31. R. G. Winkler, J. Elgeti and G. Gompper, J. Phys. Soc. Jpn., 2017, 86, 101014.
32. D. Rogel Rodriguez, F. Alarcon, R. Martinez, J. Ramírez and C. Valeriani, Soft Matter, 2020, 16, 1162–1169.
33. S. Joo, X. Durang, O.-C. Lee and J.-H. Jeon, Soft Matter, 2020, 16, 9188–9201.
34. A. Militaru, M. Innerbichler, M. Frimmer, F. Tebbenjohanns, L. Novotny and C. Dellago, Nat. Commun., 2021, 12, 2446.
35. J. Zhang, T. Huang, G. Xu and Y. Chen, Commun. Theor. Phys., 2022, 74, 075601.
36. Y. Kim, S. Joo, W. K. Kim and J.-H. Jeon, Macromolecules, 2022, 55, 7136–7147.
37. Z. Gong, T. Hueckel, G.-R. Yi and S. Sacanna, Nature, 2017, 550, 234–238.
38. S. G. Stuij, M. Labbé-Laurent, T. E. Kodger, A. Maciołek and P. Schall, Soft Matter, 2017, 13, 5233–5249.
39. H. J. Jonas, S. G. Stuij, P. Schall and P. G. Bolhuis, J. Chem. Phys., 2021, 155, 034902.
40. 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.
41. T. Nguyen, A. Newton, D. Kraft, P. G. Bolhuis and P. Schall, Materials, 2017, 10, 1265.
42. S. G. Stuij, H. J. Jonas, Z. Gong, S. Sacanna, T. E. Kodger, P. G. Bolhuis and P. Schall, Soft Matter, 2021, 17, 8291–8299.
43. P. J. M. Swinkels, S. G. Stuij, Z. Gong, H. Jonas, N. Ruffino, B. van der Linden, P. G. Bolhuis, S. Sacanna, S. Woutersen and P. Schall, Nat. Commun., 2021, 12, 2810.
44. S. Stuij, J. Rouwhorst, H. J. Jonas, N. Ruffino, Z. Gong, S. Sacanna, P. G. Bolhuis and P. Schall, Phys. Rev. Lett., 2021, 127, 108001.
45. W. Poon, Science, 2004, 304, 830–831.
46. E. Bianchi, J. Largo, P. Tartaglia, E. Zaccarelli and F. Sciortino, Phys. Rev. Lett., 2006, 97, 168301.
47. J. Russo, F. Leoni, F. Martelli and F. Sciortino, Rep. Prog. Phys., 2022, 85, 016601.
48. J. Rouwhorst, P. Schall, C. Ness, T. Blijdenstein and A. Zaccone, Phys. Rev. E, 2020, 102, 022602.
49. S. Dikshit and S. Mishra, Eur. Phys. J. E: Soft Matter Biol. Phys., 2022, 45, 21.
50. B. Hrishikesh and E. Mani, Soft Matter, 2023, 19, 225–232.
51. E. Lattuada, D. Caprara, R. Piazza and F. Sciortino, Sci. Adv., 2021, 7, 1–7.
52. J. R. Howse, R. A. L. Jones, A. J. Ryan, T. Gough, R. Vafabakhsh and R. Golestanian, Phys. Rev. Lett., 2007, 99, 048102.
53. H. Seyforth, M. Gomez, W. B. Rogers, J. L. Ross and W. W. Ahmed, Phys. Rev. Res., 2022, 4, 023043.
54. H. J. Jonas, P. Schall and P. G. Bolhuis, Soft Matter, 2024, 20, 2162–2177.
55. N. Oikonomeas, H. J. Jonas, P. Schall and P. G. Bolhuis, 2025, submitted.
56. P. J. M. Swinkels, Z. Gong, S. Sacanna, E. G. Noya and P. Schall, Soft Matter, 2023, 19, 3414–3422.
57. H. J. Jonas, P. Schall and P. G. Bolhuis, J. Chem. Phys., 2022, 157, 094903.
58. D. Frenkel and B. Smit, Understanding molecular simulation, Academic Press, San Diego, CA, 3rd edn, 2023.
59. I. M. Ilie, W. K. den Otter and W. J. Briels, J. Chem. Phys., 2014, 141, 065101.
60. S. G. Stuij, PhD thesis, University of Amsterdam, 2020.
61. L. Fang, L. Li, J. Guo, Y. Liu and X. Huang, Phys. Lett. A, 2022, 427, 127934.
62. F. Sciortino and E. Zaccarelli, Curr. Opin. Colloid Interface Sci., 2017, 30, 90–96.
63. R. Soto and R. Golestanian, Phys. Rev. Lett., 2014, 112, 068301.
64. R. Soto and R. Golestanian, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 1–9.
65. L. Angelani, C. Maggi, M. L. Bernardini, A. Rizzo and R. Di Leonardo, Phys. Rev. Lett., 2011, 107, 138302.
66. A. Vijaykumar, P. G. Bolhuis and P. R. ten Wolde, Faraday Discuss., 2016, 195, 421–441.
67. M. Von Smoluchowski, Z. Phys. Chem., 1917, 92, 129.
68. V. A. Parsegian, van der Waals Forces, Cambridge University Press, 2005, pp. 1–380.

---