The role of temperature in the rigidity-controlled fracture of elastic networks

We study the influence of thermal fluctuations on the fracture of elastic networks, via simulations of the uniaxial extension of central-force spring networks with varying rigidity, i.e. connectivity. Studying their failure response, both at the macroscopic and microscopic level, we find that an increase in temperature corresponds to a more homogeneous stress (re)distribution and induces thermally activated failure of springs. As a consequence, the material strength decreases upon increasing temperature, the damage is spread over larger lengthscales and a more ductile fracture process is observed. These effects are modulated by network rigidity and can therefore be tuned via the network connectivity and the rupture threshold of the springs. Knowledge of the interplay between temperature and rigidity improves our understanding of the fracture of elastic networks, such as (biological) polymer networks, and can help to refine design principles for tough soft materials.

: Fracture of thermal spring networks. (a) Portion of spring network under 1.5% extensional strain. Several snapshots (corresponding to different shades of gray) are overlaid to indicate the effect of thermal motion. (b) Example of stress-strain response of a diluted spring network (p = 0.65, λ = 0.03, T * = 10 −4 , L = 128). We highlight the peak-stress σ p , the corresponding peak strain ε p and the maximum drop in stress ∆σ max .
networks (see Fig. 1), similar to previous studies [22,9,8,11]. To introduce thermal fluctuations into these systems we perform Langevin Dynamics simulations, in which an implicit solvent keeps the temperature of the system constant. We find that the strength of the networks is dependent on temperature and that the effect of the thermal fluctuations is coupled to the rigidity of the network. The simple structure of the model allows us to highlight the interplay between rigidity and temperature, and to provide insight in the underlying microscopic mechanisms of stress homogenization and diffuse failure.

Model and methods
We consider diluted spring networks with a 2D triangular topology consisting of L × L nodes separated by a distance ℓ 0 . Nearest neighbors are connected by bonds, which gives a maximum network connectivity z max = 6. The network is subsequently randomly diluted by removing a fraction 1 − p of the bonds, such that the average connectivity becomes z = p z max . Periodic boundary conditions are employed in all directions. The bonds are harmonic (linear) springs with spring constant µ and rest length ℓ 0 . Excluded volume interactions are not present in the system. During fracture simulations, bonds break irreversibly when their deformation ∆ℓ exceeds a rupture threshold λ, that is the same for all the springs. We will focus on networks with λ = 0.03.
Simulations are performed using LAMMPS [26] and nodes follow Langevin dynamics: where F = µ∆ℓ, m is the mass that is set to unity, ζ is the friction coefficient related to the (implicit) solvent viscosity, k B the Boltzmann's constant, T the temperature and R(t) white noise with zero-mean. The integration time step is set to δt = 0.001τ , where τ = mℓ 2 0 /E is the unit time of our simulations, and the energy scale is set to E = 1. The spring stiffness (in reduced units) is set to µ = 1000. In our analysis, we will use the reduced temperature T * = k B T /µℓ 2 0 , indicating the ratio between thermal and elastic energies. We fix ζ = 10 (by setting the damping factor in LAMMPS to 0.1). In addition, we can define a relaxation time scale τ relax = ζ/µ, the time that a node requires to travel a distance ℓ 0 if it is subjected to a force µℓ 0 . For our typical set of simulations τ relax = 5 · 10 −3 τ .

Measuring linear modulus and non-affinity
We calculate the linear Young's modulus E from the difference in average stress at 0% strain and 1.5% strain. At both strain values the system is equilibrated for 100τ and averaged over 1900τ . The stress σ is defined as the yy-component (along the deformation axis) of the virial stress tensor normalized by µ. We also calculate the non-affinity parameter (at 1.5 % strain) defined as where r is the time averaged position of the individual nodes after an applied deformation of 1.5% strain and r aff the position assuming an affine displacement of 1.5% strain with respect to the time averaged position of individual nodes at rest. ε is the strain and . indicates the average over all nodes. For these simulations, both non-percolating clusters and primary dangling ends (i.e. nodes that are only connected to one bond after dilution) are iteratively removed from the network.

