Scale-free bursting activity in shrinkage induced cracking

Based on computer simulations of a realistic discrete element model we demonstrate that shrinkage induced cracking of thin layers of heterogeneous materials, generating spectacular crack patterns, proceeds in bursts. These crackling pulses are characterized by scale free distributions of size and duration, however, with non-universal exponents depending on the system size and shrinking rate. On the contrary, local avalanches composed of micro-cracking events with temporal and spatial correlation are found to obey a universal power law statistics. Most notably, we demonstrate that the observed non-universality of the integrated signal is the consequence of the temporal superposition of the underlying local avalanches, which pop up in an uncorrelated way in homogeneous systems. Our results provide an explanation of recent acoustic emission measurements on drying induced shrinkage cracking and may have relevance for the acoustic monitoring of the electro-mechanical degradation of battery electrodes.


Scale-free bursting activity in shrinkage induced cracking
Roland Szatmári 1 , Akio Nakahara 2 , So Kitsunezaki 3 & Ferenc Kun 1,4* Based on computer simulations of a realistic discrete element model we demonstrate that shrinkage induced cracking of thin layers of heterogeneous materials, generating spectacular crack patterns, proceeds in bursts.These crackling pulses are characterized by scale free distributions of size and duration, however, with non-universal exponents depending on the system size and shrinking rate.On the contrary, local avalanches composed of micro-cracking events with temporal and spatial correlation are found to obey a universal power law statistics.Most notably, we demonstrate that the observed non-universality of the integrated signal is the consequence of the temporal superposition of the underlying local avalanches, which pop up in an uncorrelated way in homogeneous systems.Our results provide an explanation of recent acoustic emission measurements on drying induced shrinkage cracking and may have relevance for the acoustic monitoring of the electro-mechanical degradation of battery electrodes.
Under steady external driving a broad class of disordered systems exhibits an intermittent response evolving through a bursting sequence of the nucleation and growth of cracks 1 , collapse of voids 2,3 , motion of magnetic domain walls 4 or dislocations 5 , rearrangement of particles 6 , firing of neurons 7 , etc. Experimental and theoretical investigations have revealed a scale-free statistics of crackling noise with a high degree of robustness across a broad range of length scales [8][9][10] .As a specific type of slowly driven systems, shrinkage induced cracking is ubiquitous in nature typically occurring when a material layer gradually dries or cools while attached to a substrate.The slowly evolving stress field leads to spectacular crack patters which can be observed from drying lake beds 11 , through the polygonal ground patterns in permafrost regions of Earth and Mars 12 , to the formation of columnar joints in cooling lava 11,[13][14][15] .For industrial applications shrinkage induced cracking, on the hand, has a high potential making possible the structural control of crack patterns 16,17 , on the other hand, however, it may be undesired since it can lead to quality reduction of products 18 .In particular, restrained shrinking of solidifying dense pastes such as concrete is accompanied by the accumulation of cracks which may lead to strength reduction, with additional long lasting consequences due to the penetration of corrosive agents into construction components 19,20 .Recently, the application of the acoustic emission technique revealed that the crackling noise generated in drying concrete has scale-free statistics with a changing power law exponent as shrinking proceeds 19,20 .A similar cracking mechanism occurring in the electrode material has been found to be responsible for the capacity fading of e.g.Li-ion batteries during repeated charging-discharging cycles 21 .The acoustic monitoring of the bursts of electrode cracking proved to be a promising non-destructive technique to follow the electro-mechanical degradation process 22,23 .Experimental and theoretical studies of the past decades focused mainly on the structure of the emerging crack pattern 11,13,24,25 , however, understanding the emergence of shrinkage induced crackling noise, and its statistical and dynamical properties still remained an open problem of fundamental importance.
Here we study this problem by means of a realistic discrete element model (DEM) of a thin material layer which undergoes restrained shrinking while attached to a rigid substrate.Based on computer simulations we show that cracking of the layer proceeds in bursts whose magnitude and duration obey scale-free statistics in agreement with the experimental findings.However, we demonstrate that power law exponents are not universal, in particular, they depend on the sample size and on the driving rate.Our calculations revealed that the nonuniversal behaviour of crackling noise exponents can be understood in terms of the temporal superposition of local avalanches of growing cracks, which are characterized by universal distributions.Our results provide a possible explanation of the experimental findings with important implications for testing and assessing of www.nature.com/scientificreports/construction materials whose manufacturing involves solidification accompanied by shrinking, and for the acoustic monitoring of battery degradation.

