Aging of living polymer networks

Microrheology experiments show that viscoelastic media composed by wormlike micellar networks display complex relaxations lasting seconds even at the scale of micrometers. With a novel coarse-grained model of soft matter networks we are able to simulate a thermal quench spanning many decades, from microseconds up to tens of seconds. After mapping our model to real units in agreement with experimentally determined scission energies, we characterize the relaxation process in terms of power-law decays of branches and dangling ends within the network of living polymers. The apparent scale-free nature of this stagnant process complicates the quantification of experimentally relevant time scales as the Weissenberg number.

Introduction -Many soft materials with important industrial and medical applications are formed by thermodynamic self-assembly of elementary constituents dispersed in aqueous solution [1].This process gives rise to networks of linear or branched fibres that can continuously break, rejoin or get entangled.These "living" materials [2] have remarkable viscoelastic properties that are typically studied either by imposing a shear stress with a rheometer [1,3] or by dragging a micro-sized probe through the medium as in active microrheology [4,5].The response and relaxation dynamics of living polymers are the result of large-scale reorientation and mutual reptation of the fibres that, in turn, depend on small scales mechanisms such as scission and rewiring of branches.This multiscale process can give rise to global relaxation dynamics that, for systems as wormlike micelles, has been measured to be of the order of seconds [4,5].
An important issue is the relation between the macroscopic dynamical response of the system, measured experimentally, and its less accessible local structural rearrangements.From the theoretical side, mean-field theories and microscopic constitutive models have been successfully developed [2,[6][7][8][9] to rationalize the outcomes of rheological experiments.However, these studies often dealt with small temperature jumps, which may not lead to the significant disturbance of the micellar networks as, e.g., in microrheological experiments at large Weissenberg numbers where a microbead is moved in the medium at high velocities [10,11].At the same time, numerical studies [12][13][14][15][16][17], performed at a coarse-grained molecular level on specific systems (cationic cylindrical micelles and cetyltrimethylammonium chloride micelles) have either looked at the competition between shear flow and internal network dynamics or estimated the strands scission and branching free energies of these systems.However, these simulations, being still quite detailed and system-specific, cannot have access to the scales spanned by rheological experiments (of the order of seconds and above µm 3 ) or address universal features of the phenomenon.
A generic understanding of the transient behavior of living polymers can be sought by addressing the following questions: which are the dominant local structures (or "mi-crostates") that compete dynamically in the rearrangement of a relaxing living random network?How does the long-time relaxation depend on the mechanisms of local rewiring and scission of the fibers?Is the time dependence of relevant observables undergoing relaxation characterized by simple exponentials or is it more similar to a dynamical scaling with power laws?Simulations of coarse grained models seem an ideal tool for answering these questions: they allow a direct investigation both at short/small and long/large scales and they may adopt a well controlled protocol for driving the system out of equilibrium.
In this Letter, we introduce a coarse-grained description of a thermally driven self-assembly process generating a living network of fibres by means of the merging, scission and rewiring of basic units that represent a portion of the tubular structure.By using large-scale Brownian dynamic simulations and analytical arguments to estimate the free-energy barriers between the competing "microstates", we characterize the underlying microscopic dynamics that leads to the global relaxation towards the formation of the random networks.To go beyond the regime of small perturbations a wide temperature quench is imposed on the system.
Model -The model is quite generic and simple enough to have access to long time scales of the relaxation process.Nevertheless, it is introduced by having in mind wormlike micelles, i.e. the most studied example of living polymers, for which a wealth of experimental data is nowadays available.
A segment of a wormlike network, including dozens of surfactants as those sketched in Fig. 1(a), is mapped to a single micellar unit composed by three sub-units in a rigid body: there is a central repulsive core of radius R (red spheres in Fig. 1(b)) and two radially opposite patchy particles (blue empty spheres) fixed at distance λR, with no excluded volume but with a mutual attraction that leads to the aggregation of units at sufficiently low temperatures.As sketched in Fig. 1(b), the units may either form linear strands or branches.By increasing the value of the branching parameter λ, one increases the propensity of the units to form "vertex" or "3 contact" states, rather than the "strands" or "2 contact" states (see Fig. 1(b)).The model is inspired by those often proposed to study aggregation and dynamical properties of models of gels self-assembled by patchy colloidal particles [18].In that case, branching points are generated by explicitly introducing a fraction of colloids decorated by three or more attractive patchy particles.In our model, only bifunctional patchy beads are considered and the local branching emerges naturally from the formation of multiple contacts between patchy particles.Hence, isolated clusters, filamentous structures and branches of a network are an emergent property of the dynamics and thermodynamics of the system.
Repulsive cores interact via a standard truncated and shifted Lennard-Jones potential with characteristic amplitude , while patchy particles experience an attractive truncated Gaussian potential with amplitude A = 40 and range σ = 0.4R [19].The dynamics of the system follows a set of coupled Langevin equations for the interacting units in an implicit solvent at fixed volume and constant temperature T (henceforth measured in units of , with Boltzmann constant k B = 1 and = 1).Numerical simulations are performed with LAMMPS engine [20].For a more detailed description of the model and for a mapping of the simulation units to those of a wormlike micellar system, see [19].
If we assume that the core radius of the elementary unit length is of the order of the average width measured in wormlike micelles, i.e.R ≈ 5nm, the volume 108R×108R×324R in the simulations corresponds to ≈ 0.5µm 3 .The simulation time unit τ is mapped to 2µs in the real system [19], hence its full duration corresponds to a time t = 10 7 τ ≈ 20s of the order of experimental times [5,21].
The initial configuration is prepared at high T and is composed of freely diffusing globular micelles.By abruptly lowering T at time t = 0, units start to aggregate into intricate networks.Their typical structures resulting in our model bear a strong resemblance of those appearing on smaller scales in less coarse-grained models [17].Yet, our runs go up to N = 6000 units (which would correspond to a much larger number of molecules of surfactants, say ≈ 10 6 ) and may explore the equilibration dynamics for a very long time.We describe the time evolution of the forming network by moni- toring the formation/disappearance of its local motifs, such as dangling ends and branched structures.These are classified by looking at the number of contacts shared by their patchy units.A patchy particle is defined to share a contact with another one if the distance between their centers is ≤ R [19].A dangling end contains a unit with k = 1 contacts, while each branch "3" contains three units each with k = 3 contacts, when they are part of a strand in the network (that is, a chain of elementary units with k = 2) which ends either with dangling ends or branching points.The length L of a strand is the number of its units.
Results -In Fig. 2 we report the time evolution of the mean strand L .After a relatively fast growth, L reaches first a plateau that lasts for a decade (cyan band) and then slowly increases as a power law towards the expected equilibrium value (orange band).The metastable state is characterized by an excess of branching points formed during the initial stage (see top-left inset in Fig. 2).Longer strands can be progressively assembled by slow rearrangements of the whole structure.
In order to rationalize this behavior we first consider the thermodynamics of the smallest motifs forming in the network, namely the strands and the branch points (respectively state "2" and "3" in Fig. 1(b)).A zero temperature energy minimization of these motifs reveals where the formation of strands is energetically more convenient with respect to branch points.This is true for sufficiently small values of λ (Fig. 3(a)), so that the attractive patchy units are so strongly hidden inside the bead core that the three-body contact is severely hindered.For the strand we expect an energy ≈ −2A due to a full overlap of patchy units while for the branching point ≈ −3A with a partial single energy overlap A (with |A | < |A|) which tends to A only when λ is sufficiently large.Hence, there is a value of λ 1.86 above which it is energetically more convenient to form branching points rather than strands (see Fig. 3(a)).
In presence of thermal fluctuations, entropy should also play a role and ultimately a comparison between the free energies F of the microstates "2" and "3" must be performed.To estimate the free energy difference ∆F 32 = F 3 − F 2 we have carried out Langevin simulations with parallel tempering sampling technique [22,23] on a system with three units confined within a small fixed volume (see [19] for details).In Fig. 3(b) we report ∆F 32 as a function of T for three values of λ.Entropy differences ∆S 32 are thereby identified by the (negative) slope ∂∆F 32 /∂T , see dashed lines.Note that the combination of the negative entropy gap ∆S 32 with a λ-dependent potential may result in a situation where the most favorable thermodynamic state is T −dependent, as for λ = 1.95,where branch points are stable only at low T .These findings match the prediction [17] that wormlike micelles may form networks with a variety of features, depending on the inclination of their tubular structures to branch.
To see how the local thermodynamic properties affect the dynamics of network formation, we have estimated the escape rates either from the strand "2" or from the branching point "3" to the detached dimer-plus-monomer state "1".They follow a standard Arrhenius law once reported as a function of 1/T .For the branches (Fig. 3(c  the current estimates of the scission energies (or enthalpies) of micellar systems [17,21,24,25].The exponential dependence of the scission times on the inverse temperature is responsible for the typical time τ s of network rearrangements: for λ = 1.75 and T = 1 we find τ s ≈ t 3 2 × 10 4 τ ≈ 10ms (Fig. 3(c)).
The macroscopic relaxation process for t > τ s , manifested by the power-law increase in L (see Fig. 2), is further understood by looking at the dynamical evolution of units classified by their contacts.Let N k be the number of units with k contacts and f k = N k /N their fraction.Due to the core repulsion, for λ 2 a unit may establish at most k = 4 contacts (two per side).
In Fig. 4 we plot the fractions f k vs time for the same system of Fig. 2.After a quick condensation to a network (f 0 → 0 for t → t 0 ≈ 3 × 10 3 τ ), the network does not contain the equilibrium proportion of network motifs.The fraction f 1 of units at dangling ends tends slowly to zero, with a scaling compatible with a power law f 1 ∼ t −1 .Eventually, at t eq ≈ 2 × 10 6 τ , mapped to ≈ 10s, the fraction f 1 reaches its asymptotic value.This very slow reduction of the number of strands with dangling ends highlights a stagnant global relaxation of the network.Also the fraction f 3 of units meeting at branches slowly decreases, f 3 ∼ t −1 (Fig. 4), till t eq and then fluctuates around its equilibrium value.
The slow scalings of f 1 , f 3 , and L suggest that the dynamics within the network proceeds through a sequence of frustrated rearrangements.For instance, a dangling end joining its free cap to the side of a strand would form a branch.This local state is however less favorable (see fig 1b), thus leading to a subsequent decay back to a dangling end and a linear strand.A direct merging of two dangling ends would better stabilize the system but this event becomes more and more unlikely with time due to the concomitant decrease of the dangling ends population.
The variety of micellar structures depends on several parameters (solvent, salt concentration, surfactant chemistry) and includes networks rich in branches.This situation is realized in our model by increasing λ.For example, with λ = 2, the relaxation proceeds from the initial metastable network with too long strands toward one with more branches, as shown in Fig. 5. Again, there is a very slow convergence of the fraction f 1 of units at the caps of dangling ends, see Fig. 5(a).In Fig. 5(b), in linear scale, it is apparent that also the fractions f 2 and f 3 slowly converge to their asymptotic average values, not reached yet at the final time of our run.
The phenomenology described above holds for several micellar densities.Indeed, by changing considerably N and keeping the volume V fixed, we have been able to vary ρ num = N/V by a factor of 3. Yet, the scaling of the mean quantities is quite similar to the previous cases.In particular, the average strand length in the initial metastable network is unaffected by the change of density.However, the density determines the time scale t 0 of the initial network condensation.By scaling arguments, at low concentrations, diffusion leads to t 0 ∼ ρ −2/3 num .At higher concentrations it is rather the "ballistic" expansion of the network core (as in other Brownian coagulation systems [26]) that determines the time scale t 0 ∼ ρ −1 num .In Fig. 3(e) we report estimates ∼ t 0 (up to a constant that is irrelevant for the scaling [19]) which are compatible with the two regimes.
Discussion -We have introduced a coarse grained model to study the relaxation dynamics of living polymer networks driven out equilibrium by a deep thermal quench.A novel relevant feature of the model is its ability to bridge the short time dynamics of its elementary motifs (dangling ends, linear strands, branches) with the long time behavior of the mesoscopic network.Although the model is sufficiently simple and generic, a detailed mapping to micellar systems matches recent experimental results, such as scission free energies up to ≈ 30k B T .
By combining large-scale simulations and analytical estimates of the free energy barriers between motifs, we rationalize the observed slow (i.e power-law) convergence of the network to equilibrium in terms of the frustrated dynamics of the local motifs.In particular we argue that the relaxation process develops through a sequence of rearrangements that are not only energetically challenging but also require the solution of topological constraints.Moreover, the far-from-equilibrium character of our setup is naturally beyond the perturbative regime, where exponential relaxation laws are expected [6,9].While we have used a standard temperature quench to produce a clean representation of the network relaxation, in practice the considerations we have drawn should be extended to transient regimes of mechanical or chemical origin.For instance, a passage of a micro-bead through a portion of a viscoelastic fluid made of wormlike micellar networks may "scramble" severely the system by generating dangling ends and by distorting the branches.This unbalanced state then should re-equilibrate by following a relaxation dynamics akin to the one we have characterized here.Such scenario is relevant for microbeads of diameter D traveling with mean velocity v at high Weissenberg number Wi = v/v r , where v r = D/τ rel is the velocity scale associated to the relaxation over D with time scale τ rel .Our findings suggest that a proper definition of τ rel should consider the possible emergence of power-law behaviors for strongly disturbed polymer networks.In fact, the relaxation times found here are of the same order of magnitude (seconds) of those measured experimentally in viscoelastic micellar networks [4,5].It is thus natural to conjecture that τ rel is the time at which each microscopic motif of the network reaches, through a power-law decay, its typical (equilibrium) fraction.In this respect, it is the ending time of this scale-free aging (not to be confused with chemical aging, as a high-temperature degradation of carboncarbon double bonds [27]) that should be considered for evaluating the standard Deborah and Weissenberg numbers used in experiments [10,11].
Finally, our findings also raise the question on which of, e.g., self-healing rubber [28,29], colloidal systems [30][31][32], or network fluids [33] do experience a relaxation process characterized by a slow rearrangement of local motifs.
Physical units.In our simulations the system temperature T is measured in units of /k B , where is the amplitude of the Lennard-Jones potential, with k B = = 1, A = 40 and σ = 0.4R.Distances are measured in units of monomer radii R. We set m/γ = τ , where τ = R √ m is the (unitary) characteristic simulation time.In order to estimate the time τ in physical units, one can proceed as follows.From the Stokes formula of the friction coefficient for spherical beads with radius R, we have γ = 6πη sol R, where η sol is the viscosity of the solvent.By eliminating m in the above formulas for the characteristic time, one has γ = τ /(2R) 2 , hence τ = 3πη sol (2R) 3 2µs In the numerical evaluation we consider the nominal water viscosity at room temperature, η sol 1 Pa s, and we set T = /k B = 300K and R = 5 nm.
For the large-scale simulations of this work, the system is evolved in a prism of sides (d, d, 3d) with d = 108R, with periodic boundary conditions.Its volume corresponds to 3d 3 0.5(µm) 3 .If N = 6000 units are introduced in this volume, their volume fraction becomes

CONTACT THRESHOLD
Given a configuration of micelles identified by the set of r iα , two distinct units i and j = i are considered in contact if there exists a pair of patchy subunits α and β such that d jβ iα ≤ θ, where θ is a threshold.To determine the optimal value of θ for the potential U specified in the previous section, we have computed the distributions of patchy particles distances d jβ iα for all i = j in generic network of micelles, see Fig. S1 for an example.The distribution displays a clear gap in the region [0.5R, 1.8R] separating the ensemble of bound patchy particles (with d jβ iα < ∼ 0.5R σ) from the ensemble of unbound ones.Accordingly, we have chosen θ = R, well inside the gap.We have verified that this threshold value is reliable in the whole range of λ and T considered in this work.

PARALLEL TEMPERING SIMULATIONS
Free-energy differences between microstates "2" (strands) and "3" (vertices) have been computed through Langevin equilibrium simulations with parallel tempering sampling technique [S2, S3].This method is employed to sample accurately the low-temperature equilibrium distribution of the system, thereby reducing the impact of very long correlation times on computational time.
Simulations have been performed for a system of N = 3 interacting micelles confined in a spherical volume with radius R 0 = 10.85R and volume fraction φ 2 × 10 −3 , which corresponds to the value of micellar density of largescale simulations with N = 2000.We have verified that in this minimal setup the confining volume is sufficiently larger than the typical micellar volume, so that it does not introduce significant finite-size effects.The free energy difference ∆F 32 is computed via the equilibrium Boltzmann distribution ∆F 32 = −T log(P 3 /P 2 ), where P 3 and P 2 are respectively the occupation probabilities of states "3" and "2" sampled over the equilibrium Brownian dynamics.Parallel tempering evolution is performed for a total time of 5 × 10 6 τ with a swapping period of 50τ .Temperature swaps are realized between adjacent temperature states within the set of temperatures considered in Fig. 3(b).

SCISSION TIMES
Scission times in Figs.3(c-d) have been computed for a system of N = 3 interacting micelles confined in a spherical volume with radius R 0 = 10.85R and evolving according to Eq. (S4).For each value of temperature T , micelles were initially prepared in a "3" configuration (Fig. 3(c)) and in a "2" configuration (Fig. 3(d)) and evolved until a scission event to a dimer-plus-monomer state (configuration "1") was observed.For each T , the average scission time has been computed by averaging over a sample of 10 independent scission events.

CONDENSATION TIME SCALES
The condensation time scale reported in Fig. 3(e) is measured as the time needed for the fraction of 0-contacts f 0 to go below the threshold value f * = 0.4 starting from a homogeneous initial condition of isolated micellar units (f 0 = 1) at time t = 0. Micelles are evolved at T = 1 according to Eq. (S4) in a cubic box with side equal to 324R and periodic boundary conditions.Data in Fig. 3(e) refer to system sizes N = 25, 50, 100, 200, 400, 800, 1600 and λ = 1.8.We have verified that f * is irrelevant for the scaling behavior of t 0 with ρ num in so far as f * is large enough to provide a sampling of the initial relaxation stages.Diffusive and ballistic regimes are further analyzed in Fig. S2, where we compare the average particle distance

FIG. 1 .
FIG. 1.(a) Sketch of a wormlike micellar network and of selected portions where each surfactant is depicted as a hydrophilic (red) head and a hydrophobic (blue) tail.(b) The coarse grained versions of the same configurations in our model, with their possible free energy landscape.The inset is an enlargement of a single rigid unit.

FIG. 2 .
FIG.2.Evolution of the average strand length L for a system of N = 6000 micelles with λ = 1.75, and T = 1 (thick black curve).The cyan band marks the metastable regime with excess of branches (see an example in the upper inset) while the orange band identifies the regime of thermodynamic equilibrium (lower inset).Vertical arrows indicate the condensation time scale t0 3 × 10 3 τ ≈ 1ms and the typical time τs 2 × 10 4 τ of network rearrangements.In the two configurations of the network, units with k = 1, 2, 3 contacts are drawn, respectively, in yellow, red and cyan.

FIG. 3 .
FIG. 3. Energetic and kinetic properties of a trimer (inset of panel (a)) (a) Ground-state energy of linear "2" (blue) and triangle "3" (red) configurations as a function of the branching parameter λ.The two lines cross at λ 1.86.(b) Free energy difference ∆F32 vs. temperature T for λ = 1.79, 1.86 and 1.95.For the same values, −∆S32 = 3.2, 3.2 and 5.1 from linear fits, see dashed lines.Open symbols reported on the vertical axis correspond to the energy difference from panel (a), at T = 0. (c) Average breaking time t3 for a trimer initially arranged in the configuration "3" and (d) average scission time t2 for a strand "2", as a function of the inverse temperature 1/T for four different values of λ.Dashed lines are Arrhenius fits of ∆E S 13 and ∆E S 12 .(e) Time scale of the network condensation as a function of the micellar density.The straight lines indicate the scalings expected at low and high densities.

4 FIG. 4 . 1 .
FIG.4.Evolution of the relative density f k of units with k-contacts for a system of N = 6000 micelles with λ = 1.75 and T = 1, as in Fig.1.Dashed lines represent a scaling ∼ t −1 .

FIG. 5 .
FIG. 5. Evolution of the fraction f k of units with k-contacts for a system with N = 2000 micelles and λ = 2, T = 1, (a) in loglog scale and (b) in linear-log scale.The inset shows that the final configuration of the network is rich in branches: gold, red, cyan, and magenta beads identify, respectively, units with k = 1, 2, 3, 4 contacts.

FIG. S1 :
FIG. S1: Distribution of distances between patchy particles for a system of N = 2000 micelles with λ = 1.75,A = 40 and σ = 0.4R.Micelles are prepared in a homogeneous state and evolved at T = 1 according to Eq. (S4) in a cubic box with side equal to 54R and periodic boundary conditions.The histogram is representative of the configuration obtained at t = 10 5 τ , see the inset.Black, red, green and blue bars refer respectively to the combination of indices (α, β) = {(2, 2), (2, 3), (3, 2), (3, 3)} in the set of distances d jβ iα for i, j = 1, • • • , N .Multiple peaks above the gap region derive from the regular spacing of particles in the strands.

FIG. S2 :
FIG. S2: Evolution of the maximal cluster size dc(t) (red solid lines) for N = 25 (panel (a)) and N = 1600 (panel (b)).Black dashed horizontal lines indicate the values of the typical particle distance d while green vertical lines identify the time at which f0 reaches the threshold value f * .