Non-linear elasticity and fracture
The network is uniaxially deformed in the y-direction up to 100% strain with a fixed strain rate ofε = 0.001τ −1 (i.e. 0.0001% deformation per unit time) while the lateral size is kept constant. We remap the node positions between time steps, temporarily enforcing affine deformation. The deformation is relatively slow compared to the relaxation time,ε = 5 · 10 −6 τ −1 relax , which indicates the system has time to respond to the affine deformation via structural rearrangements. Results for varyingε are reported in the ESI. A quick equilibration run of 50τ precedes the deformation. The network response is quantified by looking at the stress σ as a function of strain ε. In addition, we follow the instantaneous non-affine response where the affine response is calculated with respect to the equilibrium position of the nodes at 0.0% strain (averaged over 1900τ ). Please note that Γ is only based on the instantaneous positions and the non-affine response is therefore a combination of both rigidity controlled non-affine network rearrangements and instantaneous thermal fluctuations (see ESI for details on the relation between Γ and Γ mech ).
In fracture simulations bonds are broken every 100 steps (i.e. 0.1τ ) when the connected nodes are separated by a distance more than ℓ 0 +λ. From the measured stress-strain curves we extract several quantities (see Fig. 1). All quantities are averaged over several configurations and expressed in reduced units. The peak stress σ p is defined as the highest measured stress, and the peak strain ε p is its corresponding strain value. The maximum stress drop ∆σ max is calculated according to a procedure [27,11] where we i) calculate the derivative of the stress-strain curve, ii) make a list of consecutive data points which have a negative derivative and note the initial and final strain of each interval, iii) calculate the stress drops by subtracting the stress at the final strain from the stress at the initial strain, iv) identify the largest stress interval, which corresponds to maximum stress drop ∆σ max . For the stress distribution analysis, we make instantaneous histograms of the bond lengths ℓ i during the simulation at every percent strain. Based on these histograms, we calculate the excess kurtosis where s is the standard deviation of the histogram, N b the total number of bonds, and ℓ the average bond length. For all these parameters the standard error is calculated as the standard deviation divided by the number of sampled configurations (standard error of the mean). The errors are shown when they are larger than the symbols displayed in the graphs.

Results and discussion
Using Langevin Dynamics simulations, we uniaxially deform diluted triangular central-force spring networks to study both linear and non-linear network mechanics. By analyzing the network response on both a macroscopic and microscopic level, we gain insight into the effects of thermal fluctuations on the fracture process of networks as well as how the impact of thermal fluctuations is controlled by the network rigidity.

