Activity affects the stability, deformation and breakage dynamics of colloidal architectures

Abstract

Living network architectures, such as the cytoskeleton, are characterized by continuous energy injection, leading to rich but poorly understood non-equilibrium physics. There is a need for a well-controlled (experimental) model system that allows basic insight into such non-equilibrium processes. Activated self-assembled colloidal architectures can fulfill this role, as colloidal patchy particles can self-assemble into colloidal architectures such as chains, rings and networks, while self-propelled colloidal particles can simultaneously inject energy into the architecture, alter the dynamical behavior of the system, and cause the self-assembled structures to deform and break. To gain insight, we conduct a numerical investigation into the effect of introducing self-propelled colloids modeled as active Brownian particles, into self-assembling colloidal dispersions of dipatch and tripatch particles. For the interaction potential, we use a previously designed model that accurately can reproduce experimental colloidal self-assembly via the critical Casimir force [Jonas et al., J. Chem. Phys., 2021, 135, 034902]. Here, we focus primarily on the breakage dynamics of three archetypal substructures, namely, dimers, chains, and rings. We find a rich response behavior to the introduction of self-propelled particles, in which the activity can enhance as well as reduce the stability of the architecture, deform the intact structures and alter the mechanisms of fragmentation. We rationalize these findings in terms of the rate and mechanisms of breakage as a function of the direction and magnitude of the active force by separating the bond breakage process into two stages: escaping the potential well and separation of the particles. The results set the stage for investigating more complex architectures.

---

1 Introduction

Structural architectures in living cells, such as the cytoskeleton in muscle or plant tissue, are both viscoelastic and active, i.e. undergo continuous injection of energy, leading to remarkable collective, non-equilibrium properties. Networks made of these living polymer filaments can be viewed as active gels,¹ a fascinating class of materials that show rich, responsive and functional mechanical properties, reminiscent of cell motility, replication and growth, and tissue repair.^2,3 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 (gel) architectures, where active particles provide continuous energy injection.⁴ While much research is being done in this area, both on biological materials^3,5–9 as well as simplified physical systems,^10–15 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. In this work we investigate activated breakage of self-assembled colloidal architectures relevant for networks and gels.

Experimental breakthroughs in nanostructure assembly and active matter provide such prototypical systems. For instance, novel synthesis routes enable the design of colloidal particles surface-decorated with DNA^16,17 or patches of a materials with different surface properties compared to bulk enabling formation of directed bonds.^18–21 Suspending the latter type of patchy particles in a near-critical binary liquid mixture (e.g. water and lutidine), induces attractive directed bonds between the patches on the surface of neighboring particles via a solvent mediated critical Casimir force. These attractive bonds allow the controlled self-assembly into complex structures such as chains, rings, and networks.^18,22–25 Experiencing thermal motion, such patchy particles obey the Boltzmann distribution, and thus can be seen as mesoscopic analogs of atoms. As they can be directly observed with e.g., confocal microscopy, patchy particles can act as an experimental model system to explore complex self-assembled colloidal architectures analogous to their molecular counterparts.^24–27

At the same time, well-controlled self-propelling particle systems allow experimental control of microscopic energy injection.^28,29 Examples of these particles are gold coated Janus particles that are self-propelled by catalysis of e.g. hydrogen peroxide, or colloids that are 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). Previous studies have addressed the self-assembly and motility induced phase separation of ABPs with isotropic interactions.^30–32

Here, we combine the above-mentioned two experimental breakthroughs in a simulation study to explore the collective non-equilibrium response of patchy colloidal architectures to the introduction of activity. We aim to obtain microscopic insight into the two main effects of introducing activity: (1) how does the active force lead to the breaking of colloidal chains? In particular, what effect does activity have on the bond configurations before escaping the potential well, and how does activity affect separating the particles into the bulk? Specifically, we look into the breakage mechanisms, – positions and – rates. (2) how is the dynamics of a colloidal architecture, such as a ring or a network, altered by such active forces? For that, we investigate global bending modes of ring structures.

To address these questions, we combine ABPs with patchy particles interacting via the critical Casimir force in Brownian dynamics simulations. To make a connection to the experimental realization of these systems, we would like to stay as close as possible to the experimental conditions. Therefore, we employ our accurate potential model for patchy particle systems that was benchmarked on experiments.³³

We find that the response to activity yields a rich behavior, in which activity can both enhance and reduce the stability of architecture, as well as alter the mechanism of fragmentation. In addition, although activity distorts the ring configuration, it leaves the dominance of global bending modes unaffected.

The paper is organized as follows. In the next section we introduce the systems, the potential model, simulation and analysis methods. In Section 3, we present and discuss the simulation results for dimers, decamers and rings. We end with concluding remarks.

2 Methods

2.1 Patchy particle architectures

In this work, we investigate two types of self-assembled structures: chains, and rings. The chains are composed of colloidal divalent patchy, or dipatch (DP), particles that have two attractive patches at opposite sites of the particle (Fig. 1). We consider two chain lengths: dimers consisting of just two, and decamers consisting of ten self-assembled DP particles. The outer two particles in the chains are made active by considering self-propelled forces acting on the centers of the particles in three qualitatively different ways, by compressing, extending, and sliding along the patch–patch bond (see Fig. 2). In contrast to a dimer, for which the activity is directly mediated by a single bond, the decamer can and must propagate the active force along the passive particles in the chain. For analysis purposes, we give each bond in the decamer an index; due to symmetry, there are only five irreducible bond numbers (Fig. 2).

Fig. 1 TPP particles have a 120 angle and DP particles a 180 angle between their patches.

Fig. 2 Schematic illustrations of the DP particle decamer (including irreducible bond numbers) and dimer. The orange colored particles are made active using three representative directions of active forces (orange arrows) opposing the patch (black arrow), along the patch, perpendicular to the patch indicated as extending, compressing or sliding forces, respectively.

The ring structures are symmetric and composed of four or six trigonal planar patchy (TPP) particles connected via 5 or 15 DP particles. Fig. 3 illustrates the structures T4D5 and T6D5 where the number after the T indicates the number of TPP particles, and after the D the number of DP particles connecting the TPP particles in the chain. Also the bonds in the rings are labeled; due to symmetry, there are only three or eight irreducible bond numbers for the 5 and 15 DP particle-based rings, respectively. The square structures with four TPP particles naturally contain an additional tension due to the mismatching bond angle of the TPP particle. Fig. 3 shows how, depending on the experimental realisation, the self-propulsion forces acting on the TPP particle either point toward the patch and creates outward-facing forces with respect to the ring, or between two patches resulting in inward-facing forces with respect to the ring. The ring simulations are performed in the presence of a gravitational potential that makes the particles sediment to a horizontal surface with a gravitational height of 0.13 times the diameter σ similar to the typical experimental setup employing critical Casimir forces.³³

Fig. 3 Schematic illustrations of the ring structures with four, illustrated with outward-facing forces, or six, illustrated with inward-facing forces, TPP particles connected by five DP particles named T4D5 and T6D5, respectively. The orange colored particles are made active TPP particles with their force directing towards a patch or between two patches creating outward-facing or inward-facing forces in the ring, respectively.

2.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)

where V_Yukawa denotes an isotropic repulsion and the second term denotes the relevant patch–patch attraction between the particles i and j, where we assume that each particle pair can only form a single bond. This condition is easily fulfilled for our systems, as the range and width of the patchy critical Casimir interaction is relatively small, and only one patch combination will result in an effective attraction.

