Necking and failure of a particulate gel strand: signatures of yielding on different length scales

Abstract

“Sticky” spheres with a short-ranged attraction are a basic model of a wide range of materials from the atomic to the granular length scale. Among the complex phenomena exhibited by sticky spheres is the formation of far-from-equilibrium dynamically arrested networks which comprise “strands” of densely packed particles. The aging and failure of such gels under load is a remarkably challenging problem, given the simplicity of the model, as it involves multiple length- and time-scales, making a single approach ineffective. Here we tackle this challenge by addressing the failure of a single strand with a combination of methods. We study the mechanical response of a single strand of a model gel-former to deformation, both numerically and analytically. Under elongation, the strand breaks by a necking instability. We analyse this behaviour at three different length scales: a rheological continuum model of the whole strand; a microscopic analysis of the particle structure and dynamics; and the local stress tensor. Combining these different approaches gives a coherent picture of the necking and failure. The strand has an amorphous local structure and has large residual stresses from its initialisation. We find that neck formation is associated with increased plastic flow, a reduction in the stability of the local structure, and a reduction in the residual stresses; this indicates that the system loses its solid character and starts to behave more like a viscous fluid. These results will inform the development of more detailed models that incorporate the heterogeneous network structure of particulate gels.

---

1 Introduction

Even well-understood basic models can sometimes give rise to strikingly complex phenomena. Among these are so-called “sticky spheres”, which are hard spheres with a short-ranged attraction. They represent arguably the smallest increase in complexity from hard spheres, the simplest non-trivial model of atomic and colloidal systems.¹ Although their equilibrium phase behavior has been, to a considerable extent, understood since the 1960s,² out of equilibrium, sticky spheres exhibit a bewildering range of phenomena, including multiple glassy states³ and gelation, which is a complex blend of phase separation and dynamical arrest leading to a far-from-equilibrium unstable state.⁴

Such gelation caused by arrested phase separation can occur in a wide range of materials across a wide range of length scales, from phase-demixing oxides⁵ and metallic glassformers,⁶ small molecules,⁷ clays,⁸ and granular matter.⁹ Among the most well-known gel-forming systems are soft materials, including proteins,¹⁰ foods,¹¹ hydrogels,¹² tissues¹³ and perhaps most studied of all, colloids.^4,14

Gels typically consist of heterogeneous networks of connected “strands”,¹⁵ which are dynamically arrested. They exhibit complex phenomena that continue to resist scientific understanding, including complex aging behavior, and the possibility of self-induced catastrophic failure after an extended and unpredictable waiting time.^4,14 During the waiting time, the gels actually become stronger before undergoing catastrophic failure.¹⁶ Given the simple nature of the constituents, and the long-understood equilibrium behavior, the non-complexity of the non-equilibrium behavior is quite remarkable.

From empirical observation^{14,15,17–23} we know that the main control parameters for gels of sticky spheres are the strength of attractive interactions between particles, and the particle volume fraction. These parameters influence the mechanical properties of the gel, including its elastic response to small mechanical perturbations. When larger forces are applied, the strands of the gel network break, leading to macroscopic flow.^4,17,24–31 An important deformation mode in this process is stretching and eventual breakage of strands under tensile loading, leading to a stress reorganization of the surrounding gel, and hence to flow.³² (This is different from polymer gels, where entanglements dominate the rheology.) An important signature of the non-equilibrium gel state is that its properties – such as elastic moduli and yield stress – depend on its history, including its age.^{16,19,24,33–36}

A fundamental open question concerns the mechanical failure of particulate gels: how can the microscopic interactions between the particles be used to control the macroscopic failure mechanism? This question is challenging because mechanical failure is sensitive to physical processes on several length scales,³⁷ from individual particle diameters,³⁸ to the strand thickness, up to the size of a macroscopic gel.^{19,23,39–41} There is a corresponding range of time scales:^27,42 in large gel samples, there may be a long period of aging, after which the gel quickly collapses.¹⁶ This situation makes gels difficult to study. Computer simulations can follow the motion of the particles, but the computational cost of simulating the long aging period is prohibitive, especially given the large (macroscopic) system sizes that are required. Experiments suffer from similar problems: techniques for monitoring a macroscopic sample do not allow individual particles' motion to be resolved, while particle-resolved techniques are limited to moderate time scales and cannot simultaneously resolve the behaviour throughout a macroscopic sample. For theoretical analysis, continuum modeling approaches are available, but these do not resolve behaviour of individual particles, relying instead on constitutive models. Nor do they treat the heterogenous nature of the material fully.

This situation calls for an approach on multiple scales, in order to combine the useful aspects of different methods. In this work, we address the fundamental events which are expected to lead to gel failure. We focus on the breakage of an individual strand, as a fundamental process that feeds into macroscopic yielding,^26,39 and into aging (or coarsening) of gels.⁴³ Like many processes in material failure, this yielding is challenging to predict and control, even with state-of-the-art theories. For example, the thickness of the strand ranges from 1 to 10 colloidal diameters, so single-particle fluctuations can have significant impact on the whole strand: this limits the predictive power of continuum models. Also, the strength of attractive interactions is only a few times the thermal energy, so there are frequent fluctuations in which interparticle bonds are broken and re-formed, necessitating a statistical mechanical approach. In addition, the strand is itself an amorphous material, so subtle features of the particle-level structure can have significant effects on its large-scale behavior,³⁴ as is familiar in glassy materials.

In response to these challenges, we analyse the failure of a single strand, under elongation. We performed computer simulations of this process, analysing the results using three complementary methods. These are a simple continuum modelling approach based on previous established models;^44–46 a new characterisation at the single-particle level; and measurements of the local stress.⁴⁷ These analyses lead to a coherent and unified view of the failure process over different length scales. Failure occurs by a necking mechanism, which proceeds via a feedback mechanism, leading to a linear instability. Such instabilities may generally occur in several different ways:^44–46 Our results confirm that the tensile force in the strand generates an increased stress in the neck, resulting in increased plastic flow there, which causes further thinning. This increases the local stress even more, and so on, until the strand breaks. In simple terms, one may imagine that the macroscopic stress σ exceeds the yield stress σ_Y in the neck,^28,38 while remaining below σ_Y elsewhere. Plastic flow in the neck is revealed by local measurements of increased particle motion, and this is coupled with a reduction in the number of low-energy (stable) structures.

We refine this established picture in two ways. First, a top-down continuum modelling approach⁴⁶ indicates the existence of an internal time-dependent plasticity field that determines the response to local stress. Second, particle-level measurements of the local stress⁴⁷ reveal a complex pattern of residual stresses^48–50 that come from the initialisation of the strand;⁵¹ this pattern changes significantly in the neck, as it develops. This is reminiscent of other amorphous systems where mechanical heterogeneous properties can determine plasticity.^{28,31,35,52,53} These results bridge the scales between continuum modelling (rheology) and the particle level (local motion and local structure), allowing new relationships to be revealed. The connections are mediated by our direct measurements of local stress.

By combining these analyses on different scales, our results greatly extend previous work on single gel strands^41,54 and on breakage of glassy (amorphous) samples,^46,55 and other soft materials.⁵⁶ Some other approaches to gel modelling^57,58 consider the gel as a continuum, without resolving individual strands – our approach can connect such models with the microscopic gel structure. Similarly, our detailed analysis of individual strands complements alternative simulation approaches where particle interactions are justified by a top-down approach:^23,32 these models do not resolve the microscopic structure of the strands, but they do enable simulations of an entire gel sample. The depletion attraction in these gels is a few times larger than the thermal energy k_BT. This leads to gels with thick strands, in which significant plastic deformation occurs before breakage.⁵⁹ This is different from dilute gels with thin strands, where there can be visco-elastic behaviour in the overall gel, but individual strands have a sudden brittle failure.²³

Overall, our results elucidate the microscopic mechanisms for strand failure. They include the discovery of new relationships between emergent properties (like yielding) and local (microscopic) structure, showing how to bridge the gaps between microscopic and macroscopic mechanical interactions in arrested out-of-equilibrium systems. An understanding of these relationships, which can serve as a foundation for the coarse-grained modelling of macroscopic gels, is vital for the design and control of gels, as an important class of soft materials.

2 Results

2.1 Overview

We briefly describe our model system, with full details in Section 4. We simulate a gel-forming system^60,61 of particles with short-ranged attractive interactions, and two different particle sizes, to suppress crystallisation.⁶¹ Particle i has position r_i and velocity v_i = ṙ_i, it evolves by Langevin dynamics in a simulation box with periodic boundaries:where m is the particle mass, is the interaction energy, λ₀ is a friction constant (proportional to the solvent viscosity), and F_solv the random solvent force, whose variance is proportional to the ambient temperature, see Section 4.1 for further details. The mean particle diameter is denoted by l.

This system forms gels by arrested spinodal decomposition.⁶¹ A snapshot of one such gel is shown in Fig. 1(a), showing a network of strands that evolves in time. In this work, we initialise our system to mimic a single strand of such a gel, by quenching a bulk colloidal liquid to form an amorphous (glassy) solid, and excising a cylindrical sample, see Fig. 1(b). The strand has initial length L_‖ = 30l and we vary its initial radius r₀ to mimic strands of different thicknesses. An important time scale is the Brownian time τ_b (see Section 4.1), which is the typical time required for an individual particle to diffuse its own radius. Before excising a strand we allow the high-density mixture to equilibriate (for 10⁴τ_b) before excising it. We then also allowed the surface of the excited strand to relax for 10³τ_b before we started to deform the simulation box. We varied this equilibration time up to 10⁵τ_b and only found a weak dependence on the age of the gel (see Section 4.2). If this aging time was comparable with the relaxation time of the high-density mixture, the nature of the failure might be different, but we do not consider this regime here.