Results
We use a discrete element model (DEM) to study the cracking process of a thin material layer while it shrinks attached to a rigid substrate.The analysis of the temporal evolution of the system revealed an intermittent cracking activity which corresponds to the acoustic outbreaks of real experiments.

Discrete element simulations of the breakup of a shrinking layer
To mimic the typical experimental conditions, in the model a disc-shaped specimen of radius R is considered, which is discretized on a random lattice of space filling convex polygons representing mesoscopic pieces of the layer (see Fig. 1(c) for illustration).Polygons adjacent to each other in the initial tessellation are connected by breakable elastic beams of Young modulus E. Shrinking is represented by gradually reducing the natural length of beams with a fixed rate s, which results in a linearly increasing homogeneous strain ε = st as time t elapses.To capture the adhesion of the layer to a substrate, the middle point of the polygons is attached to the underlying plane by an elastic spring of stiffness D s , which is stress free in the initial configuration.As a consequence of the restrained shrinking of the layer, beams get overstretched and break according to a physical breaking rule (see Methods for details).Adhesion springs do not break so that no delamination of the layer can occur in the model.
To generate the time evolution of the particle system we carry out discrete element simulations solving the equation of motion of the polygons 26,27 (see Methods).The breaking criterion is evaluated in each iteration step and those beams which fulfil the condition are removed from the system.Cracks are formed by adjacent broken beams along the edges of polygons.The model naturally captures the simultaneous propagation and interaction of cracks, which eventually leads to the breakup of the layer into fragments.A snapshot of the evolving crack pattern is presented in Fig. 1c, where for clarity the beams are not shown.Further details of the model construction and the value of material parameters can be found in Ref. 28 .The model has been applied to study the shrinkage induced cracking of a layer, where it successfully reproduced the experimental findings on the structure of the planar crack pattern, and on the shape and statistics of the mass of fragments 28,29 .