The anisotropy of the patch interactions is captured by two switching functions S′(θ) that are in principle a function of the orientations Ω of both particles, but are simplified to a dependency on the angles θ of each particle

where 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. 4. Each particle pair can only form one bond, which is mimicked by the max function in S, resulting in the minimum attractive energy of the set of all possible patch–patch combinations in eqn (1).

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

The critical Casimir interaction is dependent on how close the system is to its critical point in terms of temperature T_c and composition c_c. In an experimental setup, typically, the composition is kept constant and the temperature is used to control the attraction between patchy particles. The temperature difference from T_c is captured by the parameter dT = T_cx − T, where T_cx is the coexistence temperature of the binary liquid (T_cx = 0.08 + T_c) and T is the temperature (T < T_c). The strength of the attraction increases upon lowering dT.³³

Fig. 5 shows the optimized patchy particle potential that is capable of reproducing the experimental system of DP particles. The exact functional forms of the Yukawa electrostatic repulsion, critical Casimir attraction, and the switching functions S′ can be found in the Appendix A and ref. 33 Appendix B.

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

The total potential energy of the system is the sum of all patchy particle pair interactions, and an external gravitational potential V_ext(z_i)

where i and j run over the N colloidal particles, and z is the vertical distance to the surface. Details about V_ext(z_i) can be found in ref. 33 Section II.D and Appendix B.

2.3 Equations of motion

The dynamical motion of colloidal particles is described by overdamped Langevin dynamics, i.e. brownian molecular dynamics (BMD), that includes the stochastic thermal motion of the particles.³⁴ The translational equation of motion of a Brownian particle iswhere r⃑ is the positional vector of the particle, and t the time. The timestep Δt = 5.0 μs was tested to be suitable for the strongest interaction strength at dT = 0.12 K, and used for all potentials. The conservative force F⃑ follows from the (gradient of the) potential V (eqn (3)), and the active force F⃑_A is only non-zero for the active particles. The translational mobility tensor with D_T = 0.0034σ² s⁻¹ as experimentally measured for dipatch particles with diameter σ = 3.2 μm, the identity matrix, β = 1/k_BT = 1 the inverse temperature with the Boltzmann constant k_B, and the temperature T.²³ The stochastic noise [small xi, Greek, vector] is a vector where each element is i.i.d., with a zero mean and unit variance over time.

The rotational equation of motion is given by

where [small tau, Greek, vector] is the torque acting on the particle coming from V (eqn (3)), and is the rotational mobility tensor with D_R = 0.05 rad² s⁻¹.³⁵ We integrate the equations of motion numerically.³⁴

A widely used model to simulate self-propelling colloidal particles is the active Brownian particle (ABP). It consists of a self-propulsion force acting on the center of mass of the colloidal particle

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.

In an experimental setup, the Janus particle aligns its motion parallel to the wall, due to the applied AC field. Furthermore, the continuous presence of gravity prevents the active particles to self-propel against the gravitational field away from the wall.³⁶ As the exact response of the Janus particles to the AC field is quite complex,³⁷ we mimic this behavior by including an effective alignment torque with magnitude^38–40

where ê_A,z denotes the z-component of the active force direction as defined above. Note that this torque acts as a restoring force that tries to minimize ê_A,z. The prefactor was chosen such that a lift-off of the particle against gravity did not occur.

The self-propulsion force F_A acting on the ABP is often expressed using the Péclet number:

with v₀ = βF_AD_T the self-propulsion velocity, D_T the translational diffusion coefficient and D_R the rotational diffusion coefficient.⁴¹ Note that the values for D_T and D_R are determined by the experimental system. In general, the ratio D_T/D_R follows the Stokes–Einstein relation. Deviations from this ratio are still possible,⁴² and would lead to different mechanisms and dynamics, but this is outside of the scope of this study.

The active forces F_A = 0, 10, 50, and 100k_BT/σ would correspond to Pe = 0, 2.6, 13.1, and 26.1, respectively, and are achievable in a typical experimental system. The Pe number is also dependent on the effective particle radius when using a soft repulsion instead of a hard-sphere model.¹¹

2.4 Monte Carlo

The BMD simulations are started from configurations that are equilibrated and decorrelated. Note that this passive initial distribution is often also taken as a starting point in experiments. Relaxation to a steady state is fast for low activity. Moreover, a non-equilibrium steady state approach prohibits the exploration of the high activity limit. To create this collection of starting configurations, we use Monte Carlo (MC).³³ For a given colloidal architecture, only single particle MC moves are performed and those MC steps that lead to bond breakage, i.e. E_pair ≥ 0, are rejected. In the single particle move, a randomly selected particle is randomly rotated (50% of the 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]. When the gravitational potential is applied, the z-direction of the random translation is divided by 10 to prevent placing the particle inside the wall or against the gravitational field which leads to unfavorable energies.

2.5 Analysis of the bond breakage process

Starting from equilibrated structures, we conducted straightforward, “brute force” BMD simulations to analyze bond breakage. Due to the random fluctuations or with the aid of the active force, a bonded particle pair can escape its attractive well. The particle pair either rebinds going back into the attractive well, or increases their interparticle distance to r = 1.5σ at which we assume the probability of rebinding is sufficiently low.⁴³ Only the latter criterion is considered a true bond breakage upon which the simulation is halted.

To analyze the bond breakage mechanism and lifetime, we separate the bond phase space into four regions indicating a strongly bound (1), weakly bound (2), diffusive (3), and truly broken (4) state as illustrated in Fig. 6. We define these four regions using the reduced energy E_λ, which is a parameter between zero and one:

E_λ = V_C(r_λ)·S/E_C,min, (9)

where V_C·S is the attractive (Casimir) part of the pair potential, and E_C,min is the value of the attractive part of E_min = V_pair(r_min,S = 1), which is the minimum of the pair potential (see eqn (1)), with r_min the corresponding interparticle distance. The distance parameter is set to r_λ = max(r,r_min), to prevent the λ-boundaries to lie in the repulsive part of the potential. See Table 1 for the criteria that define the four regions.

Fig. 6 Contour map of the bond breakage order parameter λ(r,S). An example trajectory of a breaking bond is shown with the white solid line for which its bond breakage mechanism is measured when crossing λ₂₃ as highlighted by the black circle. The four states strongly bound (1), weakly bound (2), free diffusion (3), and fully broken (4) as defined in Table 1. The breaking trajectory undergoes two stages: I. Escaping the potential well from state 1 → 3, and II. Separating the particles from state 2 → 4.

Table 1 The energy (E_λ in Eq. 9) and distance criteria defining the four region as depicted in Fig. 6. These values are chosen such that the flux from region 1 to 2 is reasonable, and region 2 is not yet in the free diffusion stage

	Region	Energy	Distance
1	Strongly bound	E λ ≤ 0.70
2	Weakly bound	0.70 < Eλ ≤ 0.01
3	Diffusive	E λ > 0.01
4	Truly broken		r λ > 1.50σ

---

We used the label λ as a (fictitious) order parameter that determines whether a configuration is inside a certain region, and, more importantly, when crossing a boundary between these regions. We denote the boundary (a.k.a. an interface) between region i and j as λ_ij. Initiated from the strongly bound region 1, trajectories can only escape via region 2 and 3 to finally reach region 4, crossing the interface λ₃₄, after which the bond is considered truly broken and the simulation is halted. This strict condition is necessary, as even after escaping the attractive well and entering region 3, the particles diffuse around and actually have a significant chance of rebinding. Counting trajectories that enter region 3 as being broken thus would severely underestimate the lifetime. However, since the time to diffuse the particles from the λ₂₃- to the λ₃₄-interface is very dependent on the location of λ₃₄, the instantaneous bond lifetime τ itself is measured as the last timestamp the trajectory crossed the λ₂₃-interface. We stress that the exact location of the interfaces does itself not influence the lifetime estimate.