Fig. 1 (a) A 2D slice of the interior of a 3D gel simulation, showing part of the simulation box containing a thick strand (red) and a strand in the process of breaking (turquoise). (b) A gel strand after initialization (γ = 0) with uniform thickness, red arrows show the direction of elongation. Parameters: , r₀ = 4l and [small gamma, Greek, dot above] = 5 × 10⁻⁶. (c) The same strand at strain γ = 0.4, showing formation of a neck. (d) and (e) The number of particles in the thinnest region of the strand N_m and the tensile force in the strand F^z, as a function of strain γ. Thin lines show 8 runs from the same starting conditions (thin blue lines). Thick lines are an average over 100 runs. (f) The cross-sectional stress [small sigma, Greek, macron]^zz, averaged over z, for different varying strain rates [small gamma, Greek, dot above]. (Shaded regions indicate the standard deviation over 6 separate runs.) The dashed lines in (e) and (f) are homogeneous solutions of the continuum theory (eqn (3) and (4)), which are valid for times before significant necking has occurred (see main text). (g) Distribution of the position of the neck, for 100 runs with identical initial conditions. Dashed line indicates the thinnest segment at γ = 0.

We deform the strand by affine elongation of the simulation box, so the length L_‖(t) increases with time t. This stretches the strand, which eventually breaks. Fig. 1(b and c) illustrates the resulting behavior, which is governed by four dimensionless parameters (see Methods 4.1): the attraction strength between the particles a rescaled elongation rate [small gamma, Greek, dot above], the strand thickness r₀/l, and a solvent damping parameter λ. Gelation is associated with metastable gas–liquid phase separation,¹⁴ whose critical point is at .⁶² The data of Fig. 1 have , corresponding to , well inside the regime of gelation. We fix λ = 10 throughout, consistent with a gel in a high-friction solvent environment. For particulate gels that form in the absence of any solvent,⁷ one would consider instead the limit of small damping; we would expect similar results in this case too. For the other parameters, we mostly focus on the representative values used in Fig. 1(b–e). The qualitative picture is robust to changing these parameters, this will be discussed below.

Fig. 1(c) illustrates the neck that forms as the strand is stretched. It becomes increasingly thin and eventually pinches off, so the strand breaks. This behaviour is illustrated in Movies S1–S3 (ESI†). To focus on this effect, we divide the simulation box into 20 segments along the z-direction, indexed by their rescaled positions Z = z/L_‖(t). We count the number of particles in each segment, and write N_m for the smallest such number, which serves as a proxy for the thickness of the neck. Fig. 1(d) shows N_m, as the strand is elongated. To illustrate the variability in the failure mechanism, we show results for 8 representative simulation runs, all from the same initial structure, with different random forces, (such isoconfigurational ensembles have been used before in studies of glassy materials⁶³). We also show the average of N_m, obtained over 100 such runs. The formation of the neck is clear, as is the gradual reduction in its thickness, leading to failure. This corresponds to a ductile failure mechanism: all the amorphous strands considered here break in this way, we do not see any signs of brittle behaviour.

We measure the stress in the strand by the method of Irving and Kirkwood (IK).^47,64–67 The resulting stress estimates are noisy, due to rapid thermal fluctuations, so all such measurements are averaged over a time period of 250τ_b, to reduce the statistical uncertainty. The effects of this averaging are discussed in Appendix A.3. Throughout this work, all stresses (and elastic moduli) are quoted in units of k_BT/l³.

The strand is under tension and the corresponding tensile force is easily obtained from the total (global) stress: its average is shown in Fig. 1(e), together with several representative trajectories. There are three clear regimes, as the strand is elongated. Initially, the response is elastic and the stress increases. This is followed by a short plateau – a signature of plastic flow – before the stress smoothly decreases, until breakage occurs. This decrease mirrors the reduction in the neck thickness [Fig. 1(d)]. Physically, the force required to maintain a constant strain rate decreases as the neck gets thinner. However, there is no sudden drop in stress, as would be expected for brittle failure. In all measured cases, the neck gets gradually thinner, until the entire tensile force has to be sustained by a cross-section containing only 1–2 particles. At this point the tension is very small, and the strand ruptures.

While mechanical stability requires that the tensile force F^z is constant along the arm, the IK method also provides a local measurement of the stress tensor σ(r), which reveals additional information about the elongation process. To make contact with continuum models of rheology,^44–46 we consider the zz component of the stress, averaged over the cross-section of the strand (numerical estimation of this quantity is discussed in Section 4). The resulting quantity is denoted by [small sigma, Greek, macron]^zz(Z), it is related to the force as [small sigma, Greek, macron]^zz(Z) ≈ F^z/A(Z) where A(Z) is the cross-sectional area. After averaging the stress over Z, we plot the result in Fig. 1(f), for several different values of the elongation rate [small gamma, Greek, dot above]. These results are compared with a simple rheological model (dashed lines), see below for details.

Recall that we performed multiple simulation runs from the same initial configuration. For this initial condition, Fig. 1(g) shows that the neck is more likely to occur in particular locations. In the simplest theories of this instability, the neck should always form at the thinnest point on the strand. This is not consistent with the data, which is a first indication that the internal structure of the arm is playing a role in the rheology, as predicted by more advanced theories.^45,46

Before analysing these effects in more detail, we identify a surprising aspect of neck formation: Fig. 2 shows that the particles in the neck just before breakage (coloured blue) have arrived in the neck region from a range of other locations in the strand. See also Movie S1 (ESI†). This demonstrates significant mobility of particles, especially for those near the surface of the strand, and for those in the neck (see also Fig. 6 and Section 4.5, below).

Fig. 2 Snapshots of the strand from a single trajectory, as the neck forms. (a) Strain γ = 0.06, (b) γ = 0.20, (c) γ = 0.36. The 16 particles that form the neck in (c) are coloured in blue, to show their movement.

2.2 Continuum-scale description

We develop a simple description of the necking process at the continuum level. On this scale, two distinct types of necking instability in amorphous cylindrical strands have been identified:^44–46 gradual plastic deformation or sudden failure. The behaviour found in our system corresponds to the plastic (gradual) case. (The sudden mechanism – not seen here – is driven by the build up of elastic stress).

Recall that Z = z/L_‖(t) is a rescaled co-ordinate along the strand. Modeling this strand as a thin filament, the theory^44–46 is based on four Z-dependent quantities: the filament's cross-sectional area a(Z); the local strain rate [small gamma, Greek, dot above]_L(Z); a scalar field W(Z) that coincides (in this case) with the stress [small sigma, Greek, macron]^zz(Z); and an internal field χ(Z) that controls the plastic flow rate (see below). The necking instability can be captured by the following toy model. We follow the model proposed in the literature,⁴⁶ but with the simplification that stress can be captures with a scalar field⁴⁵ and a simplified internal field relaxation.

The elongation rate is slow enough that advective effects are negligible during elongation, and local force balance holds at all times. Then conservation of mass requires ∂_ta = ([small gamma, Greek, dot above] − [small gamma, Greek, dot above]_L)a, and the force balance condition is div σ = 0 (see Appendix C). For the filament, this means that

We identify aW as the tensile force in the strand, which is constant along its length.

Our equation of motion for W assumes that stress increases due to elastic loading and relaxes by plastic flow:

where G is the elastic modulus and p is the rate of plastic relaxation,⁴⁶ which depends on the stress and the plasticity χ. We assume that plastic flow only takes place above the yield stress σ_Y, taking p(W,χ) = τ_p(χ)⁻¹[(W/σ_Y) − 1]θ(W − σ_Y) where τ_p is a plastic time scale and θ is the Heaviside (step) function. Following ref. 46, the plasticity χ is analogous to a local temperature that determines the probability of an activated rearrangement, inspired by the theory of shear transformation zones.⁶⁸ Hence, τ_p(χ) = τ_ref exp(1/χ) where τ_ref is a reference time scale.

The general approach of ref. 44 and 45 starts with homogeneous solutions to the equations of motion (that is, a, [small gamma, Greek, dot above]_L, W, χ independent of Z). Necking is a linear instability of this solution. For the homogeneous solution, we first assume that the field χ is a simple constant, independent of both t and Z. One finds W(t) = G[small gamma, Greek, dot above]t + W(0) for short times (such that W < σ_Y) while for long times W(t) = σ_Y/(1 − [small gamma, Greek, dot above]τ_p). This theory was used to fit the data of Fig. 1(e), in the early-time regime before the neck forms, and is only valid for small strains where we can assume the homogeneous solution of the hydrodynamic equations. However, the data for faster elongation rates in Fig. 1(f) shows a stress overshoot, which cannot be described at this level of theory.

To account for this effect, we allow the plasticity field χ to depend on time with a relaxational dynamics that is controlled by the plastic relaxation itself:⁴⁶

The additional parameter χ_∞ fixes the steady-state value of χ. This model (with time-dependent χ) was used to fit the stress in Fig. 1f: there is good agreement, given the simplicity of the model. (As noted above, it is essential that χ is time-dependent, otherwise the model cannot capture the stress overshoot at the largest strain rate).

