Thermally activated intermittent dynamics of creeping crack fronts along disordered interfaces

We present a subcritical fracture growth model, coupled with the elastic redistribution of the acting mechanical stress along rugous rupture fronts. We show the ability of this model to quantitatively reproduce the intermittent dynamics of cracks propagating along weak disordered interfaces. To this end, we assume that the fracture energy of such interfaces (in the sense of a critical energy release rate) follows a spatially correlated normal distribution. We compare various statistical features from the obtained fracture dynamics to that from cracks propagating in sintered polymethylmethacrylate (PMMA) interfaces. In previous works, it has been demonstrated that such an approach could reproduce the mean advance of fractures and their local front velocity distribution. Here, we go further by showing that the proposed model also quantitatively accounts for the complex self-affine scaling morphology of crack fronts and their temporal evolution, for the spatial and temporal correlations of the local velocity fields and for the avalanches size distribution of the intermittent growth dynamics. We thus provide new evidence that an Arrhenius-like subcritical growth is particularly suitable for the description of creeping cracks.


Propagation model
Constitutive equations. We consider rugous crack that are characterised by a varying and heterogeneous advancement a(x, t) along their front, x being the coordinate perpendicular to the average crack propagation direction, a the coordinate along it, and t being the time variable (see Fig. 1). At a given time, the velocity profile along the rugous front is modelled to be dictated by an Arrhenius-like growth, as proposed by Cochard et al. 8 : where V (x, t) = ∂a(x, t)/∂t is the local propagation velocity of the front at a given time and V 0 is a nominal velocity, related to the atomic collision frequency 41 , which is typically similar to the Rayleigh wave velocity of the medium in which the crack propagates 42 . The exponential term is a subcritical rupture probability (i.e., between 0 and 1). It is the probability for the rupture activation energy (i.e., the numerator term in the exponential) to be exceeded by the thermal bath energy k B T 0 , that is following a Boltzmann distribution 41 . The Boltzmann constant is denoted k B and the crack temperature is denoted T 0 and is modelled to be equal to a constant room temperature (typically, T 0 = 298 K). Using this constant temperature corresponds to the hypothesis that the crack is propagating slowly enough so that no significant thermal elevation occurs by Joule heating at its tip (i.e., as inferred by Refs. 39,40 ). Such propagation without significant heating is notably believed to take place in the experiments by Tallakstad et al. 28 that we here try to numerically reproduce, and whose geometry is shown in Fig. 1. Indeed, www.nature.com/scientificreports/ their reported local propagation velocities V did not exceed a few millimetres per second, whereas a significant heating in acrylic glass is only believed to arise for fractures faster than a few centimetres per second 40,44 . See the supplementary information for further discussion on the temperature elevation. In Eq. (1), the rupture activation energy is proportional to the difference between an intrinsic material surface fracture Energy G c (in J m −2 ) and the energy release rate G at which the crack is mechanically loaded, which corresponds to the amount of energy that the crack dissipates to progress by a given fracture area. As the front growth is considered subcritical, we have G < G c . We here model the fracture energy G c to hold some quenched disorder that is the root cause for any propagating crack front to be rugous and to display an intermittent avalanche dynamics. This disorder is hence dependent on two position variables along the rupture interface. For instance, at a given front advancement a(x, t), one gets G c = G c (x, a) . The coefficient α 2 , in Eq. (1), is an area which relates the macroscopic G and G c values to, respectively, the microscopic elastic energy U = α 2 G stored in the molecular bonds about to break, and to the critical energy U c = α 2 G c above which they actually break. See Vanel et al. 36 , Vincent-Dospital et al. 40 or the supplementary information for more insight on the α 2 parameter, which is an area in the order of d 3 0 /l , where d 0 is the typical intra-molecular distance and l is the core length scale limiting the stress divergence at the crack tip.
Finally, the average mechanical load that is applied on the crack at a given time is redistributed along the evolving rugous front, so that G = G(x, t) . To model such a redistribution, we here use the Gao and Rice 3 formalism, which integrates the elastostatic kernel along the front: In this equation, G is the mean energy release rate along the front and PV stands for the integral principal value. We, in addition, considered the crack front as spatially periodic, which is convenient to numerically implement a spectral version of Eq. (2) 45 as explained by Cochard et al. 8 .
Equations (1) and (2) thus define a system of differential equations for the crack advancement a, which we have then solved with an adaptive time step Runge-Kutta algorithm 46 , as implemented by Hairer et al. 47 . The complete code for the crack simulation is available as a Software Heritage archive 48 . Further details on the code can be obtained by contacting the authors.