The lifetime of the bond is strongly related to the rate constant for the breaking process. Considering bond-breaking as a two state dissociation process, we can express the corresponding rate as the inverse of the bonds' average lifetime (a.k.a. residence time) k_break = 1/〈τ〉. An alternative way to compute the rate constant for breakage is to compute the flux for trajectories starting in the strongly bound state to escape:⁴⁴

k_break = Φ_λ12P(λ₂₃|λ₁₂)P_sep, (10)

where Φ_λ12 is the (effective positive) flux through the λ₁₂-interface, P(λ₂₃|λ₁₂) is the conditional probability of escaping the potential well, i.e. reaching λ₂₃, given that the system comes from region 1, and P_sep is the probability the particles successfully separate, i.e. reaching λ₃₄. Respectively, they are defined as:where is the total duration of all the simulations used for the estimate, N_λ12 the number of positive crossings through the λ₁₂ interface, i.e. from region 1 toward 2, and N_1→3 (N_2→4) is the number of times region 3(4) is reached given that one came from 1(2), or, in other words, the number of first crossings of λ₂₃ (λ₃₄) subsequent to λ₁₂ (λ₂₃). One can quickly verify that the breakage rate is then, as expected, with τ_sim the average simulated bond lifetime.

Measuring the quantities in eqn (10)–(13) per bond provides microscopic insight into the influence of activation on the breakage dynamics. The first stage of bond breakage, i.e. escaping the potential well, is defined by Φ_esc = Φ_λ12P(λ₂₃|λ₁₂) which expresses the rate at which the particles escape their potential well. The second stage of bond breakage, i.e. successfully separating the particles to r = 1.5σ, is expressed as the probability P_sep. Separating the breakage process into two stages helps to identify the role of activity as stage I is mainly defined by the shape and depth of the pair potential and buckling of the bond, while stage II by, for example, the position of the bond in the structure, or crowding effects by other particles.

As patchy particles have rotational and translational degrees of freedom, one can characterize the mechanism of breaking a critical Casimir induced bond between two particles as a combination a two qualitatively different limiting cases. Fig. 7 schematically illustrates the two breakage mechanisms along the λ(r,S)-order parameter: (1) a pure radial mechanism, where particles move away from each other along the bond vector, keeps the particles' patches aligned and S = 1. And (2) a purely rotational, angular mechanism in which particles rotate along each other perpendicular to the bond vector, until the overlap area between the patches is vanishing, makes S = 0. In practice, both routes are available to the particles and the breakage mechanism is naturally a combination of these limiting routes. Accordingly, we can follow the breakage mechanism along the value of S as it can distinguish between radial (S = 1) or angular (S = 0) trajectories.

Fig. 7 Two distinct breakage mechanisms illustrated along the λ(r,S)-order parameter. Starting from the bound state (top left corner), the bond breaks via the radial mechanism keeping the patches aligned (S-value close to 1), or via the angular mechanism where the patches rotate away (S-value close to 0). Both routes may lead to the broken state or rebind back into the bound state.

For analyzing the breakage mechanism, we again use the separation of the two stages of breakage. In the first stage, the system escapes the potential well and reaches the λ₂₃-interface for the first time at a particular S-value, which we denote S_mech. These S_mech-values are counted as n_1→3(S_mech) which' probability distribution is given by eqn (14). Note that the counts N_i→j from eqn (12) and (13) equal the sum over all S_mech-values as . Continuing the simulation, the trajectory either rebinds, adding a count to n_3→1(S_mech), or breaks, adding a count to n_2→4(S_mech). Then, noting that n_1→3(S) = n_2→4(S) + n_3→1(S), the probability of breaking at given S is written in eqn (15). Only upon true breakage, the last n_2→4(S_mech) is counted and normalized as the breakage mechanism probability P_mech in eqn (16). One can see that the breakage mechanism probability also follows from P_mech(S) = P_1→3(S)P_2→4(S)P_sep⁻¹ where P_sep⁻¹ is the normalization.

By only considering the values of S at λ₂₃ under the condition the system came from region 1, λ₂₃-recrossing events are filtered out of the breakage mechanism analysis. Recrossing events, as observed for the example trajectory in Fig. 6 where the trajectory recrosses λ₂₃ at S ≈ 0.1 after escaping the well at S ≈ 0.3, are mainly determined by free diffusion and do not contain information on how the system escaped the potential well.

Many of the observables, e.g. the bond lifetime, show an exponential decay over time. Therefore, we report averages with a 95% confidence interval of the mean calculated using block averaging. As the start configurations are equilibrated before the measurements start, and the breakage itself is a rare event, all breakage events are uncorrelated and we take blocks of minimal 100 simulations.^45,46

3 Results and discussion

In this section, we present and discuss the numerical simulations. We start with the simplest colloidal architectures in Sections 3.1 and 3.2: the dimers and decamers. We performed 1000 brute force BMD simulations at various attraction depths, tuned by the critical Casimir interaction via the temperature dT. The activity is tuned by the force magnitude F_A and direction (see Methods and Fig. 2).

In Section 3.3, we increase the complexity of the architectures toward rings, and evaluate the influence on strained and relaxed rings of different sizes. Here we include the gravitational field to mimic the experimental setup.

3.1 Dimer

Before the bond breaks, activity naturally affects the exploration of phase space and we show this effect with probability histograms in Fig. 8. Measurements are taken at each timestep during BMD simulations in a volume or isosurface of interest. For the effect on the bond orientation of the dimer, we measure S inside the bonding volume, i.e. region 1 and 2, at which the bond is intact (Fig. 8a) and we observe that activity does not significantly buckle the bond as the histogram of S of passive and active dimers overlap.

Fig. 8 Exploring the potential well of passive and active dimers at dT = 0.12 K. The probability histogram of S (eqn (2)) of a passive dimer (red solid line) and active dimers with F_A = 100k_BT/σ and force directions: extending, compressing, and sliding. The measurements are taken in (a) the full bonding volume, i.e. region 1 + 2, and (b) when crossing (positively and negatively) the λ₂₃-interface. There are 50 bins for S ∈ [0.0,1.0]. As the λ₂₃-interface has a minimal value of S = 0.01 (at E_λ = 0.01), the first bin in (b) is approximately twice as low as the second bin.

The histogram shown in Fig. 8b shows the distribution of the S-values at crossing the λ₂₃-interface positively and negatively. In equilibrium, this histogram, if taken as − ln(P), represents the free energy in units of k_BT and we observe that bonds mostly sample values with θ > 0°. Each particle in the bond can rotate around the interparticle vector to maintain S′(θ) and r at a constant value (Fig. 4), resulting in a bonding volume proportional to sin θ. Such behavior leads to an entropic contribution to the free energy, and configurations with θ > 0° are found to be more prevalent.

Activity has a specifically large effect where the gradient of the potential is comparable or smaller than the active force, such as when crossing λ₂₃ (Fig. 8b). Under extending forces, the probability density shift to higher S-values compared to the passive case. This has two reasons: first, the particles start out in the strongly bound region 1 where S is close to one and activity pulls the particles apart, keeping S high. Second, recrossing λ₂₃ at low S-values will be suppressed for extending forces as they direct away from the λ₂₃-interface. Both effects are indirectly measured using P_1→3 and P_2→4, see below. Conversely, under compressing forces, the active forces' direction is reversed and thus more density is observed at low S-values.