This model supports two main regimes. For very slow elongation with [small gamma, Greek, dot above]τ_p ≪ 1, the field χ is constant until W reaches the yield stress, at which point χ relaxes quickly to χ_∞. The slowest elongation rate in Fig. 1(f) is consistent with this regime, fitting gives σ_Y ≈ 3.3 (recall the units are k_BT/l³). For faster elongation rates, the relaxation of χ competes with the elongation rate, this is responsible for the stress overshoot in Fig. 1(f). Fitting simultaneously to all the curves in Fig. 1(f), we estimate G = 420, χ_∞ = 80 and τ_ref = 180τ_b, as well as the initial condition χ(t = 0) = 0.65, see Appendix C for further details. We note that the essential physical features that are required to fit the data are the existence of a yield stress that depends weakly on [small gamma, Greek, dot above] (which ensures a stress plateau as seen in Fig. 1(f)), and a non-trivial time dependence of χ (which accounts for changes in structure on elongation and allows the model to capture the stress overshoot). The model proposed here is just one way to capture these features, and there are other models in the literature that capture similar features.⁴⁵ More detailed atomistic simulations of amorphous solids have been done in the past with success to get relevant parameters for continuum models based on shear transformation zones.^69,70 In the present context, it is not obvious a priori that any kind of continuum description is relevant for these very thin strands, which are only a few particles in diameter. Hence we focus on a simple illustrative model, which can capture the main features of necking and failure.

Having characterised the homogeneous solutions for W, we can now use this model to perform a linear stability analysis by allowing Z-dependent perturbations, see Appendix C. The results indicate that the system is always unstable once W rises above σ_Y, so necking should occur in all cases, as observed.

A simple physical picture of this instability is that since the tensile force is constant along the strand [eqn (2)], the stress is largest at its narrowest point: this will be the location where W first exceeds σ_Y, leading to plastic flow near this point, and hence to further thinning. However the observation of Fig. 1(g) – that the neck does not always form at the thinnest point – indicates that heterogeneities in the internal structure of the strand also play a role in the instability. These features enter the model as inhomogeneities in χ, enabling the theory to explain the decoupling of the neck position and the initial thickness. In addition to this role of internal structure in determining the neck position, note that thermal fluctuations are also relevant because they give rise to the non-trivial distribution in Fig. 1(g). Such fluctuations are suppressed as increases and the system approaches the athermal limit.

2.3 Particle-level description

Having analysed the behaviour at the level of the entire strand and found that they match well regardless of the large fluctuations on these small scales, we turn to the microscopic (particle-level) structure.

We begin with a comparison of dynamical quantities. We define a correlation function C_B(Z;γ,Δγ) that measures how much particles' local environments have relaxed in segment Z, for the period between strain γ and strain γ + Δγ (see Section 4.5).

Fig. 3(a and b) shows that particles in the neck region have significantly faster relaxation than those in the bulk. This is characterized in two ways: Fig. 3(a) fixes γ at a value where the neck has already formed, showing that C_B decays faster for the neck and slower for the bulk, as a function of Δγ. On the other hand, Fig. 3(b) fixes the strain lag Δγ = 0.05 and varies γ, showing how the difference in relaxation rate grows, as γ increases and the neck develops. This is a microscopic signature of the prediction of the continuum theory, that plastic flow events occur preferentially in the neck and is an indication the local stress can be used as the varying internal field variable χ proposed in the continuum theory.

Fig. 3 (a) Dynamical correlation function C_B(γ,Δγ) on varying Δγ, comparing the neck (Z = 0.5) and the bulk (Z = 1). The correlation function decays as the particles move away from their neighbours. (b) Dynamical correlation on varying γ at fixed lag Δγ. Increasing separation of the curves indicates dynamical contrast between the neck and the bulk. (c) The ensemble-averaged number of neighbours n_n. (d)–(f) The ensemble-averaged number of tetrahedra n_tet, triangular bipyramids n_tb and pentagonal bipyramids n_pb in which a particle participates. The legend in (a) is common to all panels. Data was averaged over 10 trajectories, parameters are the same as shown in Fig. 1. Shaded regions indicate the standard deviation. The dotted line indicates the onset of necking (the strain beyond which N_m has decreased significantly).

We also analyze the local structure of the strand using the topological cluster classification (TCC).⁷¹ This provides a detailed characterisation of local packing, by identifying specific geometrical structures within the system. We write n_n(Z;γ) for the average number of neighbouring particles around a particle, at position Z and strain γ. Similarly we write n_tet for the average number of fully-bonded tetrahedra in which a particle participates (the total number of structures divided by the number of colloids in a Z-segment), and n_tb and n_pb for numbers of triangular and pentagonal bipyramids,⁷² normalised in the same way. (These local packing motifs are illustrated in Fig. 3.)

Fig. 3(c–f) shows these quantities for both the neck and the bulk, as γ increases. The average number of neighbors n_n and the average number of tetrahedra n_tet in which a particle participates remain almost constant with γ, until the system starts to break around γ = 0.36. (The large error estimates in this region arise because breakage occurs at different values of γ for different runs.) Both these quantities follow the number of neighbors, which means that the neck differs significantly from the bulk only when the strand is close to breaking.

However, we find different behavior for the average numbers of triangular bipyramids n_tb and pentagonal bipyramids n_pb in which a particle participates. These measurements correspond to larger structures of 5 and 7 particles respectively, which are sensitive to details of the packing, and tend to be more common in amorphous materials that are deep in the energy landscape.^73,74 Hence n_tb and n_pb are larger in materials that are well-annealed and stable. The neck shows a deficit in these low-energy structures, compared with the bulk (where the number is almost constant): these differences appear as the neck begins to form (γ ≈ 0.1). This effect can be understood in terms of the faster dynamics in the neck, which tends to break up low-energy local structures, as happens in sheared bulk samples.^{20,30,37,73–75} This occurs as particles keep flowing away from the neck, and hence the blue curves in Fig. 3(e and f) follow the neck thinning shown Fig. 1(d).

2.4 Local stress

In order to bridge scales between the particle-level and the whole strand, we measure the local stress, which is the fundamental object for controlling rheology.

The stress [small sigma, Greek, macron]^zz(Z) corresponds to the tensile force in the strand divided by its cross-sectional area, at position Z. This quantity was already discussed in Fig. 1(f), which shows its behaviour after averaging over Z. The Z-dependence of this quantity is shown in Fig. 4(a and b) and Movie S2 (ESI†): one sees that the formation of the neck is accompanied by a local increase in stress. This result was already anticipated in eqn (2): the tension is constant along the strand, so the stress must be larger in segments where the cross-sectional area is smaller. Fig. 4(c) shows the development of the neck, and Fig. 4(d) compares the stress in the neck with the stress in “the bulk” (far from the neck). The central observation is that the stress in the neck reaches a plateau at the yield stress of the strand, so plastic events continue to occur there, allowing the strand to elongate. As the neck develops, the reduction in area at (almost) constant stress explains the reduction in the tensile force [Fig. 1(e)]. This decreasing force leads to a reduction in the bulk stress, which falls below σ_Y. Hence these measurements show directly the mechanism by which plastic events concentrate in the neck.

Fig. 4 (a) The Z-dependent stress [small sigma, Greek, macron]^zz at initialization. (b) The stress at γ = 0.4. All particles in a discretized Z-segment are coloured to the corresponding [small sigma, Greek, macron]^zz. (c) The number of particles in different segments of the strand, as the strain increases, for the single trajectory shown in (a) and (b). (d) The average stress profile [small sigma, Greek, macron]^zz(Z) at the neck (Z = 0.5) and in bulk (Z = 1) as a function of strain γ. (Data averaged over 6 runs.) Colored area is the standard deviation of different runs. (e) and (f) The local stress σ^zz, for the same configurations shown in (a) and (b). (g) The anisotropic stress measurement S̄(Z) as a function of strain γ. (h) The stress distribution curves just before breakage (γ = 0.35–0.4) for the neck and the bulk, compared with the unstrained distribution. Measurement is done over a single run.

A more detailed analysis of the stress reveals additional complexity. Fig. 4(e) and (f) and Movie S3 (ESI†) show the zz component of the local stress σ^zz, averaged over small regions of size (1.5l)³. Before elongation, one might have imagined that the local stress in the strand [Fig. 4(e)] would be close to the average behaviour [Fig. 4(a)]. Instead, σ^zz has strong inhomogeneities: its average is positive (corresponding to tension), but there are significant regions where σ^zz < 0. These inhomogeneities are residual stresses, arising from the (non-equilibrium) process of initialising the strand. We verified that the resulting spatial structures are long-lived and tied to the amorphous structure of the material (see Appendix A), they are not the result of fast thermal fluctuations.

We also checked that the strands satisfy local force-balance (div σ = 0), the residual stresses are divergence-free and arise from the inhomogeneous (amorphous) microstructure of the strands. Since σ^zz depends quite strongly on z, this requires that the stress also has significant off-diagonal components. Such residual stresses were not included in the simple rheological models considered above, where the stress W depended only on Z. However, bridging between rheological models and the local motion of individual particles requires some analysis of these stresses. For example, the plasticity χ presumably depends on the local structure of the strand and the residual stresses there – this offers the opportunity for a deeper understanding of the physical meaning of this field, if it could be connected to particle-level observables.

To explore this, Fig. 4(f) shows the behavior as the strand is elongated and the neck forms. An interesting feature is that the local stress in the neck region appears more homogeneous, compared to the bulk. To quantify this observation, we measure the anisotropy of the stress tensor σ: we define S̄(Z) by averaging the (square root of the) second stress invariant over a given segment of the strand (see Section 4.4). This quantity tends to be large when the stress has large fluctuations away from its average, within a segment; it is smaller when the stress is homogeneous. After averaging over many trajectories, Fig. 4(g) compares S̄(Z) in the neck with the bulk of the strand. It shows that the neck does indeed have a more homogeneous stress field.