The effect of thermal fluctuations on linear elasticity
First, we study the effects of thermal fluctuations on the linear elasticity. In Fig. 2(a) we plot the linear modulus E of the network as a function of the network connectivity factor p for several reduced temperatures T * . The linear modulus E describes the resistance of a network to deformation and we observe that networks below the isostatic point of mechanical stability [12] (i.e. networks below p iso ≈ 0.66) display a finite linear modulus E, which would be absent for athermal systems (in the limit of T = 0). This finite E is an effect of entropic stiffness, a temperature-dependent phenomenon. As reported in literature [22,9], the scaling of the linear modulus with temperature E ∝ T α depends on both connectivity and temperature itself. By plotting E as a function of T * we extract the scaling exponent α from a power-law fit for three different values of p, as shown in Fig. 2(b). For a sub-isostatic network with e.g. a connectivity parameter p = 0.55 the linear modulus scales with α = 0.84, which roughly corresponds to the dependence found in the anomalous regime as defined in Ref. [22], where a shear deformation was instead considered. It was argued that the disordered network structure causes this sub-linear dependence. Whilst there is a clear dependence of the linear modulus on the temperature below the isostatic point, the curves for the different temperatures start to converge when approaching a structurally rigid network ( Fig. 2(a)). Accordingly, stiff networks display temperature insensitivity (α ≈ 0), as can be seen for a network with p = 0.70 in Fig. 2(b). As T * increases, however, the network connectivity becomes less important as the energetic contribution arising from the structural rigidity becomes negligible compared to the entropic elasticity. This is noticeable in Fig. 2(a) where the curve for T * = 10 −2 is roughly flat for the entire p-range, and also in Fig. 2(b) where for p = 0.70, E increases for T * > 10 −3 . As predicted in Ref. [22], we also find a different scaling for networks close to the isostatic point, see e.g. curve for p = 0.62 in Fig. 2(b). Although the exponent α is slightly different from the findings of Ref. [22] (where shear deformation and different simulation methods were employed), we were also able to obtain critical rescaling as shown in Fig. 2(c). We can conclude that there are different regimes of dependence for the linear modulus on the temperature based on both rigidity and temperature.
Furthermore, we find similar rigidity dependent behaviour of the thermal fluctuations in the non-affinity parameter Γ mech (Eq. 2), reported in Fig. 2(d) as a function of p for different T * . The non-affinity of the network describes how much the time-averaged local deformation differs from the global (externally imposed) deformation. At low T * we find a peak in non-affine deformation around the isostatic point (p ≈ 0.66). This peak arises from the tendency of the spring network to minimize internal stress upon deformation. If the spring network is far below the isostatic point, the stress can be reduced significantly by a small amount of non-affine rearrangements while at the isostatic point many non-affine rearrangements are required. At the isostatic point, an increase in T * decreases Γ mech , which suggests that thermal fluctuations act as a stabilizing field, similar to the bending rigidity in fiber networks [15]. However, we note that the effect of thermal fluctuations is always present, even without external deformations, leading to structural rearrangements in the rest state (see ESI). Above the isostatic point, we observe that the non-affinity converges for most values of T * (see ESI for details), which indicates that above the isostatic point the network rigidity dominates the non-affine response. Only if T * > 10 −4 , we see that thermal fluctuations affect the non-affine response, increasing Γ mech . This is in contrast to fiber networks, where the non-affinity decreases with an increase in bending rigidity. We hypothesize that this difference occurs because in the case of fiber networks the fibers have a preference to remain straight to minimize stress caused by fibre bending, while in the case of thermal fluctuations an affine displacement of the nodes will not minimize the stress caused by the randomly oriented thermal fluctuations. Below the isostatic point, the effect of thermal fluctuations on the non-affine response is significant. We observe that at T * = 10 −8 the non-affine response is the smallest and that a moderate increase in T * up to T * = 10 −6 leads to an increase in the non-affine response, corresponding to what is observed for fiber networks. However, we also observe a decrease in Γ mech if the temperature is increased beyond T * = 10 −6 , which is not observed in fiber networks. It is unclear if this deviation is caused by a fundamental difference between thermal fluctuations and bending rigidity as a stabilizing field or that longer equilibration times are required to gain quantitative information on the non-affine response in this regime (see ESI for details).
In general we find that at a global level thermal fluctuations act as a stabilizing field, dampening rigidity dependent behaviour around the isostatic point. However, our results suggest that the random nature of the thermal fluctuations causes significant differences in the local response with respect to stabilizing fields in athermal systems such as bending.

