Phase transitions and correlations in fracture processes where disorder and stress compete

We study the effect of the competition between disorder and stress enhancement in fracture processes using the local load sharing fiber bundle model, a model that hovers on the border between analytical tractability and numerical accessibility. We implement a disorder distribution with one adjustable parameter. The model undergoes a localization transition as a function of this parameter. We identify an order parameter for this transition and find that the system is in the localized phase over a finite range of values of the parameter bounded by a transition to the non-localized phase on both sides. The transition is first order at the lower transition and second order at the upper transition. The critical exponents characterizing the second order transition are close to those characterizing the percolation transition. We determine the spatiotemporal correlation function in the localized phase. It is characterized by two power laws as in invasion percolation. We find exponents that are consistent with the values found in that problem.

We study the effect of the competition between disorder and stress enhancement in fracture processes using the local load sharing fiber bundle model, a model that hovers on the border between analytical tractability and numerical accessibility. We implement a disorder distribution with one adjustable parameter. The model undergoes a localization transition as a function of this parameter. We identify an order parameter for this transition and find that the system is in the localized phase over a finite range of values of the parameter bounded by a transition to the non-localized phase on both sides. The transition is first order at the lower transition and second order at the upper transition. The critical exponents characterizing the second order transition are close to those characterizing the percolation transition. We determine the spatiotemporal correlation function in the localized phase. It is characterized by two power laws as in invasion percolation. We find exponents that are consistent with the values found in that problem.

I. INTRODUCTION
It has been known for a long time that heterogeneities make materials more resilient against failure under load by offsetting the point at which a microfracture becomes unstable [1]. Examples of materials where heterogeneities play an important role are concrete [2,3] and carbonfiber composites [4] which possess larger strength and toughness than the individual components [5,6]. Heterogeneities introduce spatial disorder in the local material strength and create spatially dependent dynamic stress field. Under sufficiently large external stress, a micro-crack appears at the weakest point of the material and creates higher stress at the crack tips. Whether this microcrack will grow to a catastrophic failure or different micro-cracks will appear with increasing stress, depends on the dynamic competition between local stress enhancement and local material strength. Therefore, understanding the role of disorder in detail is crucial for technological purposes, in order to optimize the material strength.
The interplay between disorder and dynamical effects during the breakdown process of brittle materials was the focus of much research within the statistical physics community during the eighties and nineties [7]. This research was based on a number of lattice-based models where the links would have a maximum load before they would fail drawn from some spatially uncorrelated distribution. For example the fuse model [8][9][10] consisted of a network of electrical resistors that would act as fuses, the central force model consisted of a network of freely rotating Hookean springs that would fail if the load they carried crossed a threshold [11], and the beam model [12] where the links in the network would be elastic beams that would fail if certain criteria was fulfilled [13]. We summarize the picture that emerged qualitatively in the following. There are two reasons for a failure to appear locally in a disordered material. Either it is due to the material being weak at that local point or it is due to the stress being high there. Imagine now loading the disordered material. At the beginning of the breakdown process, the material will fail where it is weak as there are no -or few -points with high stress. When local failures develop, spots with intense stress appear at the tips of microcracks. As the failure process proceeds, these high-stress spots will start dominating. As it was argued by Roux and Hansen, [14] disorder makes local failures repulsive whereas the stress field makes them attractive. Imagine that there is a first local failure. Draw a sphere around this first local failure and identify the weakest spot within the sphere. The larger the sphere is, the weaker the weakest spot within it will be. As a result, if the next local failure is due to the local weakness of the material, it will occur as far as possible from the first one. Hence, the failures are repulsive when they are caused by the disorder in the strength of the material. The effect of the stress field is opposite. The closer one is to a local failure, the higher the largest stress will be. This makes it more likely that the next failure will be close to the one that just appeared. Hence, local failures attract each other when they are due to high stress field. The result of this is that there is a competition between disorder and stress concentration throughout the failure process. Early in the failure process the disorder tends to dominate, resulting in local failures appearing distributed throughout the material. However, as the stress field takes over, there is localization ending with a single growing crack beginning to dominate. Depending on the disorder, localization may occur early or later in the process: The wider the disorder, the later on the onset of localization would occur. In the limit of infinite disorder [15][16][17], localization never sets in and the failure process is a screened percolation -i.e. a process where links fail at random given that they are connected in such a way that they carry stress. When the disorder is weak enough to cause a competition with the stresses, a phase diagram may be constructed showing the onset of localization as a function of the disorder [18][19][20].
It is the aim of this work to study the effect of the competition between disorder and stress enhancement using the fiber bundle model [21,22]. The advantage of using this model over other models is that it is not computationally demanding, leading to good statistics for large samples. Furthermore, it has a high level of analytical tractability. We describe the fiber bundle model in detail in Sec. II. It should be noted that Stormo et al. [23] studied the interplay between disorder and stress enhancement in the soft clamp model which is a more complex version of the fiber bundle model than the version we study here. Their conclusions differ from those we present here; this due to a very different way of analyzing the fracture process. We return to this work in the concluding section. We implement a disorder distribution that has one adjustable parameter. By tuning this parameter, we produce disorders with power law tails towards either zero strength or infinite strength. We focus in Sec. III on the localization transition that occurs for some value of the disorder parameter, see Fig. 2. For the range of the values of the disorder parameter that produces power law tails towards zero strength, we find a second order phase transition, whereas for the range of values for the disorder parameter that produces a power law tail towards infinite strength, we find a first order phase transition. We determine the values of critical exponents associated with the second order phase transition. They are close to those found in percolation. In Sec. IV we study the spatiotemporal correlation function first introduced by Furuberg et al. [24] in connection with invasion percolation. We determine the scaling exponents characterizing the correlation function and find them to be close to the values observed for invasion percolation. We relate the exponents to other exponents characterizing the transition using theory developed by Roux and Guyon [25], Maslov [26] and Gouyet [27]. We end by drawing our conclusions in Sec. V.