Breakage mechanism. Fig. 9 presents the probability distributions P_1→3(S), P_2→4(S) and P_mech(S) (eqn (14)–(16)) for passive and active dimers. For the passive case, the attraction strength is varied between dT = 0.12–0.22 K. For the active case, the active force magnitude is set to F_A = 100k_BT/σ and the dimers experience extending, compressing and sliding force directions at dT = 0.12 and 0.22 K.

Fig. 9 Analysis of the breakage mechanism of the dimers by P_1→3, P_2→4, and P_mech from eqn (14)–(16), respectively. (a)–(c) Passive dimers at dT = 0.12–0.22 K. (d)–(i) Active dimers at dT = 0.12 (d)–(f) and 0.22 K (g)–(i) with extending, compressing, and sliding active forces at F_A = 100k_BT/σ. Each row shows its legend on the right. Colored bands indicate 95% confidence interval.

Starting with the passive case, Fig. 9a shows that to reach λ₂₃, the patches rotate significantly more in the strong attraction potential (dT = 0.12 K) compared to the weak attraction potential (dT = 0.22 K). At the same time, the conditional probability to break P_2→4(S) peaks around S ≈ 0 (i.e. the patches are rotated away), decreases quickly, and plateaus at a low value for S > 0.5 (Fig. 9b). The peak at low S is reduced for stronger attraction. Since P_1→3(S) and P_2→4(S) show opposing trends, the effect of dT on the breaking mechanism is reduced (Fig. 9c). Nevertheless, the majority of bond breakage events clearly occur via the angular mechanism. The maximum of P_mech(S) shifts from S_max ≈ 0.2 to S_max ≈ 0.1 with increasing attraction strength.

Although the translational and rotational diffusion constants have no influence on equilibrium distributions, they do play a critical role in bond breakage dynamics as demonstrated in ref. 42 and 47. Hence, the introduction of active particles that experience directional enhanced translation is expected to significantly impact the breakage mechanism. As hypothesized above, when using extending forces, P_1→3(S) shows a shift toward high S-values (Fig. 9d and g), indicating radial breakage, as the active force pulls the particles from high S in region 1 toward the λ₂₃-interface. Subsequently, the probability of successfully separating the particles increases for all values of S (Fig. 9e and h), as extending forces are directed away from the bonding volume (see Fig. 10 for a schematic illustration). The breakage mechanism probability distribution in Fig. 9f and i shows that for extending forces both angular and radial routes are likely, with only a small bias toward smaller S, i.e. angular breakage. In contrast, for compressing forces, P_1→3(S) shows a shift toward lower S-values, indicating angular breakage, with respect to the passive case. Moreover, the corresponding P_2→4(S) indicates that configurations with S > 0.5 are completely blocked from particle separation as the compressing active force restores the particles back into the bound state (Fig. 10). The breakage mechanism probability distribution P_mech(S) shows now a significant shift toward lower S-values indicating that all bonds must break via the angular route. Finally, sliding forces results in a distribution P_1→3 similar to that of compressing forces, while P_2→4 shows a large increase at S < 0.5. Therefore, one would expect a shift of P_mech(S) toward angular breakage. However, as sliding forces also increase P_sep = N_2→4/N_1→3 (see Fig. 11d). the breakage mechanism probability P_mech(S) is similar to the passive dimers, even if the two stages of breakage are very different from passive breakage.

Fig. 10 For a configuration at λ₂₃, its probability to rebind depends on the overlap between the bonding volume (purple region between patches), and the positional density distribution (orange halo around the particles) which' probability is amplified in the direction of the active force(orange arrow). This is a schematic representation and not based on simulations.⁴⁸ See Fig. 9 for measured P_2→4(S) = 1 − P_rebind(S) at the various force directions.

Fig. 11 Passive and active dimers. (a) The measured average bond lifetime τ_sim in the passive case at dT = 0.12–0.22 K and in the active cases at dT = 0.12 and 0.22. The scaled rate constant (b), stage I (Φ_esc) (c), and stage II (P_sep) (d) at dT = 0.12 and 0.22 K and F_A = 10, 50, 100, and 250k_BT/σ. The color coding indicates the dT, similar to Fig. 5, and the line coding is shown in figure on the right.

This analysis clearly indicates that the breakage mechanism of dimers is strongly affected by the separation probability P_2→4, both in the passive or active case, and at strong or weak interaction strength.

Bond lifetime. Fig. 11a–d presents, respectively, the measured average bond lifetimes in the simulation τ_sim, the ratio of the active and passive rate constants the flux out of the potential well (stage I) Φ_esc, and the separation probability P_sep of stage II. Naturally, in the passive case, the bond lifetimes τ_sim depend strongly on dT, as a shallower potential leads to exponentially faster breakage (Fig. 11a).^49,50 Under the influence of activity, τ_sim can either be enhanced or suppressed depending on the direction of the applied force.

As expected, only stage I of breakage (Fig. 11c), i.e. the rate of escaping the potential well, depends on the potential depth as determined by dT. Compressing and extending forces, respectively, exponentially decrease and increase the rate constant of stage I proportional to the magnitude of F_A. Because in the dimer the particles are bound with directed patches, forming a relatively stiff bond, and which maintain the active forces at a ±180 angle, the effect of activity inside the potential well of the dimer resembles an equilibrium problem in an effective tilted potential along the direction of r.⁵¹ This effect can be qualitatively captured with an Arrhenius type expression

k_tilt(F_A) = νe^{−βV(r)+βαF_AΔr} (17)

with α = −1, 1 for compressing and extending forces, respectively, and Δr the distance from the potential minimum to the barrier location. As this distance is roughly Δr ≈ 0.007σ, the expected increase/decrease is roughly a factor of 2 for F_A = 100, qualitatively in agreement with the observation (see Fig. 11c). Note that this analysis does not hold for the sliding force case, where the force direction is not aligned with the interparticle distance vector.

The effect of activity on stage II of breakage, i.e. successful particle separation, was qualitatively explained using Fig. 10 and quantitatively measured in Fig. 9. Fig. 11d confirms that extending and sliding forces enhance the probability of separation P_sep, while compression reduces it. Additionally, P_sep does not show a potential depth dependency.

Considering that activity can be viewed as an effective tilt on the potential, and P_sep is dT independent, the contribution of the potential depth can be separated from the activity by dividing the breakage rate constant k_break by the rate constant at zero activity, as shown in Fig. 11b as . Consequently, for the active dimer, activity with the same magnitude and direction have an equivalent effect on the breakage rate.

3.2 Decamer

In contrast to a dimer, a longer chain can transmit stresses resulting from active forces exerted on its outer particles through the other bonds in the chain (Fig. 2). Although the breakage mechanism of the bonds in decamers remains qualitatively similar to that of dimers (see Fig. 18 in Appendix B), activity has a significant impact on the dynamics and conformational distributions of the chain, leading to opposing effects in bond lifetime compared to a dimer.¹⁰

Before any bond has broken, activity already affects the bond distributions, causing the chain to buckle or straighten.⁵²Fig. 12 shows the probability distribution of S for each bond in a passive chain and active chains with F_A = 100k_BT/σ at dT = 0.12 K. In the freely moving passive chain, all nine bonds are statistically identical, in line with Boltzmann's equipartition theorem, and appear as a single (red) curve.