The effect of thermal fluctuations on non-linear elasticity
In the previous section, we have shown that thermal fluctuations rigidify sub-isostatic networks. Here, we analyze the non-linear elasticity of unbreakable networks to quantify the effect of temperature when networks become more and more strained. The strain-stiffening observed for sub-isostatic networks in the athermal limit has been extensively studied [28,29,30,17,31,16,18,14]. In Fig. 3(a), we report the stress-strain curves for a network with p = 0.56 at different temperatures. It is evident that the network strain-stiffens for all the T * investigated. We observe that the onset of strain-stiffening is barely dependent on temperature. Furthermore, the stress response becomes independent of T * at high strains, similarly to what has been observed for other stabilizing fields, e.g. bending [29,17].
Signatures of strain-stiffening can also be observed in the non-affine response of the network. In Fig. 3(b), we report the instantaneous non-affinity parameter Γ, that intrinsically includes both the non-affine contributions from instantaneous thermal fluctuations and structural rearrangements. As a result, high non-affinity values can be observed at low strains, where the size of the non-affine thermal fluctuations is large compared to the applied strain. At low temperatures, a peak can be observed in the non-affine response around the onset strain. At high temperatures, this peak is overshadowed by the non-affine thermal fluctuations. At high strain, the network elasticity is controlled by stretching of the bonds and the network response becomes increasingly affine for most temperatures. Only at T * = 10 −2 , the non-affine fluctuations are still visible.
To disentangle the effects of temperature and network connectivity, we normalize the stress-strain curve with the ones obtained in the athermal energy-dominated limit. In particular, we plot the stress ratio σ/σ ath in Fig. 3(c), where we used the data obtained at T * = 10 −8 for σ ath . A ratio of σ/σ ath ≈ 1 implies that the mechanical behaviour is basically insensitive to variations in temperature. As can be seen in Fig. 3(c) for a network with p = 0.56, there is a regime of strain in which the stress ratio depends on T * that decreases upon stretching the network more and more. At increasing temperature, this stress ratio is both higher at the start and approaches temperature insensitivity at a higher strain. The start of decrease in stress ratio for all temperatures occurs at approximately the same strain value, corresponding to the onset of strain-stiffening. This transition could therefore be interpreted as a transition between a regime dominated by thermal fluctuations to a regime dominated by bond stretching. This is analogous to the bending-to-stretching transition observed in fiber networks [28,32,29]. We summarize these observations in a mechanical phase diagram sketched in Fig. 3(d), where we can distinguish two regimes: a mechanically-dominated regime (blue) where structural rigidity overpowers the effect of thermal fluctuations and a temperature-controlled regime (orange) where thermal fluctuations play a more important role in the elastic behaviour. The transition between these regimes depends on the reduced temperature T * (and therefore both on the actual temperature T and the bond stiffness µ), the connectivity parameter p and the strain ε. This transition is in general very gradual as can be seen from the two cross-sections of the mechanical phase diagram reported in Fig. 3(e-f), where we show the stress ratio obtained by some of our simulations. When T * is fixed (Fig. 3(e)) and we increase p, we observe a steep decrease in the strain associated to the thermal-stretching transition. Above the isostatic point, the mechanics of the rigid networks is barely affected by thermal fluctuations at this temperature. In Fig. 3(f), we observe that with increasing temperature the stress ratio increases but the strain characterizing the transition seems to reach a limiting value. This limiting value is a result of the onset of strain-stiffening, which is independent of temperature and corresponds to the transition to the elastic regime.
In summary, we identified a rigidity-dependent transition between two regimes where thermal fluctuations are or are not important. In the following sections, we will investigate whether this underlying transition also influences the fracture of these elastic networks.

The effect of thermal fluctuations on macroscopic fracture
Under athermal conditions, bonds break only after the onset strain, as only at this stage the bonds are under tension. Furthermore, Fig. 3 indicates that the contribution of the thermal fluctuations to the stress in the system is significantly reduced beyond the onset strain. Does this mean that there is only a minor influence of temperature on the failure response?
We first focus on macroscopic descriptors and characterize the stress-strain curves obtained from fracture simulations. In Fig. 4(a-b), we show the response of two representative networks with small rupture threshold λ = 0.03 and different connectivity at several temperatures T * . For the network with p = 0.65 ≃ p iso (panel a), a clear decrease in peak stress σ p for increasing T * is observed, while a variation in the peak strain ε p is less evident as the fracture becomes more ductile and the decrease in stress after the peak is less pronounced. For the very rigid network (p = 0.90, panel b), the decrease in both σ p and ε p is clearly observed. Similarly, the fracture becomes more ductile for higher T * , even though a clear stress drop is still recognizable at the highest temperature simulated. In both cases, the networks become weaker with increasing T * . Furthermore, when approaching the athermal limit (T * → 0) the peak stress becomes less sensitive to variation in temperature. In Fig. 4(c), we show the temperature dependence of σ p for several connectivities with λ = 0.03. The common trend is little variation at low temperatures, almost a plateau that is indicative of approaching the athermal limit, followed by a decrease when temperature is increased, with σ p eventually dropping to zero when a temperature of T * ≃ 10 −4 is reached. On one hand, for low T * the peak stress is evidently controlled by the network rigidity, as previously investigated in the athermal limit [9,11]. On the other hand, when the thermal energy is of the order of 1 2 µ(λℓ 0 ) 2 the network structure is irrelevant, as springs spontaneously break and the system shows melting behaviour. We will later describe the melting point using the reduced quantity k B T /[ 1 2 µ(λℓ 0 ) 2 ] = T * /( 1 2 λ 2 ). In between these limits, there is a cross-over regime. To better assess the role of rigidity in this intermediate regime, we normalize σ p by its value in the athermal limit σ p,ath and plot this ratio in Fig. 4(d). The transition between the athermal limit, where σ p /σ p,ath = 1, and the melting limit, where such a ratio goes to zero, depends on a subtle coupling between connectivity and temperature itself. Far below (p < 0.60) and far above the isostatic point (p > 0.80) the connectivity plays a small role since at every temperature the curve exhibits two plateaus (at small and large p). However, around the isostatic point rigidity and thermal fluctuations are coupled, since at all the intermediate temperatures we can observe a sharp increase in σ p /σ p,ath upon increasing p, connecting the two limiting plateaus. On passing, we note that the plateau for small p is lower, suggesting that temperature starts to affect failure of very diluted networks earlier than for networks with large p. Furthermore, we speculate that the complex temperature-dependence around the isostatic point arises from locally floppy regions that are rigidified by thermal fluctuations (whose magnitude depends on temperature itself) and are therefore able to sustain and concentrate stress, and break. Since the isostatic point marks the onset of mechanical stability, such effect is largest for networks close to it.