Discretization.
In this section, we discuss the main principles we have used in choosing the numerical accuracy of our solver. The related parameters are illustrated in Fig. 2.
In attempting to correctly reproduce the experimental results of Tallakstad et al. 28 , this solver needs to use space and time steps, here denoted x s and t s , at least smaller than those on which the experimental fronts were observed and analysed. Thus, x s needs to be smaller than the experimental resolutions in space (the camera pixel size) x = a of about 2 to 5 µ m and 1/�t s needs to be higher than the experimental camera frame rate 1/�t . This frame rate was set by Tallakstad et al. 28 to about (100V )/�x , where V is the average front velocity of a given fracture realisation. The propagation statistics of our simulated fronts, henceforward shown in this manuscript, have, for consistency, always been computed on scales comparable to the experimental x , a , t steps. Thus, as �x s < �x and �t s < �t , we have first decimated the dense numerical outputs on the experimental observation grid, by discarding smaller time scales and by averaging smaller space scales to simulate the camera frame rate and pixel size.
As the camera resolution was 1024 pixels, the lengths L of the crack segments that Tallakstad et al. 28 analysed were 1024 x = 3 to 7 mm long, and we have then analysed our numerical simulations on similar front widths. Yet, these simulations were priorly run on longer front segments, L s > L , in order to avoid any possible edge effects in the simulated crack dynamics (for instance in the case where L would not be much bigger than the typical size of the G c quenched disorder).

Figure 2.
Illustration of the discretization principles and of the solver and observation grids. Three crack fronts at three successive times are shown, over which the parameters discussed in the "Discretization" section are defined. www.nature.com/scientificreports/ Overall, we have checked that the numerical results presented henceforward were obtained using a high enough time and space accuracy for them to be independent of the discretization (as is shown in the supplementary information).

Physical parameters values.
For the model dynamics to be compared to the experiments 28 , one must also ensure that the V 0 , α , T 0 , G and G c parameters are in likely orders of magnitude.
As V 0 is to be comparable to the Rayleigh velocity of acrylic glass, we have here used 1 km s −149 . Lengliné et al. 37 furthermore estimated the ratio α 2 /(k B T 0 ) to be about 0.15 m 2 J −1 and they could approximate the where G c is the average value of G c in the rupturing interface. With our choice on the value of V 0 , we then deduce G c ∼ 250 J m −2 (note that the trade-off between V 0 and G c should be kept in mind when comparing our results with those by Cochard et al. 8 , as both papers use a different V 0 ). The value thus inverted for the fracture energy ( 250 J m −2 ), that is to represent the sintered PMMA interfaces, is logically smaller but comparable to that inferred by Vincent-Dospital et al. 40 for the rupture of bulk PMMA (about 1300 J m −2 ). Qualitatively, the longer the sintering time, the closer one should get from such a bulk toughness, but the less likely an interfacial rupture will occur. Experiments in two different regimes were run 28 : a forced one where the deflection of the lower plate (see Fig. 1) was driven at a constant speed, and a relaxation regime, where the deflection was maintained constant while the crack still advances. In both scenarii, the long term evolution of the average load G(t) and front position a(t) was shown 8,37 to be reproduced by Eq. (1). In the case of the experiments of Tallakstad et al. 28 , the intermittent dynamics measured in the two loading regimes were virtually identical. Such similarity likely arises from the fact that the avearge load G was, in both cases, computed to be almost constant over time, in regard to the spatial variation in G, described by Eq. (2) (see the supplementary information). Here, we will then consider that the crack is, in average along the front, always loaded with the same intensity (i.e., G(t) = G).
The actual value of G , together with the average surface fracture energy of the medium G c , then mainly controls the average crack velocity V . This average velocity was investigated over five orders of magnitude in Ref. 28 , from 0.03 to 140 µm s −1 , which, in our formalism, shall correspond to values of (G c − G ) between 145 and 85 J m −2 , respectively, which is actually consistent with the values of G measured by Lengliné et al. 37 for cracks propagating at similar speeds. The intermittency of the crack motion was experimentally inferred to be independent on V and we show, in the supplementary information, that it is also the case in our simulations. The velocity variation along the front shall then only arise from the disorder in G c and from the related variations of G due to the roughness of the crack front. Further in this manuscript, we will use G = 120 J m −2 , which corresponds to an average propagation velocity of about 1.5 µm s −1 .