II. DESCRIPTION OF THE MODEL
The fiber bundle model is a model of fracture where one can control the range of stress enhancement upon the appearance of a crack. In this model, N fibers (Hookean springs) are placed between two clamps under external force F . Each fiber carries a force where κ and x are respectively the elastic constant and extension of the fiber. The extension of each fiber has a threshold (x t ), beyond which it fails and the load it was carrying is distributed among surviving fibers according to some pre-set rule. The stress distribution scheme models the local stress enhancement whereas the disorder in x t models the local material heterogeneity. If the load of the failed fibers is distributed over all the surviving fibers, there is no local stress enhancement. This is the equal load sharing (ELS) scheme. The ELS fiber bundle model was introduced by Peirce in 1926 [28] as a simple model for failure in fibrous materials. Daniels approached the ELS fiber bundle model as a problem in statistics in a seminal paper in 1945 [29]. If the load is distributed only to the nearby fibers, we are dealing with the local load sharing (LLS) model [30,31]. Here local stress enhancement competes with the local heterogeneity leading to localization. Sornette introduced the ELS fiber bundle model to the statistical physics community in 1992 [32]. Soon, the focus of this community was on the rich avalanche statistics that this model offer [33][34][35][36], which in the ELS case is analytically tractable. The sequential -or time -correlation between failure events in the fiber bundle model makes it possible to explore the brittle to ductile transition [37] within it, a well-studied phenomenon in material science. The spatial correlations between the failures, on the other hand, have not been studied in detail.
In the LLS fiber bundle model, the load carried by the failed fibers is distributed equally among their nearest intact neighbors. We define a crack as a cluster of s failed nearest-neighbor fibers. The perimeter of the crack is the set of p intact fibers that are nearest neighbors to the failed fibers constituting the crack. These nearest neighbors define the hull of the cluster [38]. The force on an intact fiber i at any instance is then give by where σ = F/N , the force per fiber. The summation runs over all cracks J(i) that are neighbors to fiber i. This redistribution scheme is independent of the lattice topology and it is also independent of the failure history In the inset, we plot the scaling relation given in Eq. (8) where we find |D| + c (∞) = 0.47 ± 0.04 and β = 0.14 ± 0.05. [39], that is, the complete stress field can be calculated from the present arrangement of intact and broken fibers without having to take into account the order in which the fibers failed.
Suppose t fibers have failed. We determine which fiber will fail at time t + 1 in the following way [22]. Let f 1 i be the force on fiber i if we set the average force on the fibers σ = 1. We then calculate which denotes the fiber j(t + 1) that fails at time t + 1. The force σ = σ(t + 1) at which this force fails is given by .
The failure thresholds x t of the fibers are assigned by generating a random number over unit interval and raising it to power D which corresponds to the cumulative distribution [13,20] This threshold distribution allows us to control the disorder by the value of D, higher value of |D| implies higher disorder. Furthermore, D > 0 and D < 0 respectively correspond to the distributions with power law tails towards weaker and stronger fibers, which, as we will see, make the failure dynamics very different.