The effect of thermal fluctuations on microscopic fracture
Clearly, the failure response is temperature dependent across the entire rigidity range, but the influence of temperature indeed seems to depend on the distance with respect to the isostatic point. Are these differences also apparent at the microscopic level? To investigate this, we monitor the number of broken bonds during the simulations. As shown in Fig. 5(a-b), the fraction of broken bonds φ as a function of deformation indicates that higher temperature leads to earlier and overall increased damage. However, the effect of temperature is more significant close to the melting temperature T * ≃ 1 2 λ 2 = 4.5 · 10 −4 , whereas for lower temperatures the system response is still very much influenced by rigidity. By focusing on the fraction of bonds broken at the peak strain φ p , counting also the bonds broken during the peak event, as shown in Fig. 5(c), the diverging behaviour when approaching melting is evident. This increase in broken bonds could explain the decrease in material strength σ p . Also the fraction of bonds that break above ε p (Fig. 5(d)), the post-peak response, increases close to the melting point, which points towards a prolonged post-peak response, i.e. higher ductility. A direct inspection of the simulation snapshots (Fig. 6) suggests that bonds that break up to the peak strain ε p (red bonds) are dispersed more homogeneously throughout the sample at a higher temperature. The snapshots also reveal a big difference in the response to temperature between networks around (p = 0.65) and far above the isostatic point (p = 0.90). Around the isostatic point, the damage up to ε p is already diffusive in the athermal limit, and its delocalization is enhanced when the temperature is increased. In contrast, the failure response far above the isostatic point shows a clear transition from crack nucleation in the athermal regime to a more diffuse failure response close to the melting point. However, the post peak response at p = 0.90 is clearly still dominated by the propagation of cracks. Nevertheless, at high temperatures we observe the development of multiple cracks, sometimes even not perpendicular to the deformation direction, and evidence of crack merging.
In summary, we show that an increase in temperature leads to an increase in diffuse failure, implying suppression of stress concentration before the peak stress. These observations suggest that thermal fluctuations are responsible for two apparently contrasting effects: on the one hand, they create "instantaneous defects" resulting in more regions with broken bonds, that reduce material strength; on the other hand, the fluctuations allow to delocalize stress away from such defects, delaying the propagation of large cracks. As a result, the damage pattern is diffuse throughout the system.  1024). (e) Peak stress in the thermodynamic limit σ ∞ p normalized by the corresponding value in the athermal limit σ ∞ p,ath , as a function of reduced temperature T * normalized by its melting value 1 2 λ 2 . Error bars represent the standard error in the fit for σ ∞ p . (f) Estimate of the system size where ∆σ max is minimal as a function of T * . In both (e) and (f) the marker shape corresponds to the value of p as introduced in panels (a-d).