Crackling noise during shrinkage induced cracking
As the layer shrinks cohesive contacts break and result in the nucleation of micro-cracks along the common edge of adjacent polygons.The breaking thresholds have fixed values for all the beams, however, the random tessellation introduces structural disorder into the material, which gives rise to a random nucleation of microcracks in the homogeneous deformation field of the layer.To characterize the temporal evolution of the cracking process, for each broken beam we record its position r b i and breaking time t b i from which the entire damage evolution of the layer can be reconstructed.In experiments crackling noise sensors like acoustic emission (AE) devices accumulate the effect of elastic waves generated by simultaneous cracking over extended spatial regions of the specimen 30 .Hence, we characterize the overall evolution of the breaking activity of the shrinking layer by the number of micro-cracks n b that occurred in equally sized time bins, which is plotted in Fig. 1a as a function of time t.The sudden rise and peak of the n b (t) curve mark the point of the dynamics where the deformation becomes sufficiently high to drive the growth of nucleated cracks in random directions.Growing cracks gradually merge and form a connected crack network along which the layers falls apart into a large number of fragments.Further shrinking results in cracks inside fragments splitting them into two pieces.For snapshots of the time evolution of the cracking layer see Fig. S1 of Supplementary Note 1. Figure 1c demonstrates that the emerging crack network has an isotropic cellular structure resulting in a nearly polygonal shape of the fragments 28 in agreement with measurements on various types of shrinking induced cracking phenomena 31,32 .The stochastic nature of the n b (t) curve in Fig. 1a indicates that in spite of the slow homogeneous increase of the deformation of the layer, damage accumulation through the nucleation and growth of cracks proceeds in an intermittent way.In order to identify bursts of the breaking activity we introduce a correlation time t c and assume that two consecutive micro-cracks (beam breakings) of time t b i and t b i+1 belong to the same trail of breakings if they follow each other within t c , i.e. if the condition t b i+1 − t b i < t c holds.This way the temporal sequence of micro-cracks is decomposed into bursts, i.e. pulses of cracking activity with a large number of breakings separated by silent periods where no breaking occurs.This is demonstrated in Fig. 1b where pulses of a small segment of the time series of Fig. 1a are highlighted.The size of a pulse is the number of beams failed in the trail of breaking events, the pulse duration T follows as the difference between the time of the last t b last and first t b first beam breaking T = t b last − t b first , while the waiting time t W between consecutive pulses is the duration of the silent period between them.These pulses of the global cracking activity are analogous to the acoustic outbreaks of real experiments on restrained shrinking so that the statistics of pulse quantities can be compared to the corresponding results of AE measurements 19,20,22,23 .The probability distribution p(�) of pulse sizes is presented in Fig. 2a for four different system sizes R at a fixed shrinking rate s along with two additional data sets at a lower and a higher value of s for R = 90.
A highly complex behaviour of p(�) is observed: at the smallest system size R = 30 the size distribution p(�) exhibits a power law decrease followed by an exponential cutoff.However, for larger systems pulses grow to larger sizes and a crossover occurs to a second power law regime of a lower exponent for large pulses, while the slope of the curve remains practically constant in the regime of small .Increasing the shrinking rate s at a fixed system size R has a similar effect, i.e. for higher s larger pulses are obtained accompanied by a crossover to a second power law of a lower exponent.To give a quantitative description of the evolution of the pulse size distribution p(�) we use the functional form to fit the curves in Fig. 2a.In Eq. (1) we assume that the size distribution p(�) is the sum of two power laws with exponential cutoffs, where τ l and τ h are the exponents in the low and high pulse size regimes, the corresponding cutoff sizes are denoted by l and h , while the parameters A and B control the relative contributions of the two terms in the distribution.The continuous lines in Fig. 2a represent fits with Eq. (1) obtained in such a way that τ l was fixed to the value τ l = 1.95 ± 0.05 of R = 30 , where the distribution is uni-modal without crossover, using only the other exponent τ h as a fitting parameter.It can be observed that the value of τ h decreases from the vicinity of τ h = 1.95 ± 0.06 towards 0 as R and s increase.The exponent δ controls the shape of the cutoff of p(�) having the same value δ = 1.50 ± 0.06 in all cases.Figure 2b demonstrates that the probability distribu- tion of the duration p(T) undergoes a similar evolution, i.e. in the regime of large durations, p(T) can always be approximated by a single power law followed by an exponential cutoff however, increasing R and s the cutoff duration T 0 shifts to higher values, while the power law exponent α gradu- ally decreases from α = 1.4 ± 0.06 to α = 0.8 ± 0.04 in the parameter regime considered.The value of the cutoff www.nature.com/scientificreports/exponent δ varies between 1.5 and 2.0.The waiting time distribution p(t W ) characterizes the statistics of silent periods between consecutive pulses.Figure 3 in Supplementary Note 3 presents that p(t W ) has a functional form similar to the duration distributions, i.e. a power law functional form is obtained with an exponential cutoff.
Increasing the shrinking rate s and the system size R the time gap between consecutive pulses decreases, hence, the exponent z of the power law regime of p(t W ) increasis while the characteristic waiting time controlling the cutoff decreases (see Supplementary Note 3).

