# Scale Invariance and Lack of Self-Averaging in Fragmentation

## Abstract

We derive exact statistical properties of a class of recursive fragmentation processes. We show that introducing a fragmentation probability $`0<p<1`$ leads to a purely algebraic size distribution in one dimension, $`P(x)x^{2p}`$. In $`d`$ dimensions, the volume distribution diverges algebraically in the small fragment limit, $`P(V)V^\gamma `$ with $`\gamma =2p^{1/d}`$. Hence, the entire range of exponents allowed by mass conservation is realized. We demonstrate that this fragmentation process is non-self-averaging. Specifically, the moments $`Y_\alpha =_ix_i^\alpha `$ exhibit significant fluctuations even in the thermodynamic limit.

PACS numbers: 05.40.+j, 64.60.Ak, 62.20.Mk

Numerous physical phenomena are characterized by a set of variables, say $`\{x_j\}`$, which evolve according to a random process, and are subject to the conservation law $`_jx_j=\mathrm{const}`$. An important example is fragmentation, with applications ranging from geology and fracture to the breakup of liquid droplets and atomic nuclei. Other examples include spin glasses, where $`x_j`$ represents the equilibrium probability of finding the system in the $`j^{\mathrm{th}}`$ valley, genetic populations, where $`x_j`$ is the frequency of the $`j^{\mathrm{th}}`$ allele, and random Boolean networks.

In most cases, stochasticity governs both the way in which fragments are produced, and the number of fragmentation events they experience. For example, in fragmentation of solid objects due to impact with a hard surface fragments may bounce several times before coming to a rest . The typical number of fragmentation events may vary greatly as it depends on the initial kinetic energy. Another seemingly unrelated example is DNA segmentation algorithms where homogeneous subsequences are produced recursively from a inhomogeneous sequence until a predefined homogeneity level is reached . Here, the number of segmentation events is determined by the degree of homogeneity of the original sequence.

In this study, we examine fragmentation models with two types of objects – stable towards fragmentation and unstable. We show that the size distribution is algebraic, and that the entire range of power-laws allowed by the underlying conservation laws can be realized by tuning the fragmentation probability. Additionally, such processes are characterized by large sample to sample fluctuations, as seen from analysis of the moments of the fragment size distribution.

Specifically, we consider the following recursive fragmentation process. We start with the unit interval and choose a break point $`l`$ in $`[0,1]`$ with a uniform probability density. Then, with probability $`p`$, the interval is divided into two fragments of lengths $`l`$ and $`1l`$, while with probability $`q=1p`$, the interval becomes “frozen” and never fragmented again. If the interval is fragmented, we recursively apply the above fragmentation procedure to both of the resulting fragments.

First, let us examine the average total number of fragments, $`N`$. With probability $`q`$ a single fragment is produced, and with probability $`p`$ the process is repeated with two fragments. Hence $`N=q+2pN`$, yielding