Thermal fluctuations increase the length scale of stress redistribution
It is possible that some characteristic lengthscale associated to stress concentration (or equivalently to stress delocalization) exists. For example, we have recently [11] shown that in athermal systems, brittle (abrupt) fracture always occurs for networks above a certain system size L * . This critical size can be tuned by the network rigidity, i.e. by varying p and λ. The onset of the size-induced brittleness can be determined by looking at the non-monotonic size-dependence of ∆σ max . Here we investigate whether temperature affects this critical system size, which can be interpreted as a characteristic lengthscale over which stress concentrates. This analysis is therefore another way to further assess the role of temperature on stress concentration. Therefore, we examine how thermal fluctuations affect the macroscopic fracture descriptors for different system sizes, focusing on the maximum stress drop ∆σ max that quantifies fracture abruptness. In Fig. 7(a-d), we plot the size-scaling of the maximum stress drop ∆σ max (closed symbols) together with the peak stress σ p (open symbols) for four combinations of p and λ at different temperatures. In all cases, we observe a monotonic decrease of σ p as a function of the system size. These trends can be fitted by a power law where σ ∞ p is the failure stress in the thermodynamic limit (infinite system size), β the size scaling exponent and α a fitting constant. In Fig. 7(e), we plot σ ∞ p normalized by its athermal value σ ∞ p,ath as a function of T * normalized by 1 2 λ 2 (see ESI for the other fitting parameters). The observed trend underlying a transition from low T * to melting is consistent with the data at fixed system size and fixed rupture threshold λ presented earlier in Fig. 4. Here, we can also appreciate the effect of varying λ in the intermediate temperature regime. For example, the largest λ = 0.30 (downward triangles, networks with p = 0.90) shows a steeper decrease in the normalized fracture stress, suggesting that thermal effects kick in at higher temperatures for these very rigid networks.
Finally, we focus on the maximum stress drop ∆σ max . From Fig. 7(a-d), we observe that a non-monotonic trend is observed in basically all cases, consistent with our previous results in the athermal limit [11]. We speculated that the initial decrease, implying a more ductile fracture upon increasing system size, is associated to the rupture and reformation of locally stressed regions (often consisting of aligned springs, and sometimes called force chains [33,34,35,9]). However, upon increasing the system size ∆σ max starts to increase, suggesting that stress concentration is present in the system, since it fractures in a more abrupt way. At even larger L, ∆σ max decreases again, now following the same trend for the peak stress σ p that sets the upper bound to the possible stress drop. In Fig. 7(c) the entire trend is visible for the system sizes explored in this work, whereas in the other panels only parts of it are captured. Importantly, for all systems, the trend depends on temperature. In particular, in Fig. 7(f) we quantify the effect of temperature by plotting the system size L min corresponding to the minimum ∆σ max as a function of T * . We observe that thermal fluctuations increase the value of L min , which can be interpreted as a lengthscale for stress concentration. The role of temperature seems particularly relevant at low connectivity, where the stress is already very delocalized in the athermal limit.
In summary, we find that also in the thermodynamic limit there is a crossover from an athermal regime to a melting regime where the failure behaviour is determined by both rigidity and thermal fluctuations. Moreover, thanks to the analysis of ∆σ max , we find evidence that temperature increases the region over which stress is delocalized.