Temporal superposition of local avalanches
In order to understand the origin of the complex evolution of pulse statistics with the system size R and driving rate s, it is important to emphasize that arranging micro-cracking events into pulses of activity is strictly based on the temporal sequence of micro-cracks.This is to mimic that in real experiments acoustic sensors integrate the signal of nucleating and propagating cracks over extended spatial regions which may involve even the entire specimen 30,33,34 .Hence, it can be expected that pulses observed in the shrinking layer emerge as the superposition of a large number of uncorrelated local avalanches of micro-cracks which are randomly scattered all over the sample.To prove that this mechanism is indeed responsible for the observed non-universality, we identify local avalanches as temporally and spatially correlated trails of micro-cracks, which have to be strictly distinguished from global activity pulses.Local avalanches always start from a single breaking bond and evolve as a spatially connected set of micro-cracks which follow each other within the correlation time t c and length l c according to the conditions t b i+1 − t b i < t c and r b i+1 − r b i < l c .In the calculations l c was set to be the double of the average poly- gon size l 0 of the discretization, while t c is the characteristic time scale of load redistribution determined by the material properties of the layer.Computer simulations revealed that local avalanches identified by the algorithm are intermittent steps of growing cracks spanning between junction points of the crack network (see Fig. 1c).
It can be observed in Fig. 3(a,b) that both the size and duration T of local avalanches are power law distributed with an exponential cutoff.It is important to emphasize that local avalanches exhibit a high degree of robustness, i.e. for both distributions p(�) and p(T) the results obtained at different system sizes R and driving rates s fall on a master curve.Although, the ranges spanned by and T are limited, i.e. the cutoff avalanche size 0 and duration T 0 are smaller than their counterparts of pulses, the curves can be very well fitted by the func- tional form of Eq. ( 2) with the power law exponents τ l = 1.95 ± 0.06 and α l = 1.4 ± 0.04 .Local avalanches are determined by the material properties of the shrinking layer such as the ratio of the stiffness of adhesion springs and beam elements D s /E , which controls the characteristic length scale of stress heterogeneities in the system.www.nature.com/scientificreports/ Figure S2 of Supplementary Note 2 demonstrates that at lower values of D s /E local avalanches grow to larger sizes and durations, however, the exponents τ l and α l remain the same.Due to the quenched structural disorder of the material, the homogeneous deformation field gives rise to a random pop up of local avalanches without any temporal correlation as the layer shrinks.It follows that the temporal occurrence of local avalanches can be approximated by a Poissonian process, where the effect of the system size R and of the driving rate s is that in a larger system a higher number of weak spots are available where cracks can nucleate, and at higher strain rates s the deformation increments s • t needed to trigger the next local avalanche is faster reached, respectively.As a consequence, in a homogeneous system the average nucleation rate r of avalanches must be proportional to the volume, here the area, of the system R 2 , and to the shrinking rate s so that �r� ∼ sR 2 holds.Figure 3c demonstrates that rescaling the waiting time distributions p(t W ) of local avalanches with sR 2 , distributions obtained at different system sizes R and shrinking rates s collapse on a master curve.The scaling function can be well approximated by an exponential as it is expected for Poissonian processes.
The results imply that the non-universal statistics of global activity pulses is caused by the superposition of temporally overlapping local avalanches of cracking, which randomly pop up in the system without any correlation.The relevance of overlapping can be quantified by the ratio of the average waiting time t W and of the average duration T of local avalanches.Figure 3d demonstrates that plotting the ratio t W / T as a function of the product sR 2 , results of all the simulations fall on a decreasing power law of exponent 1 implying �t W � ∼ 1/sR 2 ∼ 1/ �r� , which confirms the homogeneous Poissonian statistics of the temporal occurrence of local avalanches.It can also be inferred from the figure that at sufficiently small system sizes R and strain rates s the average waiting time can be significantly larger than the avalanche duration �t W � / �T� ≫ 1 , which implies that local avalanches form a sequence of well separated events without hardly any possibility of overlap in time.However, the value �t W � / �T� ≪ 1 obtained for large R and s indicates a strong temporal overlap of local events.For very large values of the product sR 2 this overlap can be so dominating that practically all local avalanches may get superimposed into a single or just a few pulses.This is supported by the evolution of the waiting time distribution p(t W ) of pulses with the system size R and driving rate s, see Fig. S3 of Supplementary Note 3. Based on our numerical results, we assume a linear dependence of the characteristic exponents τ h and α of the distribu- tions of pulse quantities on the nucleation rate r of local avalanches where a and b are constants that can be obtained by fitting.Note that Eqs.(3, 4) capture the fact that in the limit of zero driving �r� → 0 the exponents τ h and α must converge to the corresponding exponents τ l and α l of local avalanches.Figure 4 demonstrates that the linear form Eqs. (3, 4) provides a reasonable description of the joint effect of the system size R and shrinking rate s on the exponents of the integrated signal with a = 0.0015 ± 0.0003 and b = 0.0010 ± 0.0004 so that the ratio of the two parameters a and b is a/b ≈ 1.5.