Also, Fig. 4(h) shows the distribution of σ^zz, within the neck and the bulk. One sees that the probability to see σ^zz < 0 is suppressed in the neck. In fact, it seems that the reduction in this probability accounts for most of the increase in the average [small sigma, Greek, macron]^zz in the neck (the typical value of σ^zz barely increases there). In comparison with these results, we also simulated the elongation and rupture of crystalline strands (details are given in Appendix B). In this case, the residual stresses are (mostly) absent and there is no homogenization of the stress within the neck [for example, S̄(Z) does not decrease there].

Coupled with Fig. 3, these results show that plastic flow in the neck is coupled with two distinct local changes: reduction in low-energy (stable) structures, but also a homogenization of the residual stresses [smaller S(Z)]. Physically, we can think of this process as a partial fluidisation of the neck, where it loses its amorphous-solid character and starts to move as a highly-viscous fluid, under the applied stress. The solid is characterized by stable local structure and large residual stresses; the fluid is more disordered in its structure but more homogeneous in stress. Similar behavior has been observed in metallic glasses in the spreading of shear transformation zones.^76,77 Overall, these results illustrate the subtle interplay of structure, forces, and dynamics, which is required to characterize the multi-scale phenomenon of yielding and fracture in these amorphous systems.

2.5 Robustness of results for different model parameters

The results presented so far focused on representative parameter values: attraction strength , strain rate [small gamma, Greek, dot above] = 10⁻⁶ and strand thickness r₀ = 4l. However, the general picture that we present is robust. This is shown in Fig. 5 where we vary the strain rate [small gamma, Greek, dot above], the strand thickness r₀, and the cohesive energy .

Fig. 5 Dependence on model parameters. (a), (c) and (e) show the average neck thickness and (b), (d) and (f) show the tensile force in the strand. In (a) and (b) the attraction strength is varied; in (c) and (d) it is the elongation rate [small gamma, Greek, dot above]; and in (e) and (f) it is the strand radius r₀. All lines are averaged over 6 runs with identical starting conditions for a given set of parameters. The blue lines correspond to baseline parameters: , r₀ = 4l and [small gamma, Greek, dot above] = 5 × 10⁻⁶.

Varying the strength of attractive forces [Fig. 5(a and b)], there is weak dependence of the neck thickness N_m on . The tensile force in the strand is larger when the attractive forces are stronger (as expected), but the qualitative behavior of the stress remains the same. (This is also true for the local stress.) We did find that for larger , the position of the neck formation was more predictable: the distribution analogous to Fig. 1(g) was more sharply-peaked, data not shown. This may be expected since thermal fluctuations become less relevant as one approaches the athermal limit .

The effect of varying elongation rate is shown in Fig. 5(c and d). It is notable that the neck thickness N_m is almost unchanged as [small gamma, Greek, dot above] varies over two orders of magnitude. For yet larger strain rates, the thinning process of the neck is accelerated but the qualitative picture remains the same. Fig. 1(f) already showed that the stress during neck formation depends weakly on the strain rate, although faster elongation does create a stress overshoot. This picture is confirmed by Fig. 5(d) which shows the corresponding tensile force (which is simulation box-size independent as can be seen in the Appendix A.2). We also performed simulations for even faster elongation: the picture of a gradual (viscosity-driven) necking instability remains robust and we did not observe any sudden (elasticity-driven) failure mode,^44,45,78 nor solid-like brittle failure. An interesting feature of these systems is that higher elongation rates cause the strand to break at lower total strain [Fig. 5(d)], in contrast to some studies which used athermal systems.^78,79 The reason is that for slow elongation rates, thermally-activated motion allows the neck to be drawn out further before breakage. For faster elongation, the rates of thermally-activated processes drop below [small gamma, Greek, dot above] (an estimate of this is given in Section 4.5), the neck is less fluid, and the strand breaks sooner.

Finally, varying the strand thickness r₀ [Fig. 5(e and f)] one sees a general trend that thicker strands fail at larger strains (and support somewhat larger stresses), but the qualitative behaviour is again the same.

3 Discussion

We have addressed some of the challenges posed by yielding of gels of sticky spheres, by addressing the elementary unit of failure, which is the breaking of a single strand of the gel network. Such local events are important in gel rheology as they trigger collective reorganization of the macroscopic network.³² We implemented a continuum description, and we also analyze our simulations at a single-particle level using higher-order structure, and also at a more coarse-grained level with the stress tensor.

At the level of the whole strand, we find that failure occurs by a necking instability, driven by plastic flow. This can be captured by a simple continuum model, following.^45,46 This description of necking will be a useful basis for future models of failure in these challenging heterogeneous materials.

At the microscopic level, we showed that particle dynamics are faster in the neck, and that complex higher-order structures associated with rigidity are significantly suppressed prior to failure. Remarkably, the structural signature of higher-order clusters, like triangular bipyramids n_tb and pentagonal bipyramids n_pb, appears to be stronger in this system than a previous study of glasses under shear.^74,75,80 As these structures are associated with rigidity,⁷⁴ this gives another indication that the gel fluidizes significantly in the neck.

To bridge between these levels, we analysed the local stress: this is a fundamental quantity for the rheology, which can also be analysed at the microscopic level. Our results demonstrate links between the local structure and the local stress: the structural and dynamical changes in the neck are accompanied by a reduction in the residual stresses, which makes the stress field more homogeneous.

Our microscopic resolution of the stress provides mechanistic insight, in that the stress [small sigma, Greek, macron]^zz(Z) has a value that remains close to the yield stress σ_Y, as the neck develops [Fig. 4(d)]. The constant tensile force in the strand means that other (thicker) parts of the strand have [small sigma, Greek, macron]^zz(Z) < σ_Y, so that plastic events are increasingly concentrated in the neck. This positive feedback drives the (plastic) necking instability, and the mechanism is consistent with the simple continuum model.

On the other hand, the simulations reveal large residual stresses in the strand, which are not accounted for directly in the continuum approach. These results reinforce the observation that the amorphous structure of the strand is heterogeneous, both in its local structure and in the stress field. Consequences of this heterogeneity include Fig. 1(g), which shows the place where the strand is most likely to break. This issue – of predicting rupture – would seem to be vital for predicting and controlling properties of gels. It is also notable that particles have significant mobility within the strand as plastic motion takes place, recall Fig. 2. Such issues are not easily addressed by continuum modelling, their resolution will require input from particle-level data. For example, a simple von Mises criterion for yielding in soft solids^81,82 predicts that a sample will fail in places where the local stress is large and anisotropic: the residual stresses mean that such locations do exist in our samples, but we do not find them to be correlated with the failure.

An alternative approach – building on the continuum theories – would be to search for links between the continuum-level plasticity field χ and the stress anisotropy (or local structure) at the microscopic level. It is certainly true that some aspects of the microscopic structure should influence the continuum models through such a field, as well as through sample-to-sample variations in constitutive parameters such as the yield stress and elastic modulus. These issues could be usefully investigated in future work.

In addition to these questions about the failure of a single strand, our results can also inform future studies of gel behaviour. In particular, we imagine extending these measurements to gel samples under shear,²⁹ or to analyse the breakage of strands during coarsening.⁴³ In these cases, the gel will have many strands, and its overall failure will depend on which ones break first, and on how aging affects the tendency to failure. We also note that the initialization procedure for our strands does not capture the full process by which structures form in realistic gels. Instead, our model strands provide a controlled baseline for measuring behavior under elongation. It would be interesting to investigate the similarities and differences between realistic gels and these model strands, and to understand in more detail the dependence of strands' failure on their history, including effects of aging.

Nevertheless, despite the differences between single strands and macroscopic gel networks, the results here indicate that the more fluid behaviour of the strand in the neck region might be a useful early-warning signal of strand breakage. Such signals would be helpful for characterising and predicting the behaviour of gels as they undergo ageing and/or collapse. An important extension of this work would be a detailed comparison of our results with experiments, for example rheological studies⁸³ and confocal microscopy observations.⁶⁵ Finally, we note recent work where imaging has enabled forces between colloids to be inferred, which could provide a direct comparison of the stress field at the particle level.⁸⁴

4 Methods

4.1 Simulation model and dimensionless parameters

We consider a bidisperse system of N colloidal particles, the two species have diameters 1.04l and 0.96l, where l is the average size. They interact by a Morse potential, which mimics the short-ranged depletion interaction^60,61where r_ij = r_j − r_i, r_ij = |r_ij|, is the interaction strength, is the mean of the diameters of particles i and j, and α sets the attraction range (we take α = 25l⁻¹). While this is longer-ranged than the “Baxter limit” of vanishing range,² different ranges can be mapped by a law of corresponding states⁶² and the dynamics can also be scaled.⁸⁵

The equation of motion is (1), in which the solvent force is where T is the temperature, and ξ a unit white noise. The particles move in a periodic simulation box of dimensions L_⊥ × L_⊥ × L_‖ with L_⊥ = 20l and L_‖(t = 0) = 30l. The Brownian time is τ_b = λ₀l²/(24k_BT). Two natural dimensionless parameters of the model are which measures the strength of the depletion attraction, and which measures the solvent damping. In addition, the system undergoes elongation at shear rate [small gamma, Greek, dot above]₀. This yields an additional dimensionless parameter [small gamma, Greek, dot above] = τ_b[small gamma, Greek, dot above]₀.