Thermal fluctuations homogenize stress
The delocalization of stress is mediated by structural rearrangements in the network. Therefore, if temperature helps to delocalize stress as suggested by Fig. 6 and Fig. 7(f), this must be evident in the distribution of stress within the network. In Fig. 8(a-b) we show two snapshots of the same deformed network (p = 0.65) at two different temperatures, together with the associated histogram of the bond deformation8(c-d), which is equivalent to the probability distribution of the microscopic stresses since the springs are linear. Both networks are deformed up to 10% strain, which is close to ε p and well above the onset strain discussed in Sec. 3.2. At the lower temperature, the stress is distributed very heterogeneously, as indicated by an asymmetric distribution with an exponential tail containing few bonds carrying a high load. The regions of high stress, typically composed of aligned highly-stressed bonds, that we call force chains, can be readily identified in the simulation snapshot. On the contrary, for the higher temperature, the distribution is very symmetric, resembling a Gaussian distribution, and the force chains can not be identified. This indicates that even above the onset strain thermal fluctuations act as a stabilizing field as discussed in Sec. 3.1 and do affect the distribution of stress in the network. To quantify this heterogeneity, we calculate the excess kurtosis κ e of the stress distribution, an indicator of the tail heaviness of a distribution (being zero for a Gaussian). This measure has been recently used to quantify stress heterogeneities in porous materials [36]. To illustrate how the heterogeneity of the stress distributions is linked to the macroscopic stress evolution, in Fig. 8(c) we plot both κ e and σ as a function of strain for an example simulation run. As observed in most cases, the strain-stiffening of the network is accompanied by a similar increase in kurtosis. The stress distribution becomes more heterogeneous until (approximately) the first bond breaks (the dashed black line in Fig. 8(c)), after which strain softening occurs and the stress distribution becomes more homogeneous. This decrease in heterogeneity is presumably caused by redistribution of the stress after bonds are broken. Strikingly, in correspondence with the peak stress, a local minimum for the kurtosis is observed. The subsequent stress drops are instead accompanied by an increase in κ e , and therefore in the microscopic stress heterogeneity. This increase indicates stress concentration somewhere in the network leading to significant bond breakage that does not allow for larger stress response. To show how the stress heterogeneity changes with temperature, we plot the maximum and the minimum of κ e as a function of T * in Fig. 8(d). To determine the minimum kurtosis, we take the smallest value of κ e in a strain interval close to the peak strain ε p , to avoid lower values that might be found before strain-stiffening. Analogously, for the maximum kurtosis, we only look at the maximum up to and including the peak strain, to avoid post-peak values. Both quantities clearly decrease when the temperature is raised, but follow different curves. In particular, the maximum of κ e , that is associated to the network strain-stiffening, is immediately sensitive to temperature changes, in line with our previous observations on the non-linear elasticity (Sec. 3.2), while the minimum of κ e , associated to the fracture peak, exhibits an initial temperature insensitive interval, similarly to the other fracture descriptors investigated above. Furthermore, as the temperature increases, the difference between the maximum and minimum becomes smaller, and eventually both quantities reach zero (homogeneous stress distribution) at the melting temperature. Note that, while we have shown here results only for a given connectivity, the fact that a higher temperature allows for better redistribution of the stress during fracture was a consistent observation in our simulation study.

Summary and conclusions
In this work we explored the relation between rigidity and network failure under the influence of thermal fluctuations. Our results demonstrate that thermal fluctuations couple with network rigidity and affect the non-linear mechanics of elastic networks. In general, thermal fluctuations lead to a lower failure strength (Fig. 4), an increased ductility and increased fraction of broken bonds (Fig. 5). We have shown that at the microscopic level the failure response is altered with respect to the athermal case in two ways: i) bond failure can be activated by instantaneous thermal fluctuations, creating additional weak spots, and ii) stress is delocalized, suppressing the expansion of existing defects (Fig. 6). We reveal that temperature acts as a stabilizing field that resists large structural non-affine deformation within the network (Fig. 2 and Fig. 3). Specifically, the thermal fluctuations increase the lengthscale over which stress is redistributed, which can be quantified via the maximum stress drop (Fig. 7) and the excess kurtosis (Fig. 8). Although these trends can be observed for all connectivities, there are distinct damage mechanisms above and below the isostatic point. Above the isostatic point, the failure up to the peak stress shifts from crack nucleation at a single site to a more diffuse failure pattern, while around the isostatic point the failure response is already delocalized in the athermal limit and the fraction of broken bonds is enhanced approaching the melting point (Fig. 6). These distinct failure processes might explain the difference in how the peak strain depends on temperature with respect to rigidity (Fig. 4).
We note that at a first glance elastic networks subjected to thermal fluctuations behave like athermal networks in a stabilizing field. However, the instantaneous nature of the thermal fluctuations introduces important differences. It is striking that, without any applied deformation, the thermal fluctuations induce structural rearrangements of the average network structure (see ESI). Furthermore, providing enough time, the thermal fluctuations allow the failure of bonds even if they are not intrinsically under tension (activated failure), leading to diffuse damage. A final consequence of introducing thermal fluctuations is that time becomes an important parameter. In our simulations the system was deformed at a constant strain rate, i.e. it was driven at a given speed. If the driving speed is too low, the system will melt due to the process of activated failure. If the driving speed is too high, the system has no time for stress relaxation as it is held back by the viscous surroundings. Therefore, the failure response of an elastic network is generally determined by the coupling between the driving speed, viscosity, rigidity and thermal fluctuations. Our study was focused on a regime in which driving and viscosity effects were small (see ESI for discussion).
This work provides new insight into the relation between the static network structure, thermal fluctuations and the failure response of network materials. Above all, it shows that the ratio T /µ can be used to significantly change the failure response of a network. In relating this knowledge to soft matter systems, a clear challenge is to better understand the influence of the physics at the element level, such as plastic rearrangements in colloidal gel strands [37] and temperature sensitivity of elastic elements such as semi-flexible fibres [38].