$$N=\{\begin{array}{cc}q/(12p),\hfill & \text{if }p<1/2\text{;}\hfill \\ \mathrm{},\hfill & \text{if }p1/2\text{.}\hfill \end{array}$$
(1)   
The average total number of fragments becomes infinite at the critical point $`p_c=1/2`$, reflecting the critical nature of the underlying branching process.

Next, we study $`P(x)`$, the density of fragments of length $`x`$. The recursive nature of the process can be used to obtain the fragment length density

$$P(x)=q\delta (x1)+2p_x^1\frac{dy}{y}P\left(\frac{x}{y}\right).$$
(2)   
The second term indicates that a fragment can be created only from a larger fragment, and the $`y^1`$ kernel reflects the uniform fragmentation density. Eq. (2) can be solved by introducing the Mellin transform

$$M(s)=𝑑xx^{s1}P(x).$$
(3)   
Eqs. (2) and (3) yield $`M(s)=q+2ps^1M(s)`$ and as a result

$$M(s)=q\left[1+\frac{2p}{s2p}\right].$$
(4)   
The average total number $`M(1)=N`$ is consistent with Eq. (1), and the total fragment length $`M(2)=1`$ is conserved in accord with $`1=𝑑xxP(x)`$. (Here and in the following the integration is carried over the unit interval, i.e., $`0<x<1`$.) The inverse Mellin transform of Eq. (4) gives

$$P(x)=q\left[\delta (x1)+2px^{2p}\right].$$
(5)   
Apart from the obvious $`\delta `$-function, the length density is a purely algebraic function. In particular, the fragment distribution diverges algebraically in the limit of small fragments. Given such an algebraic divergence near the origin $`P(x)x^\gamma `$, length conservation restricts the exponent range to $`\gamma <2`$. In our case $`\gamma =2p`$, and since $`0<p<1`$, the entire range of acceptable exponents emerges by tuning the only control parameter $`p`$.

Interestingly, at the critical point $`p_c=\frac{1}{2}`$, the fragment length distribution becomes independent of the initial interval length. Starting from an interval of length $`L`$, Eq. (5) can be generalized to yield

$$P(x)=q\delta (xL)+2pqL^{12p}x^{2p}.$$
(6)   
Thus, the critical point may be detected by observing the point at which the segment distribution becomes independent of the original interval length.

The recursive fragmentation process can be generalized to $`d`$ dimensions. For instance, in two dimensions we start with the unit square, choose a point $`(x_1,x_2)`$ with a uniform probability density, and divide, with probability $`p`$, the original square into four rectangles of sizes $`x_1\times x_2`$, $`x_1\times (1x_2)`$, $`(1x_1)\times x_2`$, and $`(1x_1)\times (1x_2)`$. With probability $`q`$, the square becomes frozen and we never again attempt to fragment it. The process is repeated recursively whenever a new fragment is produced.

Let $`P(𝐱)`$, $`𝐱(x_1,\mathrm{},x_d)`$, be the probability density of fragments of size $`x_1\times \mathrm{}\times x_d`$. This quantity satisfies

$$P(𝐱)=q\delta (𝐱\mathrm{𝟏})+2^dp\frac{d𝐲}{y_1\mathrm{}y_d}P(\frac{x_1}{y_1},\mathrm{},\frac{x_d}{y_d}).$$
(7)   
with $`𝑑𝐲=𝑑y_1\mathrm{}𝑑y_d`$. Following the steps leading to Eq. (4), we find that the $`d`$-dimensional Mellin transform, defined by $`M(𝐬)=𝑑𝐱x_1^{s_11}\mathrm{}x_d^{s_d1}P(𝐱)`$ with the shorthand notation $`𝐬(s_1,\mathrm{},s_d)`$ obeys

$$M(𝐬)=q\left[1+\frac{\alpha ^d}{s_1\mathrm{}s_d\alpha ^d}\right],\mathrm{with}\alpha =2p^{1/d}.$$
(8)   
Eq. (8) gives the total average number of fragments, $`N=M(\mathrm{𝟏})=q/(12^dp)`$ if $`p<2^d`$ and $`N=\mathrm{}`$ if $`p2^d`$. One can also verify that the total volume $`M(\mathrm{𝟐})=1`$ is conserved. Interestingly, there is an additional infinite set of conserved quantities: all moments whose indices belong to the hyper-surface $`s_1^{}\mathrm{}s_d^{}=2^d`$ satisfy $`M(𝐬^{})=1`$. In a continuous time formulation of this process the same moments were found to be integrals of motion . The existence of an infinite number of conservation laws is surprising, because only the volume conservation has a clear physical justification.

Next, we study the volume density $`P(V)`$, defined by

$$P(V)=𝑑𝐱P(𝐱)\delta \left(Vx_1\mathrm{}x_d\right).$$
(9)   
The Mellin transform $`M(s)=𝑑VV^{s1}P(V)`$ can be obtained from Eq. (8) by setting $`s_i=s`$,

$$M(s)=q\left[1+\frac{\alpha ^d}{s^d\alpha ^d}\right].$$
(10)   
Using the $`d^{\mathrm{th}}`$ root of unity, $`\zeta =e^{2\pi i/d}`$, and the identity $`\frac{1}{s^d1}=\frac{1}{d}_{k=0}^{d1}\frac{\zeta ^k}{s\zeta ^k}`$, $`M(s)`$ can be expressed as a sum over simple poles at $`\alpha \zeta ^k`$. Consequently, the inverse Mellin transform is given by a linear combination of $`d`$ power laws

$$P(V)=q\left[\delta (V1)+\frac{\alpha }{d}\underset{k=0}{\overset{d1}{}}\zeta ^kV^{\alpha \zeta ^k}\right].$$
(11)   
One can verify that this expression equals its complex conjugate and hence, it is real. Additionally, the one-dimensional case (5) is recovered by setting $`d=1`$.

The small-volume tail of the distribution can be obtained by noting that the sum in Eq. (11) is dominated by the first term in the series, which leads to

$$P(V)AV^\gamma \text{as}V0,$$
(12)   
with $`\gamma =\alpha =2p^{1/d}`$ and $`A=\alpha q/d`$. Although the value of the exponent changes, the possible range of exponents for this process remains the same since $`0<2p^{1/d}<2`$ when $`0<p<1`$. In the infinite dimension limit, $`P(V)`$ becomes universal: $`P(V)V^2`$.

The leading behavior of $`P(V)`$ in the large size limit can be derived by using the Taylor expansion and the identity $`_{k=0}^{d1}\zeta ^{kn}=\delta _{n,0}`$ for $`n=0,\mathrm{},d1`$. One finds that in higher dimensions the volume distribution vanishes algebraically near its maximum value,

$$P(V)B_d(1V)^{d1}\text{as}V1,$$
(13)   
with $`B_d=\alpha ^d/(d1)!`$.

In fact, the entire multivariate fragment length density can be obtained explicitly. This can be achieved by expanding the geometric series

$`{\displaystyle \frac{\alpha ^d}{s_1\mathrm{}s_d\alpha ^d}}={\displaystyle \underset{n0}{}}{\displaystyle \underset{i=1}{\overset{d}{}}}\left({\displaystyle \frac{\alpha }{s_i}}\right)^{n+1},`$   
and performing the inverse Mellin transform for each variable separately. Using the transform $`𝑑xx^{s1}\left[\mathrm{ln}\frac{1}{x}\right]^n=n!s^{n1}`$ gives

$$P(𝐱)=q\left[\delta (𝐱\mathrm{𝟏})+\alpha ^dF_d(z)\right],$$
(14)   
with the shorthand notations

$`F_d(z)={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{z^n}{n!}}\right)^d,z=\alpha \left({\displaystyle \underset{i=1}{\overset{d}{}}}\mathrm{ln}{\displaystyle \frac{1}{x_i}}\right)^{1/d}.`$ (15)   
In two dimensions, $`F_2(z)=I_0(2z)`$ where $`I_0`$ is the modified Bessel function.

The small size behavior of $`P(𝐱)`$ can be obtained by using the steepest decent method. The leading tail behavior, $`F_d(z)(2\pi z)^{\frac{1d}{2}}e^{zd}`$ for $`z1`$, corresponds to the case when at least one of the lengths is small, i. e. $`x_i1`$. Returning to the original variables we see that the fragment distribution exhibits an unusual “log-stretched-exponential” behavior

$`P(𝐱)\left[{\displaystyle \underset{i=1}{\overset{d}{}}}\mathrm{ln}{\displaystyle \frac{1}{x_i}}\right]^{\frac{d1}{2d}}\mathrm{exp}\left[d\alpha \left({\displaystyle \underset{i=1}{\overset{d}{}}}\mathrm{ln}{\displaystyle \frac{1}{x_i}}\right)^{1/d}\right].`$ (16)   
The fragment distribution represents an average over infinitely many realizations of the fragmentation process; hence, it does not capture sample to sample fluctuations. These fluctuations are important in non-self-averaging systems, where they do not vanish in the thermodynamic limit. Useful quantities for characterizing such fluctuations are the moments $`Y_\alpha `$

$`Y_\alpha ={\displaystyle \underset{i}{}}x_i^\alpha ,`$ (17)   
where the sum runs over all fragments.

We are interested in the average values $`Y_\alpha `$ and $`Y_\alpha Y_\beta `$. For integer $`\alpha `$, $`Y_\alpha `$ is the probability that $`\alpha `$ points randomly chosen in the unit interval belong to the same fragment. The expected value of $`Y_\alpha `$ satisfies

$$Y_\alpha =q+pY_\alpha 𝑑y\left[y^\alpha +(1y)^\alpha \right].$$
(18)   
The first term corresponds to the case where the unit interval is not fragmented, and the second term describes the situation when at least one fragmentation event has occurred. Eq. (18) gives

$$Y_\alpha =q\left[1+\frac{2p}{\alpha +12p}\right]$$
(19)   
if $`\alpha >2p1`$, and $`Y_\alpha =\mathrm{}`$ if $`\alpha 2p1`$. As expected, Eq. (19) agrees with the moments of $`P(x)`$ obtained by integrating Eq. (5), $`Y_\alpha =𝑑xx^\alpha P(x)`$.

However, higher order averages such as $`Y_\alpha Y_\beta `$ do not follow directly from the fragment density. For integer $`\alpha `$ and $`\beta `$, $`Y_\alpha Y_\beta `$ is the probability that, if $`\alpha +\beta `$ points are chosen at random, the first $`\alpha `$ points all lie on the same fragment, and the last $`\beta `$ points all lie on another (possibly the same) fragment. This quantity satisfies

$`Y_\alpha Y_\beta `$ $`=`$ $`q+pY_\alpha Y_\beta {\displaystyle 𝑑y\left[y^{\alpha +\beta }+(1y)^{\alpha +\beta }\right]}`$ (20)   
$`+`$ $`pY_\alpha Y_\beta {\displaystyle 𝑑y\left[y^\alpha (1y)^\beta +(1y)^\alpha y^\beta \right]},`$ (21)   
yielding

$`Y_\alpha Y_\beta `$ $`=`$ $`q+{\displaystyle \frac{2pq}{\alpha +\beta +12p}}`$ (22)   
$`+`$ $`2p{\displaystyle \frac{\mathrm{\Gamma }(\alpha +1)\mathrm{\Gamma }(\beta +1)}{\mathrm{\Gamma }(\alpha +\beta +1)}}{\displaystyle \frac{Y_\alpha Y_\beta }{\alpha +\beta +12p}}`$ (23)   
when $`\alpha ,\beta ,\alpha +\beta >2p1`$, and $`Y_\alpha Y_\beta =\mathrm{}`$ otherwise.

Eq. (22) shows that $`Y_\alpha Y_\beta Y_\alpha Y_\beta `$, and in particular, $`Y_\alpha ^2Y_\alpha ^2`$. Therefore, fluctuations in $`Y_\alpha `$ do not vanish in the thermodynamic limit, and the recursive fragmentation process is non-self-averaging. While for $`p<1/2`$ non-self-averaging behavior is expected because the average number of fragments is finite, the emergence of non-self-averaging for $`p>1/2`$ is surprising. Hence, statistical properties obtained by averaging over all realizations are insufficient to probe sample to sample fluctuations.

In principle, higher order averages such as $`Y_\alpha ^n`$ can be calculated recursively by the procedure outlined above. The resulting expressions are cumbersome and not terribly illuminating. Instead, one may study the distribution $`Q_\alpha (Y)`$ which obeys

$`Q_\alpha (Y)=q\delta (Y1)`$ (24)   
$`+p{\displaystyle 𝑑l_0^Y𝑑Z\frac{1}{l^\alpha }Q_\alpha \left(\frac{Z}{l^\alpha }\right)\frac{1}{(1l)^\alpha }Q_\alpha \left(\frac{YZ}{(1l)^\alpha }\right)}.`$ (25)   
In addition to the recursive nature of the process, we have employed extensivity, i.e, $`Y_\alpha L^\alpha `$ in an interval of length $`L`$.

Clearly, $`Y_0=N`$ and $`Y_11`$, and therefore $`Q_1(Y)=\delta (Y1)`$ and $`Q_0(N)`$ can be determined analytically as well. Generally, different behaviors emerge for $`\alpha >1`$ and $`\alpha <1`$. We concentrate on the former case where the support of the distribution $`Q_\alpha (Y)`$ is the interval . The Laplace transform, $`R_\alpha (\lambda )=_0^1𝑑Ye^{\lambda Y}Q_\alpha (Y)`$, obeys

$$R_\alpha (\lambda )=qe^\lambda +p_0^1𝑑lR_\alpha \left[\lambda l^\alpha \right]R_\alpha \left[\lambda (1l)^\alpha \right].$$
(26)   
The behavior of $`Q_\alpha (Y)`$ in the limit $`Y0`$ is reflected by the asymptotics of $`R_\alpha (\lambda )`$ as $`\lambda \mathrm{}`$. Substituting $`R_\alpha (\lambda )\mathrm{exp}(A\lambda ^\beta )`$ into both sides of Eq. (26), evaluating the integral using steepest decent and equating the left and right hand sides gives $`\beta =1/\alpha `$. Consequently, we find that the distribution has an essential singularity near the origin

$$Q_\alpha (Y_\alpha )\mathrm{exp}\left[BY_\alpha ^{\frac{1}{\alpha 1}}\right],Y_\alpha 0.$$
(27)   
Extremal properties can be viewed as an additional probe of sample to sample fluctuations. We thus consider $`(x)`$, the length density of the largest fragment. For a self-averaging process with an infinite number of fragments, one expects $`(x)\delta (x)`$ in the thermodynamic limit. To see that $`(x)`$ is non-trivial for any $`p`$, let us first determine $`(x)`$ for $`x1/2`$. In this region,

$`(x)=q\delta (x1)+p{\displaystyle _x^1}{\displaystyle \frac{dy}{y}}\left({\displaystyle \frac{x}{y}}\right).`$ (28)   
If the original unit interval has not been fragmented, the largest fragment is obviously the unit interval. If the first fragmentation is performed, only one of the two resulting fragments can be larger than $`x>1/2`$. Therefore, only subsequent breaking of this fragment (of length $`y>x`$) can contribute to $`(x)`$, which explains Eq. (28).

Eq. (28) is similar to Eq. (2), and can be solved by the same technique to give

$`(x)=q\delta (x1)+pqx^p\mathrm{for}x1/2.`$ (29)   
In the complementary case of $`x<1/2`$, $`(x)`$ satisfies

$`(x)`$ $`=`$ $`p{\displaystyle \underset{1x}{\overset{1}{}}}{\displaystyle \frac{dy}{y}}\left({\displaystyle \frac{x}{y}}\right)+p{\displaystyle \underset{1/2}{\overset{1x}{}}}{\displaystyle \frac{dy}{y}}\left({\displaystyle \frac{x}{y}}\right)_{}\left({\displaystyle \frac{x}{1y}}\right)`$ (30)   
$`+`$ $`p{\displaystyle \underset{1/2}{\overset{1x}{}}}{\displaystyle \frac{dy}{y}}\left({\displaystyle \frac{x}{1y}}\right)_{}\left({\displaystyle \frac{x}{y}}\right).`$ (31)   
The first term on the right-hand side of this equation is constructed as in Eq. (28): if we first break the unit interval into two fragments of lengths $`y>1/2`$ and $`1y`$, then for $`1y<x`$ the longest fragment is produced by breaking the fragment of length $`y`$. The next two terms describe the situation when $`1y>x`$, so the longest fragment can arise out of breaking any of the two fragments. The factors $`_{}(u)=_0^u𝑑v(v)`$ guarantee that the longest fragment of length $`x`$ comes from the fragment of length $`v`$ in the first generation. Since $`(x)`$ obeys different equations in different regions, it looses analyticity on the boundaries. Namely, $`(x)`$ possesses an infinite set of singularities at $`x=1/k`$ which become weaker as $`k`$ increases. Similar singularities underly extremal properties of a number of random processes including random walks, spin glasses, random maps, and random trees .

In summary, we have found that recursive fragmentation is scale free, i.e., the fragment length distribution is purely algebraic. In higher dimensions, the volume distribution is a linear combination of $`d`$ power laws, and consequently, an algebraic divergence characterizes the small-fragment tail of the distribution. A number of recent impact fragmentation experiments reported algebraic mass distributions with the corresponding exponents ranging from $`1`$ to $`2`$ . It will be interesting to further examine whether our simplified model is appropriate for describing fragmentation of solid objects.

We have also found that the recursive fragmentation process exhibits a number of features that arise in other complex and disordered systems, such as non-self-averaging behavior and the existence of an infinite number of singularities in the distribution of the largest fragment. These features indicate that even in the thermodynamic limit sample to sample fluctuations remain, and that knowledge of first order averages may not be sufficient for characterizing the system. Our 1D model is equivalent to applying the aforementioned DNA segmentation algorithm to a random sequence. It will be interesting to study self-averaging and extremal properties of DNA sequences, which are known to have commonalities with disordered systems. Indeed, if these subtle features are found for genetic sequences as well, this would suggest that much caution should be exercised in statistical analysis of DNA.

We are thankful to S. Redner and O. Weiss for useful discussions, and to DOE, NSF, ARO, NIH, and DFG for financial support.