Fig. 12 Deformations in decamers at dT = 0.12 K. The probability histogram of S (eqn (2)) of a passive chain (red solid line) and active chain with F_A = 100k_BT/σ (dots) and the force directions: extending (a), compressing (b), and sliding (c). The measurement is taken at regions 1 and 2. The colored dots indicate the reduced bond number as shown in the legend.

Passive chain. Although the bonds are energetically identical, both breakage stages show an enhancement at the passive chain's end, leading to bond 1 breaking twice as likely compared to the bonds in the center, as the breakage rate constant of bond 1 is approximately twice as high as those of bonds 2–5 (see Fig. 13e). The breakage profile of passive semi-flexible chains has been shown to depend on the non-linearity (i.e. anharmonicity) of the interaction potential, stiffness of the chain, and ratio between bending relaxation times and bond lifetime.^53–57 The critical Casimir patchy particle chains, that have a persistence length between 250–150 particles for dT = 0.12–0.22 K (see Fig. 11 in ref. 33), are relatively stiff and thus show end-of-chain breakage.

Fig. 13 Passive and active decamer at dT = 0.12 K with F_A = 10 (left) and 100 (right) k_BT/σ. Plotted are the Stage I escape flux Φ_esc (see eqn (12)) (a) and (b), the Stage II separation probability P_sep (c) and (d) (see eqn (13)), and the breakage rate k_break (e) and (f) as a function of the bond index (see Fig. 2). Color and line coding in the legends at the bottom. Colored bands indicate the 95% confidence interval. The logarithmic plots handle the sometimes large scale differences between the curves. The passive decamer shows mostly breakage at the edges, while compressing and sliding active forces shift the breaking rate profile to bonds in the middle of the chain.

Extending forces. When applying extending forces to the decamer, the chain straightens, rendering lower S-values less likely (as can be observed in Fig. 12a). Simultaneously, the bonds at the center (bonds 2–5) are more stabilized compared to bonds at the end (bond 1), evidenced by the probability distribution of the bond energy − ln P(E_pair) (see Fig. 19 in Appendix B). This stabilization impacts stage I of breakage with F_A = 100, and to a lesser extent with F_A = 10: bond 2 to 5 diffuse slower out of the well compared to the passive chain as shown in Fig. 13a and b, respectively.

The stabilization of the chain and lowering of the rate of stage I would predict an increased lifetime. However, Fig. 14 indicates τ_sim exponentially decreases with the extending force magnitude F_A. It is the increase of stage II probability P_sep, that is tripled when F_A = 100, that is primarily responsible for the decrease in the decamer's lifetime (Fig. 13c and d).⁵⁸ The increase of the stage II probability is bond number independent, leaving the breakage profile in Fig. 13e and f largely unaffected compared to the passive decamer.

Fig. 14 The τ_sim of passive decamers (dT = 0.12–0.22 K) and active decamers (dT = 0.12, and 0.22 K) with extending, compressing, and sliding forces.

Compressing and sliding forces. Under compressing and sliding forces, the chain primarily buckles around the center and the second bond, respectively, as the distributions of S in Fig. 12b and c indicate.⁵⁹ The buckling leads to a strong destabilization of the chain, which in turn increases the rate of stage I at high activity (see Fig. 13b). At low activity (Fig. 13a), the effect on stage I is not as strong. The particle separation probability (stage II) is less severely impacted compared to stage I (Fig. 13a–d). Thus, the lowering of the (free) energy barrier of stage I is the main driving force of breakage and largely explains the breakage rate profile in Fig. 13e and f.

Weak interaction strength. When the attraction strength is weak (dT = 0.22 K), the architecture is not able to adjust its shape and mediate the effect of activity. These systems break fast and often near the active force. Data not shown.

3.3 Ring structures

Complex networks or gels, in contrast to freely moving chains, can contain tension even when energy minimized at zero temperature. The simplest architecture for which such tension can occur, is a ring structure. In our assessment, we consider two types of rings: a four TPP particle ring, which is under tension, and a six TPP particle ring, which is tension free (see Fig. 3). Additionally, we vary the length of the DP-chain from 5 to 15 particles, which can mediate the tension within the rings. We set the critical Casimir temperature to dT = 0.16 K, striking a balance between ensuring a reasonable simulation length until breakage and allowing the structures to remain intact long enough to observe the effects of activity on the deformations in the architectures.

The average bond lifetime of the rings is presented in Fig. 15. As anticipated, the presence of tension reduces the lifetime of the structures, and larger number of bonds leads breakage. The lifetimes τ_sim of the structures follow the order T4D5 ≪ T4D15 < T6D15 < T6D5. For all structures and both directions of active forces, the lifetime of the rings initially increases with increasing active force, but eventually decreases.

Fig. 15 The ring structures' lifetime τ_sim at dT = 0.16 K responds non-linearly to activity.

We use the T6D15 structure for an in-depth analysis of how the deformations (Fig. 16) and breakage rate are affected by the active forces (Fig. 17). This structure does not suffer from the misalignment of the TPP-DP bond nor from a limited number of bonds that can mediate the activity.

Fig. 16 Deformations in T6D15 at dT = 0.16 K. The probability histogram of S (eqn (2)) of a passive and active rings with F_A = 100k_BT/σ and the force directions: inward (a), and outward (b). The measurement is taken in region 1 + 2. The inner color of the dots indicate the position of the bond, and the edgecolor the passive (red) or active ring (black).

Fig. 17 Stage I (a) and (b), stage II (c) and (d), and breakage rate (e) and (f) per bond of T6D15 (a), (c) and (e), and of T4D15 (b), (d) and (f). Colored bands indicate 95% confidence interval.

Passive rings. While each bond is energetically equal, not all bonds are conformationally equal. The filled red circles in Fig. 16 indicate the TPP-DP bond (bond 1) in the ring visits lower values of S compared to bonds 2–8 (Fig. 16). The rate constant of stage I around bond 1, as well as stage II, is enhanced, making it the most preferred breakage position (Fig. 17).

Inward-facing forces. Under inward-facing forces, the dipatch particle chains buckle, as follows from the strong shift of the black curves in (Fig. 16a) toward lower S values. Remarkably, the buckling is distributed evenly over the bonds. The uniform buckling makes all bonds more prone to escape the potential well (see Fig. 17a). Stage II is affected differently, depending the magnitude of force (Fig. 17c). At low force (F_A = 10), only bond 1 is blocked from particle separation, similar to what is observed in the decamer under compressing forces. (Fig. 17e) shows this blocking leads to a flat breakage rate (and profile). At high force (F_A = 100), due to a leverage effect of the activity, the bonds close to the active force (bond 1–4) shows enhancement of particle separation, while bond 5–8 are suppressed. Overall it leads to enhanced breakage close to the focal point of the active force.

Outward-facing forces. Under outward-facing forces, the dipatch particle chains show a minor stabilization by evenly straightening the bonds, as indicated by the shift of the black curves to higher S values in Fig. 16b. Both large and small forces marginally suppress the first stage of a potential escape (see Fig. 17a), but have a more significant effect on the stage II probability P_sep. All bonds show enhanced probability of particle separation proportional to the force magnitude (Fig. 17c). For outward-facing forces, the second stage of bond breakage results in an enhancement of the breakage rate that is flat at small active force and is more pronounced at bonds close to the center of activity (i.e., low bond number) at large force.

Tension. Next, we will compare the T4D15 and T6D15 structures to identify the role of the misaligned TPP-DP bond on the response of activity using Fig. 17.