Discussion
We showed that the restrained shrinking of a homogeneous layer of disordered materials gives rise to an intermittent evolution of cracking which can be decomposed into pulses of activity separated by silent periods.The size and duration of pulses have a scale-free statistics, similarly to other crackling phenomena, however, with non-universal exponents depending on the system size and on the driving rate.Our calculations revealed an astonishing universal micro-scale process underlying the non-universal behaviour of the integrated signal.Most of cracking pulses as function of sR 2 .For both exponents an approximate linear dependence is obtained on the nucleation rate of local avalanches described by Eqs.(3, 4).Note that the number of pulses decreases with increasing nucleation rate sR 2 .As a consequence the error bars of the exponents increase from 0.04 to 0.15 from the lowest to the highest values of sR 2 .
notably, we demonstrated that the local avalanches of micro-cracking events are characterized by universal scaling laws.The homogeneity of the deformation field and the uncorrelated spatial disorder of the layer give rise to a Poissonian process of the temporal occurrence of local avalanches.The observed non-universality of the global cracking activity is caused by the superposition of temporally overlapping local avalanches of micro-cracks.The power law exponents of the distributions of pulse sizes and durations proved to decrease linearly from the corresponding exponents of local avalanches towards zero as the nucleation rate of avalanches increases.
The number of systematic acoustic emission studies on shrinking induced cracking is rather limited in the literature.AE measurements on shrinking concrete revealed non-universal power law exponents increasing as the sample dries 19,20 .Our calculations suggest that the decreasing shrinking rate during solidification decreases the nucleation rate of local cracking events, which reduces their temporal overlap, and hence, increases the exponents towards their counterparts of local avalanches.Intermittent growth of shrinking induced cracks has also been observed during the formation of columnar joints, e.g. for the Kilauean lava lakes by recording seismic waves 35 .However, no systematic data is available on the statistics of cracking events.In general, the chisel-like marks often observed on the surface of basalt columns, called striae, are traces of rapid advancement of the crack front triggered by slow shrinking 36,37 .Based on the simulation results we conjecture that the area of striae is composed of local avalanches of the crack front so that its distribution may depend on the cooling rate of lava which should be tested by field measurements.A somewhat similar situation has been observed in crack propagation measurements where global avalanches of the advancing front of a planar crack have been found to be composed of clusters of locally failing elements along the front, which can be identified e.g. by complementing AE measurements using high speed imaging techniques 38,39 .A non-universal behaviour of crackling activity was also observed in the case of Barkhausen noise where the superposition of local avalanches of the motion of magnetic domain walls was found to be responsible for the driving rate dependence of crackling noise exponents.Varying the rate of change of the driving magnetic field a linear behaviour of the size and duration exponents was found 40 .
During charging-discharging cycles of e.g.Li-ion batteries the lithiation-delithiation process causes volume change of the electrode material which eventually gives rise to cracking accompanied by acoustic outbreaks 22,23 .Our results can help to decipher the information encoded in the acoustic event series measured at different rates of charging-discharging, which may serve as a basis for the assessment of the state of batteries in operation.