Conflicts of interest
There are no conflicts to declare.
Electronic Supplementary Information for "The role of temperature in the rigidity-controlled fracture of elastic networks" Justin Tauber, Aimée R. Kok, Jasper van der Gucht, * and Simone Dussi Dependence of non-affinity on the reference coordinates We calculate the non-affine deformation with respect to the time averaged position at 0% strain ( Fig. S1(a)) r 0 . Alternatively, the non-affinity can be calculated with respect to the initial positions of the nodes (located on a regular triangular lattice) r init (Fig. S1(b)). The difference between Fig. S1(a) and Fig. S1(b) indicates that temperature has an effect on the equilibrium node positions at 0% strain. This is in contrast with athermal networks and bending stabilized networks, where the equilibrium node positions are equal to the initial position of the nodes. Fig. S1(c) shows that the average displacement of the nodes from their initial position to their equilibrium position dr = |r 0 − r init | depends on both temperature and connectivity. Especially below the isostatic point there are significant reorganizations within the network. However, at high temperatures also the network structure well above the isostatic point is affected. Clearly, the thermal fluctuations do affect the equilibrium structure at 0% strain.

Size of the thermal fluctuations
To quantify the size of the thermal fluctuations, we monitor the root mean squared displacement of the nodes with respect to their equilibrium position u 2 therm and define the size of the fluctuations as dr fluc = u 2 therm with · representing a time-average. From Fig. S2(a) it is clear that the size of the fluctuations of the nodes depends on both temperature and connectivity. In general the size of the fluctuations decreases with an increase in connectivity, indicating there is feedback between the number of constraints imposed on a node in the network and how far the nodes can move. Before measuring u 2 therm , all systems are subjected to the same calibration run of 100τ (see Fig. S2(b)/(c)). We note that for p = 0.65 ( Fig. S2(b)) the required time to reach a stable value of u 2 therm is longer for lower temperatures, indicating that the rate of equilibration depends on temperature. Furthermore, we see that at the higher connectivity value p = 0.90 (Fig. 2(c)) the fluctuations are smaller and reach their equilibrium value faster.

Relation between time-averaged non-affinity and instantaneous non-affinity
In an athermal elastic network the position of the crosslinks is determined by the applied deformation and the non-affine response of the nodes. In a thermal elastic network, the positions of the nodes are also influenced by thermal fluctuations of the nodes. The position of a node r under uniaxial extension ε can therefore be described as r(ε, T, p) = r 0 + u aff (ε) + u naff (ε, T, p) + u therm (ε, T, p) , where p is the network connectivity parameter, and T temperature. u stands for a displacement vector and r 0 is the time averaged position at 0% strain. If a system is fixed at a certain strain ε the average position of the particle over time will be, r = r 0 + u aff (ε) + u naff (ε, T, p) , assuming |u therm | = 0. If we instead monitor the average size of the fluctuations of the nodes we find that (r − r) 2 = u 2 therm (ε, T, p) Note that in case of drift in the system center of mass, this needs to be taken into account. In our case, we assume < r >= r com . Below, we will detail how these contributions are related to the measure for non-affinity Γ.