In the passive case and at low activity, we observe that the tension is evenly distributed over all bonds as T6D15 and T4D15 show the same shape of stage I and stage II versus bond number. However, T4D15 is faster in both stages.

At high activity, in the first stage with outward-facing forces, all bonds in T6D15 straighten and strengthen, while for T4D15, the straightening of the chain between the active particles expose the misalignment of the TPP-DP bond, making it prone to break, as shown in Fig. 17b, where Φ_esc at bond 1–2 is higher and bond 3–8 lower than the passive case.⁶⁰ Both structures shows similar increased stage II throughout all the bonds. Using inward-facing forces, bond 1, and 4 to 6 are faster in stage I in T4D15. However, in both structures, P_sep (stage II) of bond 4–8 is below 1% and the differences of stage I cannot be observed in the breakage rate k_break. Thus, for both force directions, the largest difference is observed at bond 1, which indeed is the misaligned TPP-DP bond.

Fluctuations. We performed analysis of the fluctuations in the intact ring structure before breakage, using principle component analysis (PCA) (see Appendix C). We find that activity mainly distorts the ring configuration, but leaves the dominant modes unaffected. Inward-facing forces have more influence than outward forces as they tend to cause buckling (see appendix for details).

4 Conclusions

In this work, we have given microscopic insight into the effect of activity, modeled as active patchy Brownian particles connected to various patchy particle architectures: dimers, decamers, and rings. To do so, we separated the contribution on the two stages of breakage: the rate of escaping the pair potential well (stage I) and the probability of successful particle separation (stage II). We unraveled the observed behavior in the bond breakage rate, breakage mechanism, and breakage profile and identify at which stage the influence of activity dominates.

Table 2 summarizes the effects of small and large active forces on dimers and decamers with F_A directions: extending, compressing and sliding.

Table 2 Summary of the breakage behavior of the dimer, decamer and rings as function of active force magnitude and direction. Dimer at dT = 0.12 K, mechanism: S (angular) r (radial). The breakage rate and breakage stages are qualitatively compared to F_A = 0 using arrows for: strong increase (⇑), increase (↑), approximately equal (− −), and decrease (↓) with respect to the passive case. For the decamer and rings, the breakage location is given as: edge: bonds close to the end of the chain, flat: uniform distribution, center: bonds in the middle of the chains, and w-shape: bonds halfway the middle

---

For the dimers, the effect of activity on the first stage of bond breakage resembles an equilibrium problem in an effective tilted potential along the radial direction, and the second stage of bond breakage was potential depth independent. Therefore, the effect of activity on the bonds' lifetime is roughly a function of the magnitude and direction of the active force only. When active forces point toward or away from the bonding volume, the bond life time τ increases or decreases, respectively.

The dimers and decamers showed qualitatively similar breakage mechanisms in both passive and active systems. Under extending (compressing) forces, the mechanism shifted to a radial (angular) route which is explained by the direction of the active force pointing away (toward) from the bonding volume. The sliding mechanism lies in between these two extremes.

The active decamers showed under all force directions and magnitudes reduced or comparable lifetimes compared to passive decamers. Extending forces straighten the chain, energetically stabilize the bonds, and suppress the rate constant of stage I. However, the chains still show an exponentially decaying lifetime proportional to the magnitude of the active force because simultaneously stage II is enhanced. For sliding and large compressing forces, the chains buckle and weaken their bonds leading to breakage peaked around the highest curvature in the chain. The weakened bond strength (stage I) dominates in explaining the effect of activity on the decreased lifetimes. However, at small compressing forces, stage II is suppressed at the ends of the chain and balances the effect of stage I. The active chain then shows an approximately equal breakage rate compared to passive chains.

All ring structures, with and without tension, show an initial stabilization of the structures upon small active forces (F_A = 10). The initial stabilization of the structure is caused by a minor suppression of stage I. Subsequently, for inward-facing forces, we observe an significant suppression of stage II in bonds close to the active force leading to a strong enhancement of the lifetime τ. However, for the outward-facing forces, the opposite is observed. The minor stabilization of stage I is followed by a minor enhancement of stage II leading to lifetimes approximately equal to the passive case.

At larger activity (F_A = 100), the bond lifetime τ has reduced significantly and converges for both directions of the forces to the same value. Additionally, this case shows most breakage near the active force at bond 1. However, the two directions of the force each target a different breakage stage. Inward facing forces induce buckling in the ring and weaken all bonds. Stage I shows an overall enhancement, but stage II is only enhanced at bonds close to the active force and suppressed at the bond in the center. Outward-facing forces show very little effect on the conformations of the ring as well as on stage I; it is stage II that has the largest effect on the breakage rate.

Lastly, we looked at the effect of the misaligned bond (bond 1) in the T4D15 structure compared to T6D15. At low activity and in the passive case, the tension caused by the misalignment can be distributed over all bonds. The lifetime of T4D15 is lower than T6D15, but there is no enhanced reactivity in bond 1; only at large forces, we observe differences. Outward- and inward-facing forces expose the weakness of bond 1 and enhance stage I of its breakage significantly. Simultaneously, as observed in T6D15, the second stage at bond 4 to 8 is very low as P_sep is less than 1%. Effectively, both directions of the force leads to enhanced breakage at the misaligned bond via stage I.

We can draw some general conclusions from our findings, summarized in Table 2. First of all, the preferred breakage mechanism is in all cases via an angular mechanism. This might at first seem surprising, since it seems only natural that a bond breaks along the radial axis. However, the entropic contribution of larger angles is substantially enlarged to counteract this, and most bonds will break at low S value, i.e. via the angular mechanism.

Second, it is a well-known fact that breakage is a combination of an activated well escape stage, and a diffusive separation stage. When applying an active force most of the time both stages will respond qualitatively identically: for the dimer extending and sliding force will enhance breakage rate, while compression will suppress it. For the decamer and the rings things are different: the extending/outward forces surprisingly lower the escape probability due to stiffening of the chain. However, the final result is still that breakage is enhanced, as the separation stage compensates. In contrast, sliding is significantly enhancing chain breakage. Compressing forces are most complex, they enhance stage I but decrease stage II. The contrasting effect of forces occurring in ring structure: inward force first stabilize rings, before they enhance breakage at higher forces.

It is clear that activity can both enhance and suppress breakage rates and mechanisms, depending on the precise magnitude and direction. We expect therefore that activity also will have a large effect on the dynamical behavior of network formation in colloidal patchy particle systems.

Our findings will thus be important for designing adaptive and responsive active physical gels, and rationalizing their behavior. In addition, if one wants to control breakage – and formation – of bonds one might target the different stages of breakage using activity. In a future study, we will report on how colloidal gels self-assembled from patchy particle mixtures respond to activity. On the other hand, one might extend the current investigation to other architectures, such as stars, multiple rings, and 3D structures.

The predictions made in this work can be experimentally tested by using particle systems that behave according to the potential used here.³³ Using microscopy the breaking behavior can be followed and tracked, see e.g. ref. 61 and 62.

While our predictions were made specifically for critical Casimir interactions, we believe that they can be generalized to a large class of patchy interaction potentials, that are qualitatively similar: strong short range attractions induced by a relatively small patch in combination with ABP dynamics. As long as the rotational and translational diffusion constants behave similarly to those of the colloids we consider here, and the ABP model applies to the self-propulsion mechanism, the resulting breakage dynamics should be qualitatively similar, irrespective of the origin of the attraction and the self-propulsion.

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. 33.

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. 33. 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 6 in ref. 33.