Heterogeneous fracture energy
Of course, the actual surface fracture energy field in which the rupture takes place will significantly impact the avalanches dynamics and the crack morphology. Such a field is yet a notable unknown in the experimental set-up of Tallakstad et al. 28 , as their interface final strength derived from complex sand blasting and sintering processes. Although these processes were well controlled, so that the rough rupture experiments were repeatable, and although the surfaces prior to sintering could be imaged 43 , the actual resulting distribution in the interface cohesion was not directly measurable. While this is, to some extent, a difficulty in assessing the validity of the model we present, we will here show that a simple statistical definition of G c is enough to simulate most of the avalanches statistics.
We will indeed consider a normally distributed G c field around the mean value G c with a standard deviation δG c and a correlation length l c . Such a landscape in G c is shown in Fig. 3a, and we proceed to discuss the chosen values of δG c and l c in the "Growth exponent and fracture energy correlation length" and "Local velocity distribution and fracture energy standard deviation" sections.
Growth exponent and fracture energy correlation length. Among the various statistical features studied by Tallakstad et al. 28 , was notably quantified the temporal evolution of their fracture fronts morphology. It was interestingly inferred that the standard deviation of the width evolution of a crack front h scales with the crack mean advancement: In this equation, x is a given position along the front, t is a time delay from a given reference time t 0 , and h writes as a being the average crack advancement at a given time. To mitigate the effect of the limited resolution of the experiments and obtain a better characterization of the scaling of the interfacial fluctuations on the shorter times, we computed the subtracted width, as proposed in Barabasi and Stanley 50 , and done by Tallakstad et al. 28 (whose experiments we here reproduce) and Soriano et al. 51. The scaling exponent β G is referred to as the growth exponent, and we will here show how it allows to deduce a typical correlation length for the interface disorder. Indeed, β G was measured to be 0.55 ± 0.08 by Tallakstad  28 . This value is close to 1/2, that is, consistent with an uncorrelated growth process (e.g., 50 ), such as simple diffusion or Brownian motion. We thus get a first indication on the disorder correlation length scale l c . To display an uncorrelated growth when observed with the experimental resolution ( x ∼ 3 µm), the fronts likely encountered asperities whose size was somewhat comparable to this resolution. Indeed, if these asperities in G c were much bigger, the growth would be perceived as correlated. By contrast, if they were much smaller (orders of magnitude smaller), the rugosity of the front would not be measurable, as only the average G c over an observation pixel would then be felt. Furthermore, and as shown in Fig. 4a, the exponent β G was observed on scales (Vt) up to 100 µ m, above which W stabilised to a plateau value of about 30 µ m. A common picture is here drawn, as both this plateau value and the typical crack propagation distance at which it is reached are likely to also be correlated with l c , as the front is to get pinned on the strongest asperities at this scale. From all these clues, we have considered, in our simulations, the correlation length of the disorder to be about l c = 50 µ m, and we show in Fig. 4a that it allows an approximate reproduction of the front growth exponent and of the plateau at high Vt . Note that the accuracy reported for the exponents in this manuscript is estimated by fitting various portions of the almost linear data points and reporting the dispersion of the thus inverted slopes. In Fig. 4b, we also show how varying l c impacts W, and, in practice, we have chosen l c by tuning it when comparing these curves to the experimental one. Noteworthily, the thus chosen l c is in the lower range of the size of the blasting sand grains ( 50 − 300 µ m) that were used 28 to generate the interface disorder. It is also comparable to the correlation length of the blasting induced topographic anomalies ∼ 18 µ m on the post-blasting/pre-sintering PMMA surfaces, as measured by Santucci et al. 43