III. LOCALIZATION TRANSITION
We implement the model on a square lattice of size N = L × L fibers. The growth of cracks for L = 256 at different disorders are shown in Fig. 1 where we see two distinct regimes, a single crack growth at low disorder (top row) and random failures at high disorder (bottom row). The color scale in this figure represents where in the breaking sequence a given fiber fails, parameterized by τ = t/N , which we will discuss in detail later. At very low disorder, the stress enhancement at the crack perimeter wins over the fiber strengths and the crack always grows from the perimeter of the existing crack. The local disorder at the perimeter makes the crack grow as an invasion percolation cluster [40].
At high disorder, the strength of fibers win over the stress enhancement and we see random clusters of failed fibers appear. These two cases appear to be same for both D > 0 and D < 0, however, a completely different picture emerges between these two regimes when we observe crack growth at intermediate disorders. In Fig.  2, we plot the crack density n c = N c /N , where N c is the number of cracks (clusters) when t fibers have broken, for different disorders as a function of τ = t/N . We see that n c stays close to 0 (N c = 1) for the whole failure process for small |D| values whereas for large |D|, n c increases with τ and reaches a maximum (n max c ), beyond which, the cracks start to coalesce as more fibers are broken. The position of the peaks at τ -axis, marked by the red dots in Fig. 2, vary continuously with disorder for D > 0 whereas the peak stays at almost same position for D < 0 after a discontinuity between the low and high |D| values. Moreover, for two disorders with same values of n max c , the peak appears much earlier for D > 0 compared to D < 0. This shows that individual cracks appear randomly for D > 0 whereas cracks grow in size together with the appearance of new cracks for D < 0. Intuitively, when the threshold distribution has a power-law tail towards strong bonds (D < 0), probability to find weak bonds at existing crack perimeters are higher, which makes existing cracks to grow. Whereas when power-law tail is towards the weak bonds (D > 0), probability to find strong bonds at the perimeter is high and new cracks appear at different positions than the perimeter.
To characterize this localization transition, we define an order parameter, where the plus and minus signs correspond to D > 0 and D < 0 respectively. Here τ where we have that γ = 0.67 ± 0.01 and τ + (|D| ∞ ) = 0.291±0.01. For D < 0, τ − m (|D|) converges to a maximum value more abruptly from which we find τ − (|D| ∞ ) = 0.27 as shown by the dashed line in the inset.
The transition for D > 0 is explored in Fig. 3, where we plot O + (|D|, L) as a function of disorder for different system sizes. Here O + vary continuously from 1 to 0 similar to a second order phase transition. We find the following scaling for O + , where |D| + c (∞) = 0.47 ± 0.04 and β ≈ 0.14 ± 0.05. The exponent β is close to 5/36, the order parameter exponent for percolation in two dimension [41,42]. However, 1/8, the exact value from Onsager solution of the twodimensional Ising model [43], is also within the precision with which we know β. The transition for D < 0 is presented in Fig. 4 where we plot O − (|D|, L) as a function of disorder. Similar to a first order phase transition, O − shows sharp transition from 0 to 1 as the system size L is increased. The plots for different L can be scaled as [44], with η = 0.7. This is shown in the inset of Fig. 4(a). The transition disorder |D| − c (L) is observed to be a decreasing function of L. To find the value of |D| − c as L → ∞, we use the scaling |D| − c (∞) − |D| − c (L) ∼ L −θ with |D| − c (∞) = 1.03 ± 0.02 and θ = 0.80 ± 0.01. This is shown in Fig.  4(b).