Methods
Discrete element model: In the model a disc-shaped specimen of radius R is discretized in terms of randomly shaped convex polygons obtained by a Voronoi construction.We use a regularised Voronoi tessellation, i.e. first a square lattice is placed on the sample and basis points of the tessellation are thrown randomly and independently into the plaquettes of the lattice.This technique results in a well-defined average polygon size l 0 , which is used as the unit length of the model.Polygons adjacent to each other in the initial tessellation are connected by breakable elastic beams of Young modulus E. The length and cross-section of beams are determined by the geometry of the polygonal lattice, which introduces randomness into their physical properties, e.g.moduli and moments 28 .
Shrinking of the layer is captured by gradually reducing the natural length of beams with a fixed rate s, which results in a linearly increasing homogeneous strain ε = st as time t elapses.To take into account the adhesion of the layer to a substrate, the middle point of the polygons is attached to the underlying plane by an elastic spring of stiffness D s , which is stress free in the initial configuration.During the restrained shrinking of the layer the beam elements get overstretched and break according to the breaking rule where ε ij denotes the longitudinal strain of the beam connecting polygons i and j, while i and j are the bending angles of the two beam ends 28 .Eq. ( 5) takes into account that stretching and bending contribute to the breaking of beams and the relative importance of the breaking modes is controlled by the breaking thresholds ε th and th .The breaking parameters ε th and th have fixed values for all the beams so that the structural randomness introduced by the random discretization is the only source of disorder in the system.Adhesion springs are not allowed to break so that no delamination can occur in the model.Those polygons which are not connected by beams may overlap during their motion, and we apply a repulsive force between them proportional to the overlap area.
To generate the time evolution of the particle system we carry out discrete element simulations solving the equation of motion of the polygons for the translational and rotational degrees of freedom 26,27 .For the numerical solution we use a 5th order Predictor-Corrector scheme which provides sufficient precision and stability.To mimic the effect of sticking to the container wall, polygons along the circular boundary of the layer are fixed during the simulations.Further details of the model construction and the value of its parameters can be found in Ref. 28 .The model has been applied to study the shrinkage induced cracking of a layer, where it successfully reproduced the experimental findings on the structure of the planar crack pattern, and on the shape and statistics of the mass of fragments 28,29 .In the present simulations the sample radius R was varied in the range 30 ≤ R ≤ 120 so that the number of polygons of discs falls between ∼ 3.000 and ∼ 45.000.

Data availability
The datasets generated and analysed during the current study are available from the corresponding author on reasonable request. (5)

Figure 1 .
Figure 1.Temporal evolution (a, b) and spatial distribution (c) of shrinking induced cracking for a system of radius R = 120 (in unit of the average polygon size l 0 ).In (a) the number of breaking beams n b is plotted as a function of time t.(b) demonstrates for a small segment of the time series that the breaking activity of the system can be decomposed into distinct pulses the size of which is plotted at the time of their occurrence.A snapshot of the cracking layer is presented in (c) after the connected crack network emerged.Micro-cracks are colored according to the local avalanche of breakings they belong to.A magnified view in the bottom-right corner reveals that local avalanches are intermittent steps of the growth of cracks.

Figure 2 .
Figure 2. Probability distributions of the size p(�) (a) and duration p(T) (b) of cracking pulses for several different system sizes R and strain rates s.Continuous lines represent fits with Eq. (1) (a) and with Eq. (2) (b).In (b) t denotes the time step of DEM simulations.

Figure 3 .
Figure 3. Probability distribution of the size p(�) (a), duration p(T) (b), and waiting time p(t W ) (c) of local avalanches for several system sizes R and strain rates s.(c) demonstrates that the waiting time distributions p(t W ) rescaled with the product sR 2 fall on the top of each other.The bold line represents an approximate exponential function.(d) Ratio of the average waiting time t W and the average duration T of local avalanches as a function of sR 2 .Results obtained at different R and s values fall on the same decreasing power law of exponent 1 indicated by the continuous line.

( 3 )Figure 4 .
Figure 4.Power law exponents τ h and α of the probability distribution of the size p(�) and duration p(T) of cracking pulses as function of sR 2 .For both exponents an approximate linear dependence is obtained on the nucleation rate of local avalanches described by Eqs.(3, 4).Note that the number of pulses decreases with increasing nucleation rate sR 2 .As a consequence the error bars of the exponents increase from 0.04 to 0.15 from the lowest to the highest values of sR 2 .