by white light interferometry.
Local velocity distribution and fracture energy standard deviation. While the crack advances at an average velocity V , the local velocities along the front, described by Eq. (1), are, naturally, highly dependent on the material disorder: the more diverse the met values of G c the more distributed shall these velocities be.
Maløy et al. 26 and Tallakstad et al. 28 inferred the local velocities of their cracks with the use of a so-called waiting time matrix. That is, they counted the number of discrete time steps a crack would stay on a given camera pixel before advancing to the next one. They then deduced an average velocity for this pixel by inverting this number and multiplying it by the ratio between the pixel size and the time between two pictures: �a/�t . Such a method, that provides a spatial map V(x, a), was applied to our simulated fronts, and we show this V(x, a) map in Fig. 3c. As to any time t corresponds a front advancement a(x, t) (recorded with a resolution a ), an equivalent space-time map V(x, t) can also be computed, and it is shown in Fig. 3b. The experimental report 28 presented the probability density function of this latter (space-time) map P(V /V ) , and it was inferred that, for high values of V, the velocity distribution scaled with a particular exponent η = 2.6 ± 0.15 28,38 (see Fig. 5a). That is, it was observed that  28 , using the waiting time matrix (see text for details). The velocity are plotted related to the average crack velocity V = 1.5 µm s −1 . All parameters used to run the corresponding simulation are summarised in Table 1 8 , who introduced the model that we here discuss, inferred that the η exponent was mainly depend- . Truly, a more comprehensive expression could also include other quantities, such as V 0 or l c . Yet, as all other parameters have now been estimated, we can deduce δG c by varying it to obtain η ∼ 2.6 . We show, in Fig. 5b, how varying δG c impacts P(V /V ) and η . We found δG c ∼ 35 J m −2 . In Fig. 5a, we show the corresponding velocity distribution for a simulation run with this parameter, together with that from Tallakstad et al. 28 , showing a good match. Note that the ability of the model to reproduce the local velocity  Table 1. The slope and plateau of the experimental data (shown in (a)) is marked by the dashed line for comparison. We chose the value of δG c by tuning it and fitting the slope and maximum of the experimental data, which are illustrated by the dashed lines. The rest of the parameters used in these simulations are as defined in Table 1.
Note that the ability of the model to reproduce the local velocity distribution was already shown by Cochard et al. 8 (see text for explanation and discussion).  8 , and this figure mainly aims at illustrating our calibration of the fracture energy field. The model we present is also slightly different to that of Cochard et al., as the interface fracture energy is, contrarily to this previous study, now described at scales below its correlation length, similarly to the observation scale of the experiments. We here verify that the reproduction of the local velocity distribution is still valid at these small scales. Satisfyingly, the inverted value of δG c is not too far from the value found by Lengliné et al. 37 for their fluctuation in the mean fracture energy G c along their sintered plates, when studying the mean front advancement (i.e., neglecting the crack rugosity) in similar PMMA interfaces, which was about 25 J m −2 .

Further statistics
We have now estimated the orders of magnitude of all parameters in Eqs. (1) and (2), including a likely distribution for an interface fracture energy representative of the experiments 28 we aim to simulate (i.e., including G c , δG c and l c ). For convenience, this information is summarised in Table 1.
We will now pursue by computing additional statistics of the crack dynamics to compare them to those reported by Tallakstad et al. 28 .