IV. SPATIOTEMPORAL CORRELATIONS
We now turn to the investigation of the spatial and temporal correlations during the breakdown process. First we measure the fractal dimension of single cracks at low disorder, defined as s max ∼ L d f where s max is size of the largest incipient infinite cluster. This is shown in Fig. 5 (a). We find d f = 1.82 ± 0.02, and note that this close to the fractal dimension of invasion percolation cluster with trapping [40,45]. The colors in Fig. 1 represent the sequence of fiber breaking and together with the positions of broken fibers, they show the spatiotemporal map of the failure process. For random failures at high disorder (bottom row), the pixels (each representing a fiber) of different colors are randomly mixed, whereas for the single crack growth (top row), the different colors are clustered, indicating localization and hence spatial and temporal correlations. To quantify this correlation, we measure pair correlation function G(r, ∆t), defined as follows: if a fiber at position r 0 breaks at a time t, then G(r, ∆t) provides the probability that another fiber at r 1 will break at time t + ∆t, where r = |r 1 − r 0 |. Here time t is measured in terms of the number of broken fibers. This correlation function was first introduced by Furu-berg et al. [24] for invasion percolation. We assume in the following that the G(r, ∆t) has the scaling form [24] G(r, ∆t) = r −1 g r z t , where z is the dynamic exponent. In Fig. 5, we plot rG(r, ∆t) for the single crack growth regime as a function of r z /t for D > 0 (b) and D < 0 (c). For z = d f , the single-parameter function g shows power-law behavior in both the large-and small-argument limits, with a peak at u = 1 or at t = r d f as shown in the figure. This implies that the most probable growth of the crack after time t occurs at a distance r d f , or, the most probable growth at a distance r occurs after a time t 1/d f . The exponent a in the long time range is found to be 1.44 ± 0.04 as shown in Fig. 5. This matches with the exponent for invasion percolation dynamics [24]. Recently, this scaling was observed experimentally for two-phase flow in porous medium during slow drainage [46] with the exponent a ranging between 1.44 and 1.73. For the power law in the small-argument limit, we obtain the exponent b = 0.47 ± 0.06, which is different from 0.68 obtained by Furuberg et al. for invasion percolation dynamics [24]. The exponents and the errors mentioned above are obtained from the least square fitting of the numerical data. The plots with multiple data file, such as figure 3 inset, figure 5 (b) and (c), the error mentioned is the average of the errors corresponding to individual data.
We may relate the exponents a and b to the exponents controlling the avalanche structure of the failure process following Roux and Guyon [25], Maslov [26] and Gouyet [27]. The minimum force applied to the fiber bundle for fiber number t + 1 to fail is σ(t + 1) given in Eq. (4). We define a forward avalanche (or burst) of size Θ starting at time t to be [26,33] σ(t + k) < σ(t) for k < Θ and σ(t + Θ) ≥ σ(t) . (12) A backwards avalanche (or burst) we define as σ(t − k) < σ(t) for k < Θ and σ(t − Θ) ≥ σ(t) . (13) Roux and Guyon [25] define two distributions, P ∆t and Q Θ (r). The first one, P ∆t , gives the distribution of the smallest avalanches Θ that pass through σ(t) and σ(t − ∆t) averaged over σ(t). It obeys the power law, where H(x) is the Heaviside function which is one when x > 0 and zero otherwise. The second one, Q Θ , gives the distribution of distances between failing fibers within an avalanche of size Θ. Roux and Guyon assume it to follow the power law The spatiotemporal correlation function G(r, ∆t) may then be constructed from these two probability distributions, which integrates to Hence, we have in Eq. (11) that a = (1 + ζ)/d f and b = c b − 1. Roux and Guyon suggest ζ = d f − 1 making a = 1. We propose here that ζ = d h , the fractal dimension of the hull [47], defined as the set of the sites that are connected to the cluster and at the neighbor to the surroundings. In the LLS fiber bundle model, the stress of a crack is re-distributed to its perimeter, the neighboring unbroken fibers of the crack. In Fig. 5, we plot this largest perimeter size p max with the system size L and from the slope, we find d h = 1.74 ± 0.02. This value matches with the value of d h for percolation hull, known to be 7/4 by the relation with ν by d h = 1 + 1/ν, where ν is the correlation length exponent which is equal to 4/3 in two-dimensional percolation [48,49]. Hence, we find that With d f = 1.82 ± 0.02 and d h = 1.74 ± 0.02, we find a = 1.51 ± 0.04, which is close to the observed value.
Maslov [26] related the backwards avalanche exponent c b to the exponent c f that governs the probability to find a forward avalanche of size Θ when the stress is σ 0 via the expression Furthermore, Goyuet in turn related c f to d h , d f and ν, which leads b = 0.46. This is also in accordance with the value we observe, b = 0.47.

V. CONCLUSION
We have here studied the effect of the competition between disorder and stress enhancement using the local load sharing fiber bundle model. We have done this by varying the disorder using a threshold distribution controlled by a single-parameter D. When D is positive, the threshold distribution is a power law towards infinitely weak elements, and when D is negative, the distribution is a power law towards infinitely strong elements. By defining an order parameter distinguishing between the localized and non-localized phases, we find two phase transitions, one at D = D + c = 0.47 ± 0.04 and one at D = D − c = −1.03 ± 0.02. The transition for D > 0 is second order and the transition for D < 0 is first order. The second order transition is governed by critical exponents that are consistent with percolation, but also those found for the two-dimensional Ising model. We then studied the spatiotemporal correlation function, finding the same behavior as first seen by Furuberg et al. [24] in an invasion percolation context. Some numerical values for the exponents controlling the correlation function are different than found by Furuberg et al., which, by following the scaling analysis of Roux and Guyon [25], Maslov [26] and Gouyet [27], we successfully related to other exponents describing the geometry of the cracks.
The conclusions we present here differ from those that Stormo et al. [23] presented. In that work, based on the soft-clamp fiber bundle model, the threshold distribution was the uniform one on the unit interval, i.e. D = 1. The parameter that was varied, was the elastic constant of the clamps to which the fibers are connected. For a given value of this parameter, the breakdown process would proceed in the beginning as an uncorrelated percolation process up to a certain point at which localization would set in. From this point on, the breakdown process would continue through the growth of a single cluster of broken fibers. This transition would not be a phase transition, but a crossover.
It would be of great interest to repeat the analysis we have presented here using the soft-clamp fiber bundle model. Only then we will be able to distinguish what is model dependent in our conclusions and what is not.