In practice, the strain is increased stepwise, with a fixed time period of Δt = 250τ_b between steps, so the strain in each step is [small gamma, Greek, dot above]₀Δt. When measuring the stress (see below), we average the results over the period Δt, to reduce the effects of (fast) thermal fluctuations. To confirm that the stepwise elongation does not affect the results, we also performed simulations with continuous elongation, which results in very similar behaviour.

The numerical simulations are implemented in the LAMMPS package⁸⁶ in which the natural time scale is , this means that τ_b/τ₀ = λ/24 ≈ 0.42 for the parameters used here. The integration time step is 10⁻⁴τ₀. The natural unit of both pressure and stress is k_BT/l³; all results for such quantities are quoted relative to this baseline.

To compare to relevant experimental units, we realize that corresponding experimental systems of sticky colloidal gels have typical diameters [small script l] between 0.1 and 3 μm.^61,87 This allows for a variation of corresponding Brownian times between τ_b ∼ 0.005–7 s and pressure units between k_BT/σ³ ∼ 10⁻⁴–10¹ Pa.

4.2 Preparation of the model gel strand

To set up the strand before elongation, we first initialize a bulk simulation of the model colloid in the NPT ensemble with a low interaction strength and a constant pressure. We set P₀ = 0.16 and slowly increase the attractive strength to , using steps of every 2τ_b. This causes the volume fraction to increase to ϕ ≈ 0.59, the system remains homogeneous because this isobaric transformation not enter the spinodal decomposition regime.⁸⁷ We allow this dense homogeneous fluid to relax for a time τ_glass.

We then switch to the NVT ensemble and instantaneously adjust to the desired value, in the range 4.5–10. The strong attractive interactions induce additional dynamical arrest and the system forms a glass-like system. After this point, we allowed the simulations to relax for another time τ_glass. The system is now inside the liquid–vapour binodal, but the glassy dynamics are slow enough that phase separation is not observed. (We performed simulations with τ_glass up to 10⁵τ_b, the results depend very weakly on this parameter. Results are shown for τ_glass = 10⁴τ_b.)

The result is a homogeneous glassy state with volume fraction ϕ ≈ 0.59. Our initial strand is obtained by excising a cylinder of radius r₀ from this system, after which we run dynamics for approximately 1000τ_b, to allow the system to relax any features that are artefacts of cutting out the cylinder. We then start the elongation process.

A different initialisation technique was developed in ref. 88, which might lead to better equilibration of the surface structure of the strand. However, we do not expect strands of particulate gels to be particularly well-equilibrated, so our simple initialization method is appropriate here.

4.3 Stress measurement

Our measurements of local stress use a volume-averaged representation of the Irving–Kirkwood stress.^64,65,89 Write p_i for the momentum of particle i, also r_ij for the vector connecting particles i and j, and f_ij the corresponding interparticle force. Now consider a spatial region Ω whose volume is |Ω|. The μν component of the IK stress for that region iswhere ϑ_i,Ω = 1 if particle i is in Ω and zero otherwise; similarly φ_ij,Ω is the fraction of the straight line connecting particles i,j that lies within Ω.

Taking Ω to be the entire simulation box Ω* gives the total stress Σ, which can also be computed from the virial. The tensile force is then .

For a local measurement of stress at point r, we take Ω to be a small cube of side l_IK, centred at r. The resulting stress is denoted by σ(r). For local measurements, we take l_IK = 1.5l, which is sufficiently small to allow a local measurement, but sufficiently large to avoid numerical uncertainties due to thermal fluctuations (see Appendix A.2).

As discussed in the main text, it is sometimes convenient to divide the system into n_seg = 20 segments along the z-direction, each of which has volume V_z = L_⊥²L_‖/n_seg. Then define N_Z as the number of particles in segment Z. Taking Ω in (6) to be one of these segments gives the tensile force in the strand, divided by the cross-sectional area of the simulation box (which is L_⊥²). However, the physically-relevant stress is the tensile force divided by the cross-sectional area of the strand. This is obtained by rescaling (6):

where Ω_Z is the segment of the box at position Z and V_s = N_Zπl³/(6ϕ) is the estimated volume occupied by the strand, within that segment. (Here ϕ = 0.59 is the volume fraction within the strand, so πl³/(6ϕ) is the mean volume per particle there.)

Note that the local IK stress is derived directly from the equations for momentum conservation. As such, it accurately reflects the fact that a locally stable gel strand (whose structure is not changing with time), satisfies local force balance div σ(r) = 0 at every point r, that is,

which holds for ν = x, y, z. See Appendix A.1 for further details. It is not possible to build a local virial stress with this property. (The IK method is not the only way to obtain such a stress tensor, but it is a convenient one.⁶⁵) The method-of-planes method⁶⁵ also measures the stress locally on a given plane, but the volume-averaged method used here is more effective at accurately reflecting the force balance condition (8). For example, if the derivatives in (8) are estimated as finite differences between adjacent measurement volumes then the volume-averaged method stress gives much better numerical agreement with (8) than the method of planes.

4.4 Stress anisotropy

The second invariant of the stress tensor measures the anisotropy of the stress, as a scalar quantity that is independent of the orientation of the coordinate system:This quantity is zero if σ is proportional to the identity, as would be expected in the bulk of a simple fluid. In a region Ω with large anisotropic stresses then J₂ will be large. The von Mises criterion^81,82 for failure of solid materials states that breakage will occur when the local J₂ exceeds a threshold.

For strands under tension with homogeneous stress, the dominant element of σ_Ω is σ^zz ≈ [small sigma, Greek, macron]^zz(Z), leading to J₂ ∝ [small sigma, Greek, macron]^zz(Z)². However, in an strand like the one in Fig. 4(e) with large residual stresses, typical elements of σ_Ω have absolute values larger than [small sigma, Greek, macron]^zz(Z), leading to a much larger value of J_2,Ω. It is convenient to average this quantity over a segment of the strand, as

where the sum runs over cubic regions of size l_IK³, within segment Z. Since J_2,Ω is a (non-negative) measure of anisotropy, one sees that S̄ captures anisotropic stress fluctuations within the strand. (Such fluctuations are averaged away in the cross-sectional stress [small sigma, Greek, macron]^zz.)

4.5 Bond-breaking correlation function

To measure local particle rearrangements, we define b_ij(Δ;γ) = 1 if particles ij are within a distance Δ of each other, when the accumulated shear strain is γ. Then the fraction of neighbours of particle i that are lost between strains γ and γ + Δγ is⁹⁰We take Δ = 1.4l throughout as that is the average cut-off range of particle interaction (other values would have given qualitatively the same results). A correlation function C_B(γ,Δγ) is then obtained by averaging c_i over all particles in a suitable region, which we take here to be the segment of the system with position Z.

As additional context, Fig. 6 shows the bond-breaking correlation function where we separate particles on the surface of the strand from those in the interior. For this, we define the 20 particles in every Z segment at every time that are closest and furthest away from the centre of the strand. Particles on the surface relax faster than those in the interior, as there are fewer bonds at the surface.^54,91 Interestingly, the surface particle relaxation time is independent on the Z-segment position. In contrast, it is clear that the interior particles relax faster in the neck, consistent with their local fluidization.

Fig. 6 Repeat of Fig. 3a, but where we distinguish interior and surface particles. These are defined as the 20 particles closest and furthest from the centre of the strand for every Z-segment.

To estimate the rate of thermally-activated surface processes, we observe from Fig. 6 that the time needed for particles to lose about half their neighbours is about 2 × 10³τ_b (or γ = 0.01 with the applied strain rate [small gamma, Greek, dot above] = 5 × 10⁻⁶). Thus we estimate our thermally-activated surface processes rate as 0.5 × 10⁻³τ_b⁻¹, comparable to the strain rate of the black line in Fig. 5.

4.6 Topological cluster classification

We analysed particles' local environments using the TCC. For each particle i, this yields: (i) its number of neighbours n_n; (ii) the number of fully-bonded tetrahedra in which it participates n_tet; (iii) the numbers of trigonal and pentagonal bipyramids in which it participates, n_tb and n_pb respectively. Details and parameters of the TCC are the same as.⁷¹ These quantities were then averaged over particles in a suitable region (typically the segment of the system with position Z).

Author contributions

The project was devised by RLJ and CPR. Simulations and analysis of numerical results were performed by KT, who also developed the analytic theory, based on discussions with TBL and RLJ. The manuscript was written by KT and RLJ, with contributions from TBL and CPR.

Data availability

The datasets generated during and/or analysed during the current study and code to reproduce the results in this manuscript are available at https://doi.org/10.17863/CAM.101482.

Conflicts of interest

The authors declare no competing interests.

Appendices

These appendices contain additional results and analysis to further justify the methods and conclusions of the main text.

Movie S1 (ESI†) shows a single trajectory of Fig. 1 in the manuscript where an strand undergoes failure. The 16 particles that will make up the neck before breaking are denoted as red. All other particles are blue.

Movie S2 (ESI†) shows the movie corresponding to Fig. 4(a and b).

Movie S3 (ESI†) shows the movie corresponding to Fig. 4(e and f).

Appendix A discusses our measurements of local stress. Appendix B discusses numerical simulations of elongation of a crystalline strand, for comparison with the amorphous gel strands considered in main text. Appendix C discusses the continuum rheological model of elongation and necking, including the fitting to numerical data.

A Stress measurements

A.1 Validation of the Irving–Kirkwood stress