A.2 External gravitational potential. The total potential is a summation over the pair potential and an external potential (eqn (3)). The external gravitational potential is composed of two terms: a hard wall represented by a steep Lennard-Jones potential V_LJ and a gravitational potential V_g that depends, among other things, on the mass of the particle. See the Appendix B2 of ref. 33 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.

B Decamers' breakage mechanism

Fig. 18 presents the P_1→3, P_2→4, and P_mech of passive and active dimers and decamers at dT = 0.12 K and F_A = 100k_BT/σ.

Fig. 18 Analysis of the breakage mechanism by P_1→3, P_2→4, and P_mech from eqn (14)–(16), respectively. Dimers (a)–(c) and decamers (d)–(f) with dT = 0.12 K and active force directions with extending, compressing, and sliding active forces at F_A = 100k_BT/σ. Note that breakage mechanism and both stages are in principle bond number dependent, but here we take the average over the whole chain.

Stage I of breakage in a decamer is affected by the shape of the chain. As shown in Fig. 12 by the distribution of S per bond, and in Fig. 19 the distribution of E_pair, under activity the decamer changes its shape. This shape change: straightening under extending forces, and buckling under compressing and sliding forces, is reflected in the P_1→3 of Fig. 18d as a shift toward high and low S meaning radial and angular breakage, respectively.

Fig. 18b and e show that the effect of activity on probability of successful particles separation in the decamer is lower than in the dimer, especially at S < 0.5. The particles in the chains are not freely diffusing in region 3, as they may be still bound to their neighbours. This effect is explicitly measured in Fig. 13c, which shows that the ends of the chain separate twice as successful compared to the center.

Comparing the breakage mechanism of the dimer and decamer (Fig. 18c and f), we observe qualitative similarities. The passive case it is practically similar, the active case shows an exaggeration of the shift in breakage mechanism. Extending forces: no preferred breakage route, and compressing and sliding forces: angular breakage.

Fig. 19 Passive and active decamers. The probability histogram of E_pair (eqn (1)) of a passive chain (red solid line) and active chain with F_A = 100k_BT/σ (dots) and the force directions: extending (a), compressing (b), and sliding (c). Due to symmetry, only the first to fifth bond are shown with symbols as depicted with the colored dots in the legend of Fig. 12. The inset shows a zoom out of the distributions. The measurement is taken at E_pair < 0 k_BT.

C Fluctuations

To be able to understand the response of the ring structures to the introduction of activity, we can analyze their dynamical fluctuations using principle component analysis (PCA).⁶⁴ We perform PCA on long straightforward Brownian dynamics simulation trajectories of rings using the MDAnalysis software package.^65,66 We focus on the ring system with 4 TPP and 15 DP particles (T4D15), as this system contains long chains that allow significant buckling, but is smaller than the 6 TPP structures and thus decorrelates significantly faster.

We start with the analysis of the influence of active forces on the fluctuation dynamics and deformations of an elastic ring using principal component analysis (PCA), as depicted in Fig. 20 and reported Table 3. To enhance the likelihood of observing intact structures during the fluctuation and deformation simulations, we set the interaction strength to dT = 0.12 K.

Fig. 20 PCA: the particles' average position (a) and (b), and the two dominating fluctuation modes: first mode (c) and (d), and second mode (e) and (f) of passive and activated T4D15 rings. For clarity, the eigenvectors of the two fluctuation modes are illustrated with 10× magnified arrows for clarity. Each column has its legend shown at the top.

Table 3 The contributions to the variance of the first two modes (shown in Fig. 20) retrieved via PCA. Column wise the magnitude and direction of the active force, the loading of the variance of modes 1 and 2 and their sum are shown

F A
Magnitude	Direction	Mode 1	Mode 2	Sum
a In the trajectory of inward-facing FA = 15 forces, buckling and breaking was observed.
0		0.55	0.26	0.81
2	In	0.59	0.25	0.84
5	In	0.56	0.31	0.87
15^a	In	0.78	0.11	0.89
2	Out	0.55	0.27	0.82
5	Out	0.49	0.28	0.77
15	Out	0.50	0.25	0.75
25	Out	0.47	0.20	0.67

---

We conducted Principal Component Analysis (PCA) on long trajectories, with all simulations lasting at least 30 000 seconds. However, for an inward-directed active force F_A = 15, the simulation duration was reduced to 5000 seconds due to observed buckling and a high probability of breaking.

The eigenvectors corresponding to the two dominant modes, which account for approximately 80% of the variance, are depicted in Fig. 20, and their respective eigenvalues are listed in Table 3. These dominant modes are particularly prevalent since they correspond to the lowest curvature in the semi-flexible polymers, resulting in minimal energy cost.

For small outward-directed forces up to F_A ≤ 25 and inward-directed forces up to F_A ≤ 5, no significant deformations were observed, and the contribution to the variance was not notably different from the passive case (see Fig. 20a, c and e). Any observed differences between the systems can be attributed to limited sampling. Conversely, larger inward-facing active forces (F_A ≥ 15) caused the DP-particle chains to buckle, as evidenced by the average positions of these trajectories (Fig. 20b). The buckling is caused by the movement of the first mode (Fig. 20d) which indeed shows a significantly higher contribution compared to the passive case (Table 3).⁶² In this deformed structure, the two dominant modes still exist, albeit with different amplitudes. Note that an even higher inward force inevitably led to buckling and breaking.

Since the active force, whether inward or outward, primarily affects the radial part of the potential through the translational displacement of particles, it is not expected to have a substantial impact on the fluctuation spectrum. The fluctuations primarily depend on the bending rigidity, as defined by the shape of the S′-function. Moreover, the active force is fixed to the particles and constrained by two bonds. Consequently, the active force cannot reorient along fluctuation modes to drive them, as observed in the actuated elastic solid in ref. 67.

Author contributions

HJJ – equal: conceptualization, formal analysis, investigation, methodology, writing – original draft, writing – review and editing. Lead: data curation, software, validation, visualization. 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

The authors have no conflicts to disclose.

Acknowledgements

The authors acknowledge the funding (Grant No. 680.91.124) from the Foundation for Fundamental Research on Matter (FOM), which is part of The Netherlands Organization for Scientific Research (NWO).

Notes and references

