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

‘‘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 ineﬀective. 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 diﬀerent 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 diﬀerent 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.


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. 1 Although their equilibrium phase behavior has been, to a considerable extent, understood since the 1960s, 2 out of equilibrium, sticky spheres exhibit a bewildering range of phenomena, including multiple glassy states 3 and gelation, which is a complex blend of phase separation and dynamical arrest leading to a far-from-equilibrium unstable state. 4ch 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 5 and metallic glassformers, 6 small molecules, 7 clays, 8 and granular matter. 9Among the most wellknown gel-forming systems are soft materials, including proteins, 10 foods, 11 hydrogels, 12 tissues 13 and perhaps most studied of all, colloids. 4,14els typically consist of heterogeneous networks of connected ''strands'', 15 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,14During the waiting time, the gels actually become stronger before undergoing catastrophic failure. 16Given 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][18][19][20][21][22][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.5][26][27][28][29][30][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. 32(This is different from polymer gels, where entanglements dominate the rheology.)][35][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?0][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. 16This 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. 43Like 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, 34 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][45][46] a new characterisation at the single-particle level; and measurements of the local stress. 47hese 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][45][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 s exceeds the yield stress s Y in the neck, 28,38 while remaining below s 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 topdown continuum modelling approach 46 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 47 reveal a complex pattern of residual stresses [48][49][50] that come from the initialisation of the strand; 51 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,53These 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. 56Some 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 B T. This leads to gels with thick strands, in which significant plastic deformation occurs before breakage. 59This 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. 23verall, 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 coarsegrained modelling of macroscopic gels, is vital for the design and control of gels, as an important class of soft materials.
with short-ranged attractive interactions, and two different particle sizes, to suppress crystallisation. 61Particle i has position r i and velocity v i = : r i , it evolves by Langevin dynamics in a simulation box with periodic boundaries: where m is the particle mass, V is the interaction energy, l 0 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. 61 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 8 = 30l and we vary its initial radius r 0 to mimic strands of different thicknesses.An important time scale is the Brownian time t 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 4 t b ) before excising it.We then also allowed the surface of the excited strand to relax for 10 3 t b before we started to deform the simulation box.We varied this equilibration time up to 10 5 t 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.
We deform the strand by affine elongation of the simulation box, so the length L 8 (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 E between the particles a rescaled elongation rate _ g, the strand thickness r 0 / l, and a solvent damping parameter l.Gelation is associated with metastable gas-liquid phase separation, 14 whose critical point is at E Ã % 2:8. 62The data of Fig. 1 have E ¼ 4:5, corresponding to E E = Ã % 1:6, well inside the regime of gelation.We fix l = 10 throughout, consistent with a gel in a high-friction solvent environment.For particulate gels that form in the absence of any solvent, 7 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 8 (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 63 ).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.
5][66][67] The resulting stress estimates are noisy, due to rapid thermal fluctuations, so all such measurements are averaged over a time period of 250t 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 B T/l 3 .
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 s(r), which reveals additional information about the elongation process.To make contact with continuum models of rheology, [44][45][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 s zz (Z), it is related to the force as s zz (Z) E F z /A(Z) where A(Z) is the crosssectional area.After averaging the stress over Z, we plot the result in Fig. 1(f), for several different values of the elongation rate _ g.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,46efore 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).

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][45][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 8 (t) is a rescaled co-ordinate along the strand.Modeling this strand as a thin filament, the theory [44][45][46] is based on four Z-dependent quantities: the filament's crosssectional area a(Z); the local strain rate _ g L (Z); a scalar field W(Z) that coincides (in this case) with the stress s zz (Z); and an internal field w(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, 46 but with the simplification that stress can be captures with a scalar field 45 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 q t a = (_ g À _ g L )a, and the force balance condition is div s = 0 (see Appendix C).For the filament, this means that @ @Z ðaWÞ ¼ 0: 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, 46 which depends on the stress and the plasticity w.
We assume that plastic flow only takes place above the yield stress s Y , taking p(W,w) For the homogeneous solution, we first assume that the field w is a simple constant, independent of both t and Z.One finds W(t) = G_ gt + W(0) for short times (such that W o s Y ) while for long times . 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 w to depend on time with a relaxational dynamics that is controlled by the plastic relaxation itself: 46 @w @t ¼ WðtÞ s Y pðW; wÞ w 1 À wðtÞ ½ : The additional parameter w N fixes the steady-state value of w.This model (with time-dependent w) 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 w 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 _ gt p { 1, the field w is constant until W reaches the yield stress, at which point w relaxes quickly to w N .The slowest elongation rate in Fig. 1(f) is consistent with this regime, fitting gives s Y E 3.3 (recall the units are k B T/l 3 ).For faster elongation rates, the relaxation of w 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, w N = 80 and t ref = 180t b , as well as the initial condition w(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 _ g (which ensures a stress plateau as seen in Fig. 1(f)), and a non-trivial time dependence of w (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. 45ore 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,70In 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 s 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 s 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 w, 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 nontrivial distribution in Fig. 1(g).Such fluctuations are suppressed as E increases and the system approaches the athermal limit.

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 (particlelevel) structure.
We begin with a comparison of dynamical quantities.We define a correlation function C B (Z;g,Dg) that measures how much particles' local environments have relaxed in segment Z, for the period between strain g and strain g + Dg (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 g 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 Dg.On the other hand, Fig. 3(b) fixes the strain lag Dg = 0.05 and varies g, showing how the difference in relaxation rate grows, as g 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 w proposed in the continuum theory.
We also analyze the local structure of the strand using the topological cluster classification (TCC). 71This provides a detailed characterisation of local packing, by identifying specific geometrical structures within the system.We write n n (Z;g) for the average number of neighbouring particles around a particle, at position Z and strain g.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, 72 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 g 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 g, until the system starts to break around g = 0.36.(The large error estimates in this region arise because breakage occurs at different values of g 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,74Hence 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 (g E 0.1).4][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).

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 s 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 s Y .Hence these measurements show directly the mechanism by which plastic events concentrate in the neck.
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 s zz , averaged over small regions of size (1.5l) 3 .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, s zz has strong inhomogeneities: its average is positive (corresponding to tension), but there are significant regions where s zz o 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 s = 0), the residual stresses are divergence-free and arise from the inhomogeneous (amorphous) microstructure of the strands.Since s 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 w 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 s: 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 s zz , within the neck and the bulk.One sees that the probability to see s zz o 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 s zz in the neck (the typical value of s 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,77Overall, 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.

Robustness of results for different model parameters
The results presented so far focused on representative parameter values: attraction strength E ¼ 4:5, strain rate _ g = 10 À6 and strand thickness r 0 = 4l.However, the general picture that we present is robust.This is shown in Fig. 5 where we vary the strain rate _ g, the strand thickness r 0 , and the cohesive energy E. Varying the strength of attractive forces [Fig.5(a and b)], there is weak dependence of the neck thickness N m on E. 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 E, 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 E ! 1.
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 _ g 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,79The 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 _ g (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 0 [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.

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. 32e 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,46This 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,80As these structures are associated with rigidity, 74 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 s zz (Z) has a value that remains close to the yield stress s Y , as the neck develops [Fig.4(d)].The constant tensile force in the strand means that other (thicker) parts of the strand have s zz (Z) o s 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 continuumlevel plasticity field w 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, 29 or to analyse the breakage of strands during coarsening. 43In 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 83 and confocal microscopy observations. 65inally, 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. 84

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,61 where is the mean of the diameters of particles i and j, and a sets the attraction range (we take a = 25l À1 ).While this is longer-ranged than the ''Baxter limit'' of vanishing range, 2 different ranges can be mapped by a law of corresponding states 62 and the dynamics can also be scaled. 85he equation of motion is (1), in which the solvent force is n where T is the temperature, and n a unit white noise.The particles move in a periodic simulation box of dimensions L > Â L > Â L 8 with L > = 20l and L 8 (t = 0) = 30l.The Brownian time is t b = l 0 l 2 /(24k B T). Two natural dimensionless parameters of the model are strength of the depletion attraction, and l ¼ l 0 l= ffiffiffiffiffiffiffiffiffiffi mkT p which measures the solvent damping.In addition, the system undergoes elongation at shear rate _ g 0 .This yields an additional dimensionless parameter _ g = t b _ g 0 .In practice, the strain is increased stepwise, with a fixed time period of Dt = 250t b between steps, so the strain in each step is _ g 0 Dt.When measuring the stress (see below), we average the results over the period Dt, 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 86 in which the natural time scale is t 0 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ml 2 =k B T p , this means that t b /t 0 = l/24 E 0.42 for the parameters used here.The integration time step is 10 À4 t 0 .The natural unit of both pressure and stress is k B T/l 3 ; 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 c between 0.1 and 3 mm. 61,87This allows for a variation of corresponding Brownian times between t b B 0.005-7 s and pressure units between k B T/s 3 B 10 À4 -10 1 Pa.

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 E ¼ 0:01 and a constant pressure.We set P 0 = 0.16 and slowly increase the attractive strength to E ¼ 2:5, using steps of DE ¼ 10 À6 every 2t b .This causes the volume fraction to increase to f E 0.59, the system remains homogeneous because this isobaric transformation not enter the spinodal decomposition regime. 87We allow this dense homogeneous fluid to relax for a time t glass .
We then switch to the NVT ensemble and instantaneously adjust E 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 t 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 t glass up to 10 5 t b , the results depend very weakly on this parameter.Results are shown for t glass = 10 4 t b .) The result is a homogeneous glassy state with volume fraction f E 0.59.Our initial strand is obtained by excising a cylinder of radius r 0 from this system, after which we run dynamics for approximately 1000t 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.

Stress measurement
Our measurements of local stress use a volume-averaged representation of the Irving-Kirkwood stress. 64,65,89Write 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 O whose volume is |O|.The mn component of the IK stress for that region is where W i,O = 1 if particle i is in O and zero otherwise; similarly j ij,O is the fraction of the straight line connecting particles i, j that lies within O.
Taking O to be the entire simulation box O* gives the total stress S, which can also be computed from the virial.The tensile force is then For a local measurement of stress at point r, we take O to be a small cube of side l IK , centred at r.The resulting stress is denoted by s(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 > 2 L 8 /n seg .Then define N Z as the number of particles in segment Z. Taking O 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 > 2 ).
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 O Z is the segment of the box at position Z and V s = N Z pl 3 / (6f) is the estimated volume occupied by the strand, within that segment.(Here f = 0.59 is the volume fraction within the strand, so pl 3 /(6f) 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 s(r) = 0 at every point r, that is, X m2fx;y;zg @ @r m s mn ðrÞ ¼ 0 which holds for n = 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. 65) The method-of-planes method 65 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 volumeaveraged method stress gives much better numerical agreement with (8) than the method of planes.

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 s is proportional to the identity, as would be expected in the bulk of a simple fluid.In a region O with large anisotropic stresses then J 2 will be large.The von Mises criterion 81,82 for failure of solid materials states that breakage will occur when the local J 2 exceeds a threshold.
For strands under tension with homogeneous stress, the dominant element of s O is s zz E s zz (Z), leading to J 2 p s zz (Z) 2 .However, in an strand like the one in Fig. 4(e) with large residual stresses, typical elements of s O have absolute values larger than s zz (Z), leading to a much larger value of J 2,O .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 3 , within segment Z. Since J 2,O 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 crosssectional stress s zz .)

Bond-breaking correlation function
To measure local particle rearrangements, we define b ij (D;g) = 1 if particles ij are within a distance D of each other, when the accumulated shear strain is g.Then the fraction of neighbours of particle i that are lost between strains g and g + Dg is 90 We take D = 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 (g,Dg) 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,91Interestingly, 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.
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 3 t b (or g = 0.01 with the applied strain rate _ g = 5 Â 10 À6 ).Thus we estimate our thermally-activated surface processes rate as 0.5 Â 10 À3 t b À1 , comparable to the strain rate of the black line in Fig. 5.

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. 71These quantities were then averaged over particles in a suitable region (typically the segment of the system with position Z).
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 s(r) = 0. To see this, we compute the (local) tensile force F z (Z) in the strand by taking O in (6) as a segment at position Z.Force balance requires that qF z /qZ = 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.
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: and In these sums, both m and n run over the three Cartesian directions x, y, z, but the term m = n = z is excluded.It is well-known that microscopic expressions for the local stress tensor are not unique. 65We briefly discuss three aspects of this issue.First, since these systems have significant residual stresses, the value of s depends on the scale at which it is measured.Our local stress tensor is measured by taking O 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 = 250t 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 t 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. 66As 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.
We first consider effects of volume-averaging, over cubes of volume l IK 3 .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 B l IK À1 for small l IK .Fig. 8(b) plots l IK % S(Z) for various sizes l IK .Increasing the box size up to l IK E 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 physicallyrelevant residual stresses.
For time averaging, the picture is simpler.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(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).

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  We assume the stress evolves with an elastic loading term and relaxes due to local plasticity events which yields @W @t ¼ G _ g À pðW; wÞ ½ ; In the simplest case we take w as a constant parameter, but this is not sufficient to capture the stress-strain relationship observed in simulations.Instead we again follow 46 in promoting w to a dynamical variable whose relaxation is controlled by the plastic time scale t p .Specifically, we take @w @t ¼ W s Y pðW; wÞ½w 1 À w: where the parameter w N sets the steady-state value of w.

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 w are independent of Z.At early times, we have W(t) = W t=0 + G_ gt, which holds for times short enough that W o s Y .For long times, we have W(t)s Y /(1 À _ gt p ) where t 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 s Y and the plastic activity starts to occur, does w relax quickly to w N , 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 s zz (averaged over Z) and we fit the data at early times to s zz = W t=0 + G_ gt, which yields G = 420 and W t=0 = 0.5 [recall that the units of both these quantities are k B T/l 3 ].We then consider the plateau of s zz (before significant necking occurs): we fit the plateau height to W = s Y /(1 À _ gt p ) which yields t p = 180t b and s Y = 3.3.Note however that this t p is the steady-state value of the plastic time scale: since t p = t ref e 1/w then this is not itself a parameter of the model.To obtain values for t ref and w N we fit the stress overshoot that occurs in Fig. 1(e) for _ g = 5 Â 10 À4 : this yields t ref = 180t b , and w N = 80, as well as the initial condition w 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 (_ g L , a, W, w).We consider a perturbation to the homgeneous solution at wavevector q, that is where the 0-subscripts indicate the solution to the homogeneous equation.Following, 46 eqn (A7) fixes _ g L in terms of a and its time derivative.Linearising eqn (A8) determines da in terms with stability matrix M ¼ W 0 @ W pðW 0 ; w 0 Þ À pðW 0 ; w 0 Þ Àa 0 @ w pðW 0 ; w 0 Þ À W 0 a 0 @ W _ wðW 0 ; w 0 Þ @ w _ wðW 0 ; w 0 Þ (A14) where the notation _ w is a shorthand for the right-hand side of (A11), interpreted as a function of (W 0 , w 0 ).
For W 0 o s Y , there is no plastic flow and M = 0: the strand responds elastically and preserves its shape under elongation.However, for W 0 4 s 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.

Fig. 1
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 (g = 0) with uniform thickness, red arrows show the direction of elongation.Parameters: E ¼ 4:5, r 0 = 4l and _ g = 5 Â 10 À6 .(c) The same strand at strain g = 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 g.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 s zz , averaged over z, for different varying strain rates _ g. (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 g = 0.

Fig. 3
Fig. 3 (a) Dynamical correlation function C B (g,Dg) on varying Dg, 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 g at fixed lag Dg.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).

Fig. 4
Fig. 4 (a) The Z-dependent stress s zz at initialization.(b) The stress at g = 0.4.All particles in a discretized Z-segment are coloured to the corresponding s 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 s zz (Z) at the neck (Z = 0.5) and in bulk (Z = 1) as a function of strain g. (Data averaged over 6 runs.)Colored area is the standard deviation of different runs.(e) and (f) The local stress s zz , for the same configurations shown in (a) and (b).(g) The anisotropic stress measurement % S(Z) as a function of strain g.(h) The stress distribution curves just before breakage (g = 0.35-0.4)for the neck and the bulk, compared with the unstrained distribution.Measurement is done over a single run.

Fig. 5
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 E is varied; in (c) and (d) it is the elongation rate _ g; and in (e) and (f) it is the strand radius r 0 .All lines are averaged over 6 runs with identical starting conditions for a given set of parameters.The blue lines correspond to baseline parameters: E ¼ 4:5, r 0 = 4l and _ g = 5 Â 10 À6 .

Fig. 7 (
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 S 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.65We briefly discuss three aspects of this issue.First, since these systems have significant residual stresses, the value of s depends on the scale at which it is measured.Our local stress tensor is measured by taking O 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 Taking an average over T ave = 250t 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 = 250t 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 \1000t b ).Hence we choose T ave = 250t b for our measurements, as a sensible compromise between time-resolution and noise reduction.

Fig. 7
Fig. 7 (a) The tensile force FðZÞ ¼ Ð x Ð y s zz dx dy calculated from the Irving-Kirkwood stress and the localized virial stress.Data correspond to a single trajectory of Fig. 1(d) at g = 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).

Fig. 8
Fig. 8 (a) The average cross-sectional tensile stress s zz , measured as a volume average over a region of size L 8 2 Â 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 E 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 s zz increases (at Z E 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 = 250t 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 E 0.35; panels (c) and (d) are from a different trajectory where the neck is at Z E 0.2.Both trajectories have the same parameters, which are those of Fig. 1(c and d).

Fig. 9
Fig. 9 Distributions of the local stress s zz .(a) Comparison of the distributions for positive and negative stress values, before any elongation has occurred.Data is found by measuring for 10 6 t b without elongation.(b) Stress distribution evolution far away from the neck (Z = 1) near the start of the simulation (g = 0-0.05),just after the system has reached its stress plateau (g = 0.03-0.05),and just before necking occurs (g = 0.35-0.4).Data has been averaged over 6 trajectories from Fig. 1 in the main text.
A9)where G is the elastic modulus and p the rate of plastic relaxation, which depends in turn on an internal plasticity field w (see below).More specifically, we follow ref.46: p is zero when W is less than a yield stress s Y , and w 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.68 ) Hence,pðW; wÞ ¼ W À s Y t p ðwÞs Y y W À s y À Á (A10)where y is the Heaviside (step) function and t p (w) = t ref exp(1/w), where t ref is a reference time scale.

Fig. 12
Fig. 12 (a) and (b) The cross-sectional stress s 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.
The general approach of ref.44and 45 starts with homogeneous solutions to the equations of motion (that is, a, _ g L , W, w independent of Z).Necking is a linear instability of this solution.
68ere t p is a plastic time scale and y is the Heaviside (step) function.Following ref.46, the plasticity w is analogous to a local temperature that determines the probability of an activated rearrangement, inspired by the theory of shear transformation zones.68Hence,t p (w) = t ref exp(1/w) where t ref is a reference time scale.