As discussed in the main text, we use the Irving–Kirkwood (IK) stress σ throughout this manuscript to measure the stress on different length scales. Derivations of the IK stress are given in ref. 47, 64 and 67 and implementation methods are shown in ref. 65. In this manuscript, we use the volume averaged stress from ref. 65. This approach ensures that the stress is self-averaging when measured over large volumes Ω, which helps to reduce statistical uncertainties. In particular taking Ω in (6) as the full simulation box, the IK stress reduces to the standard virial stress, which is

For local stress measurements, it is possible to define a local virial stress by restricting the sums in (1) to particles in a particular region. However, this choice does not ensure that stress gradients cause changes in local momentum, nor that force-balanced systems have div σ(r) = 0. To see this, we compute the (local) tensile force F^z(Z) in the strand by taking Ω in (6) as a segment at position Z. Force balance requires that ∂F^z/∂Z = 0, up to small corrections due to the thermal fluctuations. Fig. 7(a) shows that this requirement is obeyed to high accuracy for the IK stress, but it fails for the virial.

Fig. 7 (a) The tensile force calculated from the Irving–Kirkwood stress and the localized virial stress. Data correspond to a single trajectory of Fig. 1(d) at γ = 0.2, without averaging over the z direction. (b) Time-correlation function of the global and local stress. Data used in (b) was obtained from an unstrained strand at the same standard parameters as (a).

The other (non zz) components of the stress tensor also have persistent non-zero values (these are local residual stresses). To illustrate that these components come from the amorphous structure of the strand (instead of fluctuating rapidly on the time scale of velocity fluctuations), we calculate the time correlation functions of the local and global stress:

andIn these sums, both μ and ν run over the three Cartesian directions x, y, z, but the term μ = ν = z is excluded.

Fig. 7(b) compares the global and local stress correlation functions, during simulations of an strand without any elongation. The local stress correlations decay slowly, showing that the residual stresses are long-lived. We attribute this to slow structural changes in the (non-equilibrium) strand, thousands of Brownian times are required for these changes to become significant. By contrast, all elements of the global stress are very small (except zz). As a result, the dominant contributions to C_Σ are fast thermal fluctuations, so this correlation function decays quickly.

It is well-known that microscopic expressions for the local stress tensor are not unique.⁶⁵ We briefly discuss three aspects of this issue. First, since these systems have significant residual stresses, the value of σ depends on the scale at which it is measured. Our local stress tensor is measured by taking Ω in (6) to be a cubic box of side 3l/2. This length scale is chosen for numerical convenience: taking larger boxes tends to smooth out the stress, but smaller boxes lead to a noisy signal (see Appendix A.2 for further discussion of this point). Second, the terms arising from pairwise forces in the IK stress are evaluated by considering a linear interaction path between the particles. Other paths are possible but the linear path is a simple and convenient choice. For long-ranged forces, the choice of path can significantly affect the stress, but for these short-ranged Morse potentials, such effects are small (as long as a reasonable path is used). Third, one might (in principle) also shift all stress values by a global constant (the reference pressure), here we insist that a very dilute colloidal suspension has a vanishingly small (osmotic) pressure, so this constant is zero.

A.2 Averaging the stress

Throughout this work, we report IK stresses that are averaged over a time period T = 250τ_b. We recall from Fig. 7b that this time scale is small enough that the stress does not relax significantly, but we do find that the averaging process reduces statistical noise. (Specifically, we compute the stress at time points separated by τ_b/2 and we average over 500 such time points to obtain the reported values.) The justification of this averaging relies on of the slow time evolution of the local stress (recall Fig. 7), which is in turn due to the dynamically-arrested (solid) structure of the strand. In the bulk of a simple fluid, this kind of averaging would yield instead an isotropic stress tensor, although the situation would be more complicated in systems with interfaces or applied forces.⁶⁶ As noted above, we also measure a volume-averaged local stress, based on a cubic box of size l_IK = 1.5l.

Fig. 8 illustrates the effects of these averaging procedures, for representative strands. For the stress itself, larger averaging times and volumes lead to smoother signals, as expected. This is shown in Fig. 8(a and c). We also computed the stress anisotropy S̄(Z), as defined in eqn (10). In general, larger averaging times and volumes suppress the effects of anisotropic fluctuations. To faithfully capture these fluctuations (as in Fig. 4), we use intermediate length and time scales for averaging, as we now discuss.

Fig. 8 (a) The average cross-sectional tensile stress [small sigma, Greek, macron]^zz, measured as a volume average over a region of size L_‖² × l_IK. Averaging over larger regions reduces the fluctuations, but the increased stress in the neck is clear in all cases (the neck is at Z ≈ 0.35, in this case). (b) Corresponding measurement of anisotropy of the local stress, as a function of the averaging volume (see text for a discussion). [The coloring of lines is the same as in panel (a).] (c) Tensile stress, estimated by time-averaging over a time T_ave. Increasing T_ave smooths the data, but does not affect the signature of the neck, where [small sigma, Greek, macron]^zz increases (at Z ≈ 0.2 in this case). (d) Corresponding stress anisotropy S̄: time-averaging reduces the effect of fast anisotropic fluctuations: for the larger averaging times, one sees a stable signal that comes from long-lived residual stresses. [The coloring of lines is the same as in panel (c).] Notes: blue lines correspond to the parameters used for all stress measurements outside this figure; consistent with this, we take T_ave = 250τ_b in (a) and (b) and l_IK = 1.5l in (c) and (d). Panels (a) and (b) are taken from a trajectory where the neck appears at Z ≈ 0.35; panels (c) and (d) are from a different trajectory where the neck is at Z ≈ 0.2. Both trajectories have the same parameters, which are those of Fig. 1(c and d).

We first consider effects of volume-averaging, over cubes of volume l_IK³. The IK stress distributes the contribution of each pair of particles along a line connecting them. For very small volumes, this results in many cubes with no contribution to the stress, and others with large anisotropic contributions. The resulting mean stress is independent of l_IK but the anisotropic fluctuations behave as S̄ ∼ l_IK⁻¹ for small l_IK. Fig. 8(b) plots l_IKS̄(Z) for various sizes l_IK. Increasing the box size up to l_IK ≈ 1.5l suppresses a noisy contribution from thermal fluctuations, helping to reveal the reduction in S̄ near the neck. For larger averaging volumes, the anisotropic fluctuations of the stress are significantly reduced: this is because the anisotropic residual stresses illustrated in Fig. 1(e and f) are being averaged away. (Indeed, averaging over the entire cross section will eventually yield the picture of Fig. 1(a and b), where anisotropic fluctuations are much smaller.) Hence we choose l_IK = 1.5l for our measurements, which is large enough to suppress fluctuations from thermal noise, without averaging away the physically-relevant residual stresses.

For time averaging, the picture is simpler. Taking an average over T_ave = 250τ_b effectively suppresses fast thermal fluctuations in the measured stress [Fig. 8(c)]. Since these fluctuations are anisotropic, the time-averaging also suppresses the anisotropy [Fig. 8(d)]. There is a broad range of times around T_ave = 250τ_b where these signals are stable. For much larger times, one loses resolution in time due to local stress relaxation (from Fig. 7, this happens on time scales are ≳1000τ_b). Hence we choose T_ave = 250τ_b for our measurements, as a sensible compromise between time-resolution and noise reduction.

A.3 Additional information about stress distribution

Fig. 4(h) of the main text shows the distribution of the local stress (more precisely, its zz component). We emphasized the differences in this distribution between the neck region and the bulk of the strand. We present here some additional information on these distributions.

At initialization, long-lived residual stresses are quenched into the system: the tails of these distributions are roughly exponential, with different decay rates for positive and negative stress. Fig. 9(a) shows that positive stresses are more common, but the negative stress distribution has the fatter tail. Overall, the total stress (which coincides with the average of this distribution) is lightly positive at initialization, corresponding to a tensile stress.

Fig. 9 Distributions of the local stress σ^zz. (a) Comparison of the distributions for positive and negative stress values, before any elongation has occurred. Data is found by measuring for 10⁶τ_b without elongation. (b) Stress distribution evolution far away from the neck (Z = 1) near the start of the simulation (γ = 0–0.05), just after the system has reached its stress plateau (γ = 0.03–0.05), and just before necking occurs (γ = 0.35–0.4). Data has been averaged over 6 trajectories from Fig. 1 in the main text.

Fig. 9(b) shows the stress in the bulk of the strand (away from the neck), as the strand is stretched. The tensile force increases during elongation corresponding to an increase in the average of this distribution. However, this shift is weak, in comparison with the residual stresses that are already present, so it has a weak effect on the stress distribution.

We also note that the change in stress distribution for the neck is not just an effect of its reduced thickness (which leads to a change in the ratio of surface to bulk). We initialized strands with different radii: they all have similar non-Gaussian distributions of the stress (Fig. 10).

Fig. 10 The stress distribution for different starting radii with no deformation. Data is found by measuring for 10⁶τ_b without straining. Other parameter values are those of Fig. 1 of the main text.

B Elongation of a crystalline strand

As a point of comparison for the strand elongation discussed in the main text, we also performed similar experiments on a crystalline strand. We initialised a face-centred cubic (FCC) crystal from which we excised a cylindrical strand. On elongation, a few plasticity events were observed, after which the strand broke by a sudden mechanism resembling brittle failure. Results are shown in Fig. 11. Looking at the mesoscopic stress, we can see that the magnitude of the local stress σ_zz increases at the failure point (Fig. 11(a)), as expected because the strand is thinnest there.