1. J. Prost, F. Jülicher and J.-F. Joanny, Nat. Phys., 2015, 11, 111–117.
2. E. H. Barriga, K. Franze, G. Charras and R. Mayor, Nature, 2018, 554, 523–527.
3. F. Burla, Y. Mulla, B. E. Vos, A. Aufderhorst-Roberts and G. H. Koenderink, Nat. Rev. Phys., 2019, 1, 249–263.
4. G. L. Jackson, J. M. Dennis, N. D. Dolinski, M. van der Naald, H. Kim, C. Eom, S. J. Rowan and H. M. Jaeger, Macromolecules, 2022, 55, 6453–6461.
5. M. P. Murrell and M. L. Gardel, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 20820–20825.
6. J. Alvarado, M. Sheinman, A. Sharma, F. C. Mackintosh and G. H. Koenderink, Nat. Phys., 2013, 9, 591–597.
7. 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.
8. 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.
9. H. N. Verwei, G. Lee, G. Leech, I. I. Petitjean, G. H. Koenderink, R. M. Robertson-Anderson and R. J. McGorty, J. Vis. Exp., 2022, 1–27.
10. R. G. Winkler, J. Elgeti and G. Gompper, J. Phys. Soc. Jpn., 2017, 86, 101014.
11. D. Rogel Rodriguez, F. Alarcon, R. Martinez, J. Ramírez and C. Valeriani, Soft Matter, 2020, 16, 1162–1169.
12. S. Joo, X. Durang, O.-C. Lee and J.-H. Jeon, Soft Matter, 2020, 16, 9188–9201.
13. A. Militaru, M. Innerbichler, M. Frimmer, F. Tebbenjohanns, L. Novotny and C. Dellago, Nat. Commun., 2021, 12, 2446.
14. J. Zhang, T. Huang, G. Xu and Y. Chen, Commun. Theor. Phys., 2022, 74, 075601.
15. Y. Kim, S. Joo, W. K. Kim and J.-H. Jeon, Macromolecules, 2022, 55, 7136–7147.
16. X. Xia, H. Hu, M. P. Ciamarra and R. Ni, Sci. Adv., 2020, 6, 1–8.
17. M. R. Jones, N. C. Seeman and C. A. Mirkin, Science, 2015, 347, 1260901.
18. 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.
19. Z. Gong, T. Hueckel, G.-R. Yi and S. Sacanna, Nature, 2017, 550, 234–238.
20. R. Khalaf, A. Viamonte, E. Ducrot, R. Mérindol and S. Ravaine, Nanoscale, 2023, 15, 573–577.
21. Y. Shelke, S. Marín-Aguilar, F. Camerin, M. Dijkstra and D. J. Kraft, J. Colloid Interface Sci., 2023, 629, 322–333.
22. T. Nguyen, A. Newton, D. Kraft, P. G. Bolhuis and P. Schall, Materials, 2017, 10, 1265.
23. 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.
24. 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.
25. 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.
26. I. Chakraborty, D. J. G. Pearce, R. W. Verweij, S. C. Matysik, L. Giomi and D. J. Kraft, ACS Nano, 2022, 16, 2471–2480.
27. P. J. M. Swinkels, Z. Gong, S. Sacanna, E. G. Noya and P. Schall, Nat. Commun., 2023, 14, 1524.
28. J. R. Howse, R. A. L. Jones, A. J. Ryan, T. Gough, R. Vafabakhsh and R. Golestanian, Phys. Rev. Lett., 2007, 99, 048102.
29. H. Seyforth, M. Gomez, W. B. Rogers, J. L. Ross and W. W. Ahmed, Phys. Rev. Res., 2022, 4, 023043.
30. G. S. Redner, A. Baskaran and M. F. Hagan, Phys. Rev. E, 2013, 88, 012305.
31. G. S. Redner, C. G. Wagner, A. Baskaran and M. F. Hagan, Phys. Rev. Lett., 2016, 117, 148002.
32. C. F. Lee, Soft Matter, 2017, 13, 376–385.
33. H. J. Jonas, S. G. Stuij, P. Schall and P. G. Bolhuis, J. Chem. Phys., 2021, 155, 034902.
34. I. M. Ilie, W. K. den Otter and W. J. Briels, J. Chem. Phys., 2015, 142, 114103.
35. S. G. Stuij, PhD thesis, University of Amsterdam, 2020.
36. K. Dietrich, D. Renggli, M. Zanini, G. Volpe, I. Buttinoni and L. Isa, New J. Phys., 2017, 19, 065008.
37. N. Oikonomeas-Koppasis, S. Ketzetzi, D. J. Kraft and P. Schall, Power-law intermittency in the gradient-induced self-propulsion of colloidal swimmers, 2023, https://arxiv.org/abs/2310.17384.
38. B. Szabó, G. J. Szöllösi, B. Gönci, Z. Jurányi, D. Selmeczi and T. Vicsek, Phys. Rev. E, 2006, 74, 061908.
39. S. van Teeffelen and H. Löwen, Phys. Rev. E, 2008, 78, 020101.
40. S. Henkes, Y. Fily and M. C. Marchetti, Phys. Rev. E, 2011, 84, 040301.
41. L. Fang, L. Li, J. Guo, Y. Liu and X. Huang, Phys. Lett. A, 2022, 427, 127934.
42. A. C. Newton, J. Groenewold, W. K. Kegel and P. G. Bolhuis, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 15308–15313.
43. A. Vijaykumar, P. R. ten Wolde and P. G. Bolhuis, Mol. Phys., 2018, 116, 3042–3054.
44. T. S. van Erp, D. Moroni and P. G. Bolhuis, J. Chem. Phys., 2003, 118, 7762–7774.
45. G. Cowan, Statistical data analysis, Oxford University Press, USA, 1998.
46. H. Flyvbjerg and H. G. Petersen, J. Chem. Phys., 1989, 91, 461–466.
47. A. Vijaykumar, P. G. Bolhuis and P. R. ten Wolde, Faraday Discuss., 2016, 195, 421–441.
48. J. K. G. Dhont, G. W. Park and W. J. Briels, Soft Matter, 2021, 17, 5613–5632.
49. H. Kramers, Physica, 1940, 7, 284–304.
50. P. Hänggi, P. Talkner and M. Borkovec, Rev. Mod. Phys., 1990, 62, 251–341.
51. A. Geiseler, P. Hänggi and G. Schmid, Eur. Phys. J. B, 2016, 89, 175.
52. L. Natali, L. Caprini and F. Cecconi, Soft Matter, 2020, 16, 2594–2604.
53. J. Paturej, A. Milchev, V. G. Rostiashvili and T. A. Vilgis, J. Chem. Phys., 2011, 134, 224901.
54. A. Zaccone, I. Terentjev, L. Di Michele and E. M. Terentjev, J. Chem. Phys., 2015, 142, 114905.
55. C. F. Lee, J. Phys.: Condens. Matter, 2015, 27, 275101.
56. C. F. Lee, J. Phys.: Condens. Matter, 2018, 30, 315102.
57. M. Razbin, P. Benetatos and A. A. Moosavi-Movahedi, Soft Matter, 2019, 15, 2469–2478.
58. A. Ghosh, D. I. Dimitrov, V. G. Rostiashvili, A. Milchev and T. A. Vilgis, J. Chem. Phys., 2010, 132, 204902.
59. A. Mondal and G. Morrison, J. Chem. Phys., 2022, 157, 104903.
60. A. M. Lorenzo, E. M. De La Cruz and E. F. Koslover, Soft Matter, 2020, 16, 2017–2024.
61. S. Stuij, J. M. van Doorn, T. Kodger, J. Sprakel, C. Coulais and P. Schall, Phys. Rev. Res., 2019, 1, 023033.
62. S. G. Stuij, A. Biebricher, Z. Gong, S. Sacanna, E. Peterman, I. Heller and P. Schall, Phys. Rev. Mater., 2022, 6, 035603.
63. V. A. Parsegian, van der Waals Forces, Cambridge University Press, 2005, pp. 1–380.
64. S. Henkes, C. Brito and O. Dauchot, Soft Matter, 2012, 8, 6092.
65. N. Michaud-Agrawal, E. J. Denning, T. B. Woolf and O. Beckstein, J. Comput. Chem., 2011, 32, 2319–2327.
66. R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, J. Domanski, D. L. Dotson, S. Buchoux, I. M. Kenney and O. Beckstein, Proceedings of the 15th Python in Science Conference, 2016, pp. 98–105.
67. P. Baconnier, D. Shohat, C. H. López, C. Coulais, V. Démery, G. Düring and O. Dauchot, Nat. Phys., 2022, 18, 1234–1239.

---