Local velocities correlations.
In particular, we here compute the space and time correlations of the velocities along the front. That is, four correlation functions that are calculated from the V(x, t) and V(x, a) matrices (shown in Fig. 3) and defined as: where i and j are the variables of either V(x, t) or V(x, a) and δi a given i increment. V j is the mean of V(i, j) taken along j at a given i 0 . The corresponding δV j is the velocity standard deviation along the same direction and for the same i 0 . The correlation functions hence defined are the same as those used by Tallakstad et al. 28 on their own data, allowing to display a direct comparison of them in Fig. 6. A good general match is obtained.
One can notice the comparable cut-offs along the x axis (Fig. 6a,c), indicating that our chosen correlation length for the interface disorder ( l c inferred in the "Growth exponent and fracture energy correlation length" section) is a good account of the experiment. Yet, one can notice that C xt (the velocity correlation along the crack front shown in Fig. 6a) is higher in the numerical case than in the experimental one. It could translate the fact that the experimental disorder holds wavelengths that are smaller than the observation scale x , and that our modelled G c distribution, where l c > �x , is rather simplified.
To go further, Tallakstad et al. 28 modelled C xt as and inverted the values of τ x and x * to, respectively, 0.53 and about 100 µ m. Doing a similar fit on the simulated data, we found τ x ∼ 0.13 and x * ∼ 94 µ m. The related function is displayed in Fig. 6a (plain line). Our small τ x ∼ 0.13 may derive, as discussed, from the better correlation that our simulation displays at small δx ( τ x may in reality tend to zero for scales smaller than those we observe) compared to that of the experiments, while the matching x * probably relates to a satisfying choice we made for l c . Overall, the existence of a clear scaling law at  (1) and (2) (8), is rather uncertain (see the two experimental plots in Fig. 6a) so that mainly the cut-off scale is of interest. On the time correlation C tx (Fig. 6b), one can similarly define the parameters A, τ t and a * to fit Eq. (7) with a function where A is a constant of proportionality. Fitting this function to Eq. (7) with a least-squares method, we found τ t ∼ 0.3 and a * ∼ 4.3 µ m, and this fit is represented by the dashed line in Fig. 6b. Tallakstad et al. 28 reported τ t ∼ 0.43 and a * ∼ 7 µ m. Figure 6b shows the experimental and simulated correlation functions in the V δt/a * -C tx (a * /V ) τ t /A domain, as this allowed a good collapse of the data from numerous experiments 28 . We show that it also allows an approximate collapse of our modelled correlation on a same trend. Finally, the derived value of a * consistently matches the apparent cut-off length in the C a x correlation function in Fig. 6d. This length being of a magnitude similar to that of the observation scale a , the crack local velocities appear uncorrelated along the direction of propagation, which is consistent with the β G ∼ 1/2 growth exponent (e.g., 50 ). Avalanches size and shape. We pursue by characterising the intermittent, burst-like, dynamics of our crack fronts and, more specifically, the avalanche (or depinning) and pinning clusters shown by the local front . The plain points were computed from the simulation whose parameters are presented in Table 1 and the hollow stars are some of the experimental data points extracted from Figs. 6 and 7 of Tallakstad et al. 28 . In plot (a), the line overlying the numerical data set corresponds to a fit using Eq. (8). (b) is plotted in a domain that allowed a good collapse of the experimental data for many experiments 28 . The parameters A, a * and τ t were inverted from Eq. (9) ). We define an avalanche when the front velocity locally exceeds the mean velocity V by an arbitrary threshold that we denote c, that is, when Similarly, we state that a front is pinned when We then map, in Fig. 7, the thus defined avalanching and pinned locations of the crack. Following the analysis of Tallakstad et al. 28 , we compute for each of these clusters the surface S, the crossline extent l x (that is, the maximum of a cluster width in the x direction) and the inline extent l a . The definition chosen for l a varies for the avalanche clusters, where the maximal extent along the a direction is regarded, or the pinned one, where the mean extent along the a direction is rather used. This choice was made 28 because the pinning clusters tend to be more tortuous so that their maximum span along the crack direction of propagation is not really representative of their actual extent (see Fig. 7). In Fig. 8a, we show the probability density function of the cluster surface P(S) and compare it to the experimental one. One can notice that it behaves as with γ = 1.44 ± 0.15 . This value is comparable to the exponent inverted experimentally 28 , that is, γ = 1.56 ± 0.04 .
Of course, the size of the avalanche (depinning) clusters highly depends on the chosen threshold c, but we verified, as experimentally reported, that the value of γ inverted from the simulated data is not dependent on c, as shown in Fig. 8b. We also show, in     Fig. 3c, as per Eqs. (10) and (11). Two thresholds are here used to define these maps relatively to the mean velocity: c = 3 and c = 6 . The white areas are the locations of interest, of surfaces S, crossline extents l x and inline extents l a . Bottom images: Difference in definition of l a for the avalanche (or depinning) and pinning clusters shown in Fig. 7. For the former, l a is the maximum extent along the a direction. For the latter it is the average width in the same direction. In both cases, l x is the maximum extent along the x direction and S the full surface (in white) of the cluster. The square pattern marks the pixel size ( x = a = 3 µm). We also computed the probability density function of l x and l a , that are respectively compared to their experimental equivalent in Fig. 9. These functions can be fitted with The modelled S is expressed in pixels (one pixel is 9 µm 2 ) and the experimental S reported by Tallakstad et al. 28 (in their Fig. 13) is in an arbitrary unit, so that the magnitude of both should not here be compared. We have here shifted up this experimental data set for an easier comparison with the numerical one. The straight line has a slope 0.68. Figure 9. (a) Probability density function of the crossline extent l x of the modelled avalanche clusters (plain points) and of the modelled pinning clusters (crosses), for a threshold c = 3 . The straight line has a slope β x = 1.7 , as per Eq. (13). The hollow stars shows the experimental probability density function obtained by Tallakstad et al. 28 for the pinning and avalanche clusters (from their Fig. 16a, inset, c = 3). (b) Probability density function of the inline extent l a of the modelled avalanche clusters (plain points) and of the modelled pinning clusters (crosses), for a threshold c = 3 . The two straight dashed lines have a slope β x = 2.2 , inline with that of the experimental data from Tallakstad et al. 28 for the pinning clusters (hollow stars, from their Fig. 16b, inset, c = 3). It should be noted that, while we have here fitted P(S), P(l x ) and P(l a ) with plain scaling laws (i.e., with Eqs. (12) to (14)), Tallakstad et al. 28 also studied the cut-off scales above which these scaling laws vanish in the experimental data, and the dependence of these cut-off scales with the arbitrary threshold c. In our case, such scales are challenging to define, as one can for instance notice in Figs. 8 and 9, where an exponential cut-off is not obvious. This may result from a limited statistical description of the larger avalanches in our simulations. Similar cut-off scales, decreasing with increasing c should however hold in our numerical data, in order to explain the decrease of average avalanche size with c, as shown in Fig. 8c. Front morphology. Finally, we show, in Fig. 10a, the relations between the clusters surface S and their linear extent l x and l a . Here, l x and l a are the mean extents for all the observed clusters sharing a same surface (with the given pixel size limiting the resolution). We could fit these relations with l x ∝ S 0.77 and l a ∝ S 0.25 for the pinning clusters, and with l x ∝ S 0.64 and l a ∝ S 0.47 for the avalanches clusters. It is in qualitative agreement with the laws observed by Tallakstad et al. 28 : l x ∝ S 0.63 and l a ∝ S 0.34 for the pinning clusters, and S ∝ l x 0.61 and l a ∝ S 0.41 for the avalanches clusters. These exponents were experimentally reported with a ±0.05 accuracy, and we estimated comparable error bars for the numerically derived ones. Thus, the shape of our simulated avalanches and pinned locations is rather similar to the observed experimental ones. Note that, from all the previous exponents, one can easily define H such that l a ∝ l x H , and we thus have H p ∼ 0.25/0.77 = 0.32 ± 0.1 and H d ∼ 0.47/0.64 = 0.73 ± 0.01 for, respectively, the simulated pinning and depinning clusters (see Fig. 10b).
It was suggested 5,52 that H is a good indicator of the front morphology, as the front shape is to be highly dependent on the aspect ratio of its avalanches. To verify this hypothesis, we computed the advancement fluctuation along the front σ , that is While this quantity was not presented by Tallakstad et al. 28 , it was provided by other experimental works done on the same set-up 26,27 , and Fig. 11a shows σ as reported by these authors, together with that computed in the output of our simulation. One can notice that the numerical fronts are less rugous than the experimental ones. Such a mismatch is here due to the fact that the experiment from Santucci et al., shown in Fig. 11a, had more rugous crack fronts than the one from Tallakstad et al., to which the simulation is calibrated (as shown in Fig. 4). In both cases, the data sets seem to present two self-affine behaviours (e.g., 50 ) with a Hurst exponent ζ that differs at low and high length scales. Noting δx * the cut-off between these length scales we indeed have: www.nature.com/scientificreports/ We derived ζ − ∼ 0.68 ± 0.05 and ζ + = 0.4 ± 0.05 for the simulation, which compare well to the exponents that were measured experimentally, respectively, ζ − = 0.60 ± 0.05 and ζ + = 0.35 ± 0.05 and which are also close to the values we found for H d and H p . The cut-off scale between the two regimes is also similar in both the experimental and numerical cases: δx * ∼ 80 µ m, comparable to the disorder correlation length l c , and to the length scales x * , below which the local propagation velocities are correlated. For scales above this correlation length, Cochard et al. 8 showed, by analytically analysing the same model as we here study, that the front morphology is dominated by the material quenched disorder with a Hurst coefficient approximating to ζ + = 0.5 . At even larger scales, above R c ∼ πl c α 2 G c /(k B T 0 ) , they also showed 8 that the roughness of the simulated cracks ceases to be governed by the quenched disorder but is rather dominated by the thermal (annealed) noise, with σ decaying logarithmically and with a Hurst coefficient tending to ζ ∞ = 0 . With our set of parameters, R c computes to 6 mm, which is close to, yet bigger than, the total analysed length of the front. The value ζ + ∼ 0.4 , that we have here inverted, arises then likely from the transition between the two regimes, ζ + = 0.5 and ζ ∞ = 0 , as already mentioned for the experimental case, in Ref. 27 . In addition to a theoretical Hurst exponent ζ + = 0.5 , Cochard et al. 8 computed an analytical approximation for the fronts morphology power spectrum P a (Ŵ) , for the length scales Ŵ for which the effect of the quenched disorder prevails: We show, in Fig. 11b, how this approximation also fits the power spectra of our modelled front.