Fig. 11 The local stress σ^zz (a) before neck formation (γ = 0.015) and (b) after neck formation (γ = 0.038) for a system where particles are initialized on an FCC lattice. (c) The corresponding measurement of anisotropy in the stress tensor S̄. Other parameter values are those of Fig. 1 of the main text.

However, a striking difference between the crystalline and amorphous strands is the absence of residual stresses in the crystalline case. As a result, the stress is relatively homogeneous. On computing the stress anisotropy in the crystal, we find that S̄(Z) is largest in the neck. This stands in contrast to amorphous strand, where the stress was more homogeneous in the neck, and S̄(Z) was smaller there.

Comparing the tensile force between a strand and the crystal, we notice a difference in failure mechanisms (Fig. 12). The crystal strand breaks with a clear brittle-like failure as we have a significant stress drop, while the strand has a ductile failure, even though the particles in both simulations have the same parameters. The only difference is the initialization history.

Fig. 12 (a) and (b) The cross-sectional stress [small sigma, Greek, macron]^zz, averaged over z and the tensile force for a single trajectory of the strand of Fig. 1. and (c) and (d) for the crystal. Particle parameters are the same.

C Continuum model

The theory presented in the manuscript is based on the analysis of ref. 46, although it can also be interpreted in the framework of ref. 44 and 45.

C.1 Definition

The model uses a thin-filament approximation so the strand is modelled in one-dimension (oriented along the z direction). We define a local velocity, V(z,t) from which follows a local strain rate, [small gamma, Greek, dot above]_L(z,t) = ∂_zV(z,t), that will differ in general from the externally imposed strain rate [small gamma, Greek, dot above]. We also assume that only the area of the strand A(z,t) is important, when considering mass conservation. We will write our equations in the co-extending frame of imposed strain rate [small gamma, Greek, dot above]. Hence we define⁴⁵

Z = z exp(−[small gamma, Greek, dot above]t), (A4)

v(Z,t) = V(z,t)exp(−[small gamma, Greek, dot above]t), (A5)

a(Z,t) = A(z,t)exp([small gamma, Greek, dot above]t), (A6)

where we have normalized z by the strand length L_⊥(t), to keep the new spatial variable Z between 0 and 1. Transformation to this co-extending frame, the mass conservation equation becomes(Here and throughout we neglect advective terms from the change of frame, which is valid for slow elongation rates.)

To account for forces in the strand, we write W(Z,t) for the tensile stress in the z direction, averaged over the cross-section of the strand. The tensile force is then aW and the system remains force-balanced at all times, so

We assume that the stress evolves with an elastic loading term and relaxes due to local plasticity events which yieldswhere G is the elastic modulus and p the rate of plastic relaxation, which depends in turn on an internal plasticity field χ (see below). More specifically, we follow ref. 46: p is zero when W is less than a yield stress σ_Y, and χ behaves like a temperature, which controls the rate of “activated” plastic events. (As discussed in the main text, such relationships are familiar from shear transformation zone theory.⁶⁸) Hence,where θ is the Heaviside (step) function and τ_p(χ) = τ_ref exp(1/χ), where τ_ref is a reference time scale.

In the simplest case we take χ as a constant parameter, but this is not sufficient to capture the stress–strain relationship observed in simulations. Instead we again follow⁴⁶ in promoting χ to a dynamical variable whose relaxation is controlled by the plastic time scale τ_p. Specifically, we take

where the parameter χ_∞ sets the steady-state value of χ.

C.2 Homogeneous solution: qualitative behaviour and fitting to numerical data

We first consider homogeneous solutions of this model, that is, solutions where strain, area, W and χ are independent of Z. At early times, we have W(t) = W_t=0 + G[small gamma, Greek, dot above]t, which holds for times short enough that W < σ_Y. For long times, we have W(t) → σ_Y/(1 − [small gamma, Greek, dot above]τ_p) where τ_p is the plastic time in the steady state.

As discussed in the main text, an important question is: when W reaches the yield stress σ_Y and the plastic activity starts to occur, does χ relax quickly to χ_∞, or is the time scale for this relaxation compete with the relaxation of W to its steady-state value? This determines whether a stress overshoot is observed.

The resulting model has several parameters, which we fit to the data of Fig. 1(e) in a multi-step procedure. We identify W with [small sigma, Greek, macron]^zz (averaged over Z) and we fit the data at early times to [small sigma, Greek, macron]^zz = W_t=0 + G[small gamma, Greek, dot above]t, which yields G = 420 and W_t=0 = 0.5 [recall that the units of both these quantities are k_BT/l³]. We then consider the plateau of [small sigma, Greek, macron]^zz (before significant necking occurs): we fit the plateau height to W = σ_Y/(1 − [small gamma, Greek, dot above]τ_p) which yields τ_p = 180τ_b and σ_Y = 3.3. Note however that this τ_p is the steady-state value of the plastic time scale: since τ_p = τ_refe^1/χ then this is not itself a parameter of the model. To obtain values for τ_ref and χ_∞ we fit the stress overshoot that occurs in Fig. 1(e) for [small gamma, Greek, dot above] = 5 × 10⁻⁴: this yields τ_ref = 180τ_b, and χ_∞ = 80, as well as the initial condition χ_t=0 = 0.65. At this point all parameters have been determined for the fitting in Fig. 1(e).

C.3 Necking as a linear stability of the homogeneous solution

To analyse necking, we consider the linear stability of the homogeneous solution for ([small gamma, Greek, dot above]_L, a, W, χ). We consider a perturbation to the homgeneous solution at wavevector q, that iswhere the 0-subscripts indicate the solution to the homogeneous equation. Following,⁴⁶eqn (A7) fixes [small gamma, Greek, dot above]_L in terms of a and its time derivative. Linearising eqn (A8) determines δa in terms of δW, a₀, W₀. Linearising (A9), (A11) and using these two relations yields a closed equation for the perturbation aswith stability matrixwhere the notation [small chi, Greek, dot above] is a shorthand for the right-hand side of (A11), interpreted as a function of (W₀, χ₀).

For W₀ < σ_Y, there is no plastic flow and M = 0: the strand responds elastically and preserves its shape under elongation. However, for W₀ > σ_Y, the matrix always has one positive and one negative eigenvalue, indicating that the homogeneous solution is linearly unstable and a neck will form. This situation – of a single real positive eigenvalue – corresponds to the slow (or gradual) instability of ref. 44 and 45, which is consistent with our numerical simulations.

Acknowledgements

We thank Daan Frenkel, Sylvain Patinet, Camille Scalliet, Amin Doostmohammadi, Abraham Mauleon-Amieva, Rui Cheng, and Malcolm Faers for helpful discussions. This work was supported by the EPSRC through grants EP/T031247/1 (KT and RLJ) and EP/T031077/1 (CPR and TBL). In the later stages of the project, KT also received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska Curie grant agreement no. 101029079.

Notes and references