Discussion and conclusion
We studied an interfacial fracture propagation model, based only on statistical and subcritical physics in the sense of an Arrhenius law (Eq. (1)) and on the elastic redistribution of stress along crack fronts (Eq. (2)). Following the work of Cochard et al. 8 , we here showed that it allows a good representation of the intermittent dynamics of fracture in disordered media, as it approximately mimics the scaling laws dictating the propagation of experimental fronts, such as their growth exponent, their local velocity distribution and space and time correlations, the size of their avalanches and their self-affine characteristics.
To run our simulations, we had to assume a given distribution for the toughness of the rupturing interface, as this quantity is not directly measurable in the laboratory. We proposed G c to be normally distributed with a unique correlation length and, of course, this can only be a rough approximation of the actual fracture energy obtained by Tallakstad et al. 28 by sintering two sand-blasted plexiglass plates. From this approximation, could arise discrepancies between our simulations and the experiments. We have indeed shown how some of the observed exponents were strongly dependent on the definition of the material disorder. We also have assumed a perfectly elastic crack front, when the local dynamics of creeping PMMA could be visco-elastic in  15), for the simulation (plain points) and an experimental data set from Santucci et al. 27 (see their Fig. 4). Different self-affine behaviours are observed above and below the δx * cut-off, with comparable Hurst exponents ζ . The dashed lines mark the slopes fitted on the simulation data for the two cases. The experimental points are from an experiment different from those of Tallakstad et al. 28 to which the model was calibrated. (b) Power spectra of the simulated crack advancement, averaged over 10,000 consecutive fronts. It is shown both before (raw) and after (binned) binning the fronts to the experimental camera pixel size. The difference between these two plots shows an influence of the observation scale on the small-scale study of the crack morphology. The plain line is the approximation 8 from Eq. (18), which is valid between l c and R c , where the morphology is dominated by the material quenched disorder. Note that the scaling regime for scales above l c was already studied by Cochard et al. 8 , while the model match to the experiment below this cut-off scale, shown in (a), is a new result, as already discussed in the "Local velocity distribution and fracture energy standard deviation" section. www.nature.com/scientificreports/ part, particularly below the typical length scale r ∼ GE/σ 2 y ∼ 30 µ m for plasticity around crack tips (e.g., 1 ) in this material, where σ y ∼ 100 MPa is the tensile yield stress of the polymer and E ∼ 3 GPa its Young modulus 53 .
These points being stated, the vast majority of the statistical quantities that we have here studied show a good match to those from the experimental observations, so that both the considered physical model and the interface definition are likely to be relevant. A further validation of this thermally activated model could derive from the comparison of its predictions with interfacial experiments at various background temperatures T 0 . However, such experimental data is, to our knowledge, not yet available. Of course, some of our considered parameters (e.g., G c , V 0 or α ) may, in practice, be temperature dependent so that a straight transposition of the model to different background temperatures could prove to be too simple. Creep experiments in bulk PMMA at various room temperatures can however be found in the literature 54 , where only the mean front velocity versus the mean mechanical load are measured. In this case 54 , it is reported that the creep dynamics is compatible with an Arrhenius-like process. By submitting many different materials to a constant load, at various temperatures, their lifetime was also shown 33,36 to follow an Arrhenius law, with an energy barrier that decreases with the applied stress. These materials include metals, alloys, non-metallic crystals and polymers (and PMMA in particular).
It should be noted that, as stated in our introduction, other models have been considered to numerically reproduce the interfacial PMMA experiments, notably, a non-subcritical threshold based fluctuating line model by Tanguy et al. 29 , Bonamy et al. 4 or Laurson et al. 5,20 and a fiber bundle approach by Schmittbuhl et al. 6 , Gjerden et al. 30 or Stormo et al. 31 . The present manuscript does not challenge these other models per se, but rather offers an alternative explanation to the intermittent propagation of rough cracks. The former model, the fluctuating line model 4,5,20,29 , considers a similar redistribution of energy release rate G as proposed in Eq. (2), but with a dynamics that is thresholded rather than following a subcritical growth law. The fronts either move forward by one pixel 5 if G > G c , or with a velocity proportional 4 to ( G − G c ). It is completely pinned otherwise ( V = 0 for G < G c ). While reproducing several statistical features of the experiments, this non-subcritical line propagation model does not simulate the mean propagation of cracks in various loading regimes (as done by Cochard et al. 8 ) or the distribution in local velocity 55 , and, in particular, the power law tail of this distribution (i.e., Fig. 5). By contrast, the latter model 6,30,31 , the fiber bundle one, can reproduce this particular power law tail. It is not a line model: the interface is sampled with parallel elastic fibers breaking at a given force threshold. This threshold is less in the vicinity of the crack than away from it (it is modelled with a linear gradient), explaining why the rupture is concentrated around a defined front, and it holds a random component in order to model the quenched disorder of the interface. An advantage of the fiber bundle model is to be able to describe a coalescence of damage in front of the crack 56 rather than solely describing a unique front. This could likely also be achieved in a subcritical framework, but would require to authorise damage in a full 2D plane, or require a full 3D modelling (i.e., also authorise out-of-plane damage), rather than only the modelling of a 2D front. In practice, thermal activation and damage coalescence may occur simultaneously. The observation of actual damage nucleation, in the experiments that we reproduce, has however never been obvious. Instead, the experimental fronts look rather continuous. Coalescence could yet still be at play at length scales smaller than the observation resolution. This being stated, an advantage of our model is to only rely on the experimental observations, on stress redistribution and on statistical physics. Another clear advantage of the Arrhenius based model, when compared to the other ones, is to hold a subcritical description that is physically well understood and that is a good descriptor of creep in many materials 1,36 . For the record, we show in Table 2 a comparison between the different exponents predicted by the three models, that all successfully reproduce some experimental observables.
Note that, if linearizing Eq. (1) with a Taylor expansion around G c − G , that is, for propagation velocities close to the mean crack speed V = V 0 exp(−α 2 [G c − G]/[k B T 0 ]) , one obtains  − G), 0] and where the coefficient of proportionality µ was named the 'effective mobility of the crack front' . Equation (19) may give some insight in the physical meaning of µ in this alternative model 4 . While the above similitude in mathematical forms may explain the obtention of some similar exponents in the dynamics of the two models (see Table 2), Eq. (19) is only a crude approximation of our highly non-linear Arrhenius formalism, which, as discussed below, allows a more exhaustive description of the experimental intermittent creep dynamics. In our simulations, the exponential Arrhenius probability term, describing the crack velocity, ranges over more than three orders of magnitude while ( G c − G ) ranges over less than two decades. Continuing with the comparison of our model with pre-existing ones, we had, in our case, to calibrate the disorder to the experimental data, in particular to accurately reproduce the η exponent, that is, to reproduce the fat tail of the crack velocity distribution. Paradoxically, this exponent, which is not accounted for by the other line model, has been found to be rather constant across different experiments and experimental set-ups. It could indicate that, in practice, the disorder obtained experimentally from the blasting and sintering of PMMA plates has always been relatively similar. Such qualitative statement is of course difficult to verify, because there exists no direct way of measuring the fracture energy of the experimental sintered samples. From Fig. 5b, one can yet notice that the calibration of the disorder amplitude does not need to be particularly accurate to obtain a qualitative fit to the experimental velocity distribution. The spread of the η exponent, for large disorders, is not that important in our model for the range of considered δG c , which can also be seen in Fig. 5 of Cochard et al. 8 . Gjerden et al. 57 suggested that the nucleation of damages, predicted by their fiber bundle model, led to a new -percolation -universality class for the propagation of cracks, explaining in particular the robustness of the exponent η . Their studies are however also numerical and cover a finite range of disorders, and an extra analytical proof would be needed to show that a system of infinite size would lead exactly to the same exponent, for any disorder distribution shape and amplitude.
Despite the variety in models reproducing the rough dynamics of creep, the present work provides additional indications that a thermodynamics framework in the sense of a thermally activated subcritical crack growth is well suited for the description of creeping cracks. Such a framework has long been considered (e.g., [32][33][34]36,58,59 ), and, additionally to the scaling laws that we have here presented, the proposed model was proven to fit many other observable features of the physics of rupture 8,37,39,40 . It accurately recreates the mean advancement of cracks under various loading conditions 8,37 , including when a front creeps in a spontaneous (not forced) relaxation regime, which cannot be achieved with the other (non subcritical) models, predicting an immobile front. When coupled with heat dissipation at the fracture tip, our description also accounts for the brittleness of matter 40 and for its brittle-ductile transition 39 .
Indeed, for zero dimensional (scalar) crack fronts, it was shown 40 that the thermal fluctuation at the crack tip, expressed as a deviation of the temperature from T 0 in Eq. (1), can explain the transition between creep and abrupt rupture, that is, the transition to a propagation velocity close to a mechanical wave speed V 0 , five orders of magnitude higher than the maximal creep velocity V that was here modelled. It was also shown, similarly to many phase transition problems, that such a thermal transition could be favoured by material disorder 39 . Thus, a direct continuation of the present work could be to introduce such a heat dissipation for interfacial cracks in order to study how brittle avalanches nucleate at given positions (typically positions with weaker G c ) to then expand laterally to become bulk threatening events.