1. C. P. Royall, P. Charbonneau, M. Dijkstra, J. Russo, F. Smallenburg, T. Speck and C. Valeriani, arXiv, 2023, preprint, arXiv:2305.02452,  DOI:10.48550/arXiv.2305.02452.
2. R. Baxter, J. Chem. Phys., 1968, 49, 2770–2774.
3. K. N. Pham, A. M. Puertas, J. Bergenholtz, S. U. Egelhaaf, A. Moussad, P. N. Pusey, A. B. Schofield, M. E. Cates, M. Fuchs and W. C. Poon, Science, 2002, 296, 104–106.
4. E. Zaccarelli, J. Phys.: Condens. Matter, 2007, 19, 323101.
5. D. Bouttes, E. Gouillart, E. Boller, D. Dalmas and D. Vandembroucq, Phys. Rev. Lett., 2014, 112, 245701.
6. R. E. Baumer and M. J. Demkowicz, Phys. Rev. Lett., 2013, 110, 145502.
7. C. P. Royall and S. R. Williams, J. Phys. Chem. B, 2011, 115, 7288–7293.
8. S. Jabbari-Farouji, G. H. Wegdam and D. Bonn, Phys. Rev. Lett., 2007, 99, 065701.
9. J. Li, Y. Cao, B. Kou, X. Xiao, K. Fezzaa and Y. Wang, Nat. Commun., 2014, 5, 5014.
10. J. J. McManus, P. Charbonneau, E. Zaccarelli and N. Asherie, Curr. Opin. Colloid Interface Sci., 2016, 22, 73–79.
11. H. Tanaka, Faraday Discuss., 2013, 167, 9–76.
12. M. E. Helgeson, S. Moran, H. An and P. S. Doyle, Nat. Mater., 2012, 11, 344–352.
13. J. L. Drury and D. J. Mooney, Biomaterials, 2003, 24, 4337–4351.
14. C. P. Royall, M. A. Faers, S. L. Fussell and J. E. Hallett, J. Phys.: Condens. Matter, 2021, 33, 453002.
15. A. Dinsmore and D. Weitz, J. Phys.: Condens. Matter, 2002, 14, 7581.
16. P. Bartlett, L. J. Teece and M. A. Faers, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 021404.
17. W. Poon, J. Phys.: Condens. Matter, 2002, 14, R859.
18. P. J. Lu, E. Zaccarelli, F. Ciulla, A. B. Schofield, F. Sciortino and D. A. Weitz, Nature, 2008, 453, 499–503.
19. H. Tsurusawa, M. Leocmach, J. Russo and H. Tanaka, Sci. Adv., 2019, 5, eaav6090.
20. A. Boromand, S. Jamali and J. M. Maia, Soft Matter, 2017, 13, 458–473.
21. C. Patrick Royall, S. R. Williams, T. Ohtsuka and H. Tanaka, Nat. Mater., 2008, 7, 556–561.
22. V. Trappe, V. Prasad, L. Cipelletti, P. Segre and D. A. Weitz, Nature, 2001, 411, 772–775.
23. M. Bantawa, B. Keshavarz, M. Geri, M. Bouzid, T. Divoux, G. H. McKinley and E. Del Gado, Nat. Phys., 2023, 19, 1178–1184.
24. L. Cipelletti and L. Ramos, J. Phys.: Condens. Matter, 2005, 17, R253.
25. K. Masschaele, J. Fransaer and J. Vermant, J. Rheol., 2009, 53, 1437–1460.
26. J. Sprakel, S. B. Lindström, T. E. Kodger and D. A. Weitz, Phys. Rev. Lett., 2011, 106, 248303.
27. V. Grenard, T. Divoux, N. Taberlet and S. Manneville, Soft Matter, 2014, 10, 1555–1571.
28. B. J. Landrum, W. B. Russel and R. N. Zia, J. Rheol., 2016, 60, 783–807.
29. L. C. Johnson, B. J. Landrum and R. N. Zia, Soft Matter, 2018, 14, 5048–5068.
30. N. Koumakis, E. Moghimi, R. Besseling, W. C. Poon, J. F. Brady and G. Petekidis, Soft Matter, 2015, 11, 4640–4648.
31. A. Nicolas, E. E. Ferrero, K. Martens and J.-L. Barrat, Rev. Mod. Phys., 2018, 90, 045006.
32. M. Bouzid, J. Colombo, L. V. Barbosa and E. Del Gado, Nat. Commun., 2017, 8, 15846.
33. S. M. Fielding, P. Sollich and M. E. Cates, J. Rheol., 2000, 44, 323–369.
34. M. Nabizadeh and S. Jamali, Nat. Commun., 2021, 12, 4274.
35. S. Patinet, D. Vandembroucq and M. L. Falk, Phys. Rev. Lett., 2016, 117, 045501.
36. J. Pollard and S. M. Fielding, Phys. Rev. Res., 2022, 4, 043037.
37. D. Richard, M. Ozawa, S. Patinet, E. Stanifer, B. Shang, S. A. Ridout, B. Xu, G. Zhang, P. K. Morse, J.-L. Barrat, L. Berthier, M. L. Falk, P. Guan, A. J. Liu, K. Martens, S. Sastry, D. Vandembroucq, E. Lerner and M. L. Manning, Phys. Rev. Mater., 2020, 4, 113609.
38. J.-L. Barrat, Phys. A, 2018, 504, 20–30.
39. S. B. Lindström, T. E. Kodger, J. Sprakel and D. A. Weitz, Soft Matter, 2012, 8, 3657–3664.
40. J. Colombo and E. Del Gado, J. Rheol., 2014, 58, 1089–1116.
41. J. E. Verweij, F. A. Leermakers, J. Sprakel and J. Van Der Gucht, Soft Matter, 2019, 15, 6447–6454.
42. J. Song, Q. Zhang, F. de Quesada, M. H. Rizvi, J. B. Tracy, J. Ilavsky, S. Narayanan, E. Del Gado, R. L. Leheny and N. Holten-Andersen, et al. , Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2201566119.
43. V. Testard, L. Berthier and W. Kob, J. Chem. Phys., 2014, 140, 164502.
44. D. M. Hoyle and S. M. Fielding, Phys. Rev. Lett., 2015, 114, 158301.
45. D. M. Hoyle and S. M. Fielding, J. Rheol., 2016, 60, 1347–1375.
46. A. Moriel and E. Bouchbinder, Phys. Rev. Mater., 2018, 2, 073602.
47. J. Irving and J. G. Kirkwood, J. Chem. Phys., 1950, 18, 817–829.
48. M. Tsamados, A. Tanguy, F. Léonforte and J. L. Barrat, Eur. Phys. J. E: Soft Matter Biol. Phys., 2008, 26, 283–293.
49. M. Tsamados, A. Tanguy, C. Goldenberg and J.-L. Barrat, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 026112.
50. H. Vinutha, F. D. Diaz Ruiz, X. Mao, B. Chakraborty and E. Del Gado, J. Chem. Phys., 2023, 158, 114104.
51. S. Zhang, E. Stanifer, V. V. Vasisht, L. Zhang, E. Del Gado and X. Mao, Phys. Rev. Res., 2022, 4, 043181.
52. K. Yoshimoto, T. S. Jain, K. Van Workum, P. F. Nealey and J. J. de Pablo, Phys. Rev. Lett., 2004, 93, 175501.
53. H. Mizuno, S. Mossa and J.-L. Barrat, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 042306.
54. J. M. van Doorn, J. E. Verweij, J. Sprakel and J. van der Gucht, Phys. Rev. Lett., 2018, 120, 208005.
55. D. Richard, E. T. Lund, J. Schroers and E. Bouchbinder, Phys. Rev. Mater., 2023, 7, L032601.
56. E. Y. Lin and R. A. Riggleman, Soft Matter, 2019, 15, 6589–6595.
57. G. Brambilla, S. Buzzaccaro, R. Piazza, L. Berthier and L. Cipelletti, Phys. Rev. Lett., 2011, 106, 118302.
58. E. Secchi, S. Buzzaccaro and R. Piazza, Soft Matter, 2014, 10, 5296–5310.
59. T. Gibaud, T. Divoux and S. Manneville, Statistical and Nonlinear Physics, Springer, 2022, pp. 313–336.
60. J. Taffs, A. Malins, S. R. Williams and C. P. Royall, J. Phys.: Condens. Matter, 2010, 22, 104119.
61. A. Razali, C. J. Fullerton, F. Turci, J. E. Hallett, R. L. Jack and C. P. Royall, Soft Matter, 2017, 13, 3230–3239.
62. M. G. Noro and D. Frenkel, J. Chem. Phys., 2000, 113, 2941–2944.
63. A. Widmer-Cooper and P. Harrowell, J. Chem. Phys., 2007, 126, 154503.
64. J. Z. Yang, X. Wu and X. Li, J. Chem. Phys., 2012, 137, 134104.
65. E. Smith, D. Heyes and D. Dini, J. Chem. Phys., 2017, 146, 224109.
66. C. Braga, E. R. Smith, A. Nold, D. N. Sibley and S. Kalliadasis, J. Chem. Phys., 2018, 149, 044705.
67. J. Dolezal and R. L. Jack, Phys. Rev. Res., 2022, 4, 033134.
68. M. L. Falk and J. S. Langer, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57, 7192–7205.
69. A. R. Hinkle, C. H. Rycroft, M. D. Shields and M. L. Falk, Phys. Rev. E, 2017, 95, 053001.
70. K. Kontolati, D. Alix-Williams, N. M. Boffi, M. L. Falk, C. H. Rycroft and M. D. Shields, Acta Mater., 2021, 215, 117008.
71. A. Malins, S. R. Williams, J. Eggers and C. P. Royall, J. Chem. Phys., 2013, 139, 234506.
72. R. L. Jack, A. J. Dunleavy and C. P. Royall, Phys. Rev. Lett., 2014, 113, 095703.
73. R. Pinney, T. B. Liverpool and C. P. Royall, J. Chem. Phys., 2016, 145, 234501.
74. R. Pinney, T. B. Liverpool and C. P. Royall, Phys. Rev. E, 2018, 97, 032609.
75. J. Ding, S. Patinet, M. L. Falk, Y. Cheng and E. Ma, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 14052–14056.
76. D. Sopu, A. Stukowski, M. Stoica and S. Scudino, Phys. Rev. Lett., 2017, 119, 195503.
77. X. Bian, D. S opu, G. Wang, B. Sun, J. Bednarcik, C. Gammer, Q. Zhai and J. Eckert, NPG Asia Mater., 2020, 12, 59.
78. L. O. Eastgate, J. S. Langer and L. Pechenik, Phys. Rev. Lett., 2003, 90, 045506.
79. C.-C. Kuo and M. Dennin, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 052308.
80. F. Turci, T. Speck and C. P. Royall, Eur. Phys. J. E: Soft Matter Biol. Phys., 2018, 41, 54.
81. S. Timoshenko, Strength of Materials: Part 2, Advanced Theory and Problems, van Nostrand, 1956.
82. L. U. Sica, P. R. de Souza Mendes and R. L. Thompson, Soft Matter, 2020, 16, 7576–7584.
83. M. Laurati, S. Egelhaaf and G. Petekidis, J. Rheol., 2011, 55, 673–706.
84. J. Dong, F. Turci, R. L. Jack, M. A. Faers and C. P. Royall, J. Chem. Phys., 2022, 156, 214907.
85. G. Foffi, C. De Michele, F. Sciortino and P. Tartaglia, Phys. Rev. Lett., 2005, 94, 078301.
86. LAMMPS, https://www.lammps.org/.
87. C. P. Royall, S. R. Williams and H. Tanaka, J. Chem. Phys., 2018, 148, 044501.
88. Y. Shi, Appl. Phys. Lett., 2010, 96, 121909.
89. R. J. Hardy, J. Chem. Phys., 1982, 76, 622–628.
90. C. Scalliet, B. Guiselin and L. Berthier, Phys. Rev. X, 2022, 12, 041028.
91. J. M. Van Doorn, J. Bronkhorst, R. Higler, T. Van De Laar and J. Sprakel, Phys. Rev. Lett., 2017, 118, 188001.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3sm00681f

---