A model for the fragmentation kinetics of crumpled thin sheets

As a confined thin sheet crumples, it spontaneously segments into flat facets delimited by a network of ridges. Despite the apparent disorder of this process, statistical properties of crumpled sheets exhibit striking reproducibility. Experiments have shown that the total crease length accrues logarithmically when repeatedly compacting and unfolding a sheet of paper. Here, we offer insight to this unexpected result by exploring the correspondence between crumpling and fragmentation processes. We identify a physical model for the evolution of facet area and ridge length distributions of crumpled sheets, and propose a mechanism for re-fragmentation driven by geometric frustration. This mechanism establishes a feedback loop in which the facet size distribution informs the subsequent rate of fragmentation under repeated confinement, thereby producing a new size distribution. We then demonstrate the capacity of this model to reproduce the characteristic logarithmic scaling of total crease length, thereby supplying a missing physical basis for the observed phenomenon.


SUPPLEMENTARY FIGURE 2
Supplementary Figure 2. Facet area distributions from manual segmentation. Distributions of facet areas x normalized by the mean area s respectively for each sample; ξ = x/s. The solid line shows the best fit to Supplementary Equation (3b), where the parameter a is calculated from the relation a(t) = t/τ presented in the main text with single fitting parameter τ across the entire dataset. Marker colors correspond to different values of∆, as indicated by the colorbar.

SUPPLEMENTARY FIGURE 6
Supplementary Figure 6. Estimation of overall breakup rate r(x). For n = 1 (top row), the fraction of facets present at crumpling iteration n = 1 which fragment by n = 2, as a function of initial area x, for samples with∆ = 0.45, 0.36, 0.27, 0.18, 0.09, and 0.045 (across). For n = 2 (bottom row), the fraction of all facets present at crumpling iteration n = 2 which fragment by n = 3 for the same samples. The sample with∆ = 0.63 had too few facets to form a sufficient representative sample. Error bars denote the standard deviation of the fragmentation probability if the fragmentation of each facet is regarded as a Bernoulli trial, with the fraction of fragmented facets taken as the success probability within each histogram bin. The dashed line corresponds to √ x. Marker colors correspond to different values of∆, as indicated by the colorbar. As noted in the text, samples at small values of∆, or high compaction, likely undergo a succession of fragmentation events between n = i and n = i + 1, and are thus poorer indicators of the statistics of single breakup events. Samples at large values of∆ are more likely resolve single breakup events, but have a lower population of facets from which to build the distribution. The choice of overall breakup rate r(x) = x 1/2 was motivated both by the stronger power law behavior at high∆, as well as its tractability in our analytical model.

SUPPLEMENTARY FIGURE 7
Supplementary Figure 7. Estimation of conditional breakup probability f (x|y). For n = 1 (top row), the probability density function ρ(x/y) of facet areas x present in crumpling iteration n = 2 normalized by their parent facet's area y from n = 1, for samples with∆ = 0.45, 0.36, 0.27, 0.18, 0.09, and 0.045 (across). For n = 2 (bottom row), ρ(x/y) of facet areas in n = 3 normalized by their parent facet's area from n = 2. Marker colors correspond to different values of∆, as indicated by the colorbar. The sample with∆ = 0.63 had too few facets to form a sufficient representative sample and is excluded here. As noted in Supplementary Fig. 6, samples at small values of∆ likely undergo a succession of fragmentation events between crumples, and thus their distributions resemble the more mature facet distributions observed at later n, as in Supplementary Fig. 2.
Supplementary Figure 8. Convergence of asymptotic approximation for survival function S Z (w;θ). Plot of Supplementary Equation (5) (dashed line) against the asymptotic approximation of Supplementary Equation (7) valid at large number of steps k (solid line), at 2k = 4, 8, 16, 32, and 64 steps from left to right, respectively. The shape parameter a = 1 for all cases, and t is appropriately determined from the relation t = 2k(a + 1)/L0. Curves of the asymptotic approximation are colored by the value of t, as indicated by the colorbar. The approximation shows increasing agreement with the exact solution for larger k, as anticipated.

SUPPLEMENTARY FIGURE 9
Supplementary Figure 9. Numerical simulation of random walks and exact escape probability compared to analytical approximation. a Plot of Supplementary Equation (5), which presents an analytical estimate of the fraction of escaped random walks looking only at final displacements (dashed line), compared against the more accurate result obtained by assessing all intermediate displacements through numerical simulation (solid line). For the latter, simulations of 50, 000 random walks, with gamma-distributed steps sampled according to Supplementary Equation (4) with shape parameter a = 1, are performed with 2k = 4, 8, 16, 32, and 64 steps, and t appropriately determined from the relation t = 2k(a + 1)/L0. Curves corresponding to numerical simulation are colored by the value of t, as indicated by the colorbar. The analytical approximation systematically underestimates the number of escaped walks by approximately a constant factor; thus, SZ (w; a, θ) remains proportional to the change in t over a crumpling iteration in both the analytically approximate and numerically calculated forms, allowing variation in the constant of proportionality. b The error between each pair of curves presented in a.

SUPPLEMENTARY METHODS
Manual segmentation. As noted in the main text, the final segmentation of all collected crease networks was performed by hand. These segmentations and the distributions of scaled facet areas are provided in Supplementary Figs. 1 & 2, respectively. We recognize in the Methods section of the main text that creases can soften over repeated crumples, and the result of unfolding and scanning between crumples could possibly contribute to the appearance of healing. To ensure that the crease patterns studied suitably fit the framework of a fragmentation model, we perform a simple analysis to affirm that healing is indeed very minimal. We make a quantitative prediction about the extent of healing by overlaying two manually segmented crease patterns from successive crumples n − 1 and n of the same sheet, and measuring the length of creases present in n − 1 which do not appear in n. We observe that the percent of healed creases make up less than 5% of any crease pattern; moreover, the fraction is typically within 1% for moderately to highly dense patterns. Thus, we conclude that healing is a small effect which does not greatly impact the dominant trends in the data.
Watershed segmentation. Prior to manual segmentation, an automated method using the watershed algorithm was initially tested. To perform this method, the maps of mean curvature for each sheet were first thresholded to produce a binary image separating creases, a pixel value of 1, and background, a pixel value of 0. The distance of each background pixel to its nearest crease pixel was then computed. The negative distance, which can be regarded as a topographic surface of hills and basins, was used as the elevation map for the watershed algorithm. The basins correspond to local minima of the surface, and regions centered at each basin are flooded until all pixels are assigned a basin. We identified pixels belonging to the same basin of the elevation map as a single facet of the crease pattern. However, several concerns prompted a more careful labeling by hand. Firstly, the separation of creases from background was performed using a custom technique referred to as the Radon transform method, detailed in the Supplementary Discussion of Ref. [1]. This technique combines global and local thresholding to accommodate variations in the intensity (curvature) of creases; nevertheless, softening of old creases near strongly imprinted ones weakens their detection. The watershed algorithm proved sensitive to creases which scar the sheet but do not form closed contours, particularly true at low confinement (high∆). Thus, the algorithm over-partitions the crease network in these cases. Mitigating the effect of smaller, isolated creases and vertices by stricter thresholding also compromises the detection of smaller facets, impacting densely scarred samples at low∆. The results of watershed segmentation and corresponding scaled facet area distributions are presented in Supplementary Figs. 3 & 4; while there is consistency with the manually segmented data, lower resolution and weaker performance for small features impacts the range over which consistent scaling is observed. Nevertheless, the automated method allows us to more easily segment a larger number of crease patterns; thus, we process a more extensive set of crumpling iterations for each compaction ratio considered.

SUPPLEMENTARY NOTE 1
Scaling solution to the fragmentation rate equation. Facet fragmentation is modeled following the theory of fragmentation kinetics outlined in Supplementary Ref. [2]. Here we reproduce the derivation of a scaling solution shared among setups with similar families of breakup rates, as well as carry out the analytical steps unique to our specific choice of these rates. The linear integro-differential equation describing the evolution of concentration of facet areas x, c(x, t), is given by: where t measure of progression of fragmentation r(x) overall rate at which a facet of area x breaks f (x|y) conditional probability that x is produced from the breakup of y, y ≥ With this formulation, larger facets are more likely to split due to the higher rate given by r(x) assuming λ > 0. Moreover, the conditional probability must satisfy area conservation, y 0 xf (x|y)dx = y.
Plugging the scaling ansatz and general homogeneous kernels into the rate equation, and defining ξ = x/s, η = y/s yields whereṡ ≡ ds/dt. By separating the dependence on x and t we must have that = ω = constant and thus have two equations Insight from experimental facet fragmentation data reveals a suitable form for the conditional breakup rate: where ρ x y = (β + 1) x y β with β a free parameter and ρ(x/y) the probability density function of facet areas x normalized by their parent facet's area y from the previous crumpling iteration. In other words, ρ(x/y)d(x/y) is the probability that a facet breaks to produce a fragment that is x/y of its initial area. This formulation introduces the assumption that fragmentation is a scale invariant process. Next we demonstrate the agreement of our scaling function with the rate equation. We begin by re-expressing β in terms of a new parameter a as Our proposed solution φ(ξ) takes the form and thus Supplementary Equation (2a) may be solved to obtain which is solved for all ξ when Moving to Supplementary Equation (2b), we therefore have thaṫ noting that λ > 0 in the case considered. The initial condition s 0 = 1 may be substituted for a fragmentation process originating from a single facet. However, in the limit of large t, the dependence of the typical area becomes insensitive to the initial condition, and our result may be simplified to s(t) = G(a, λ) t −1/λ .
For the special case of λ = 1/2 considered in the main results of this work, Displacement of a 1-D random walk. For a random walk in one dimension comprised of random displacements R i , the displacement from the origin after 2k steps is If consecutive steps occur in opposite directions, representing a fold, then individual steps can be grouped into k right (positive) and k left (negative) steps: i=2,4,6...

|R i |
for each r i drawn from the same distribution. If segment lengths |R i | are drawn from a gamma distribution with shape parameter a + 1 and scale parameter θ, consistent with the distribution of facet lengths traversed by a one-dimensional vertical cross-section, then R + k and R − k are each distributed according to a gamma distribution with shape parameter k(a + 1) and scale parameter θ: + 1), θ Thus D 2k is the difference of two identically distributed gamma variates. We can obtain the probability density function for D 2k through a convolution of the probability density functions of R + k and −R − k . Let X = R + k , Y = R − k , and Z = D 2k ; then As f X (x) and f Y (y) both have non-negative support, where we have chosen the integration variable in the convolution to ensure the arguments of the probability density functions remain positive. With identical gamma distributions The integral above may be solved using the following identity [3]: where K ν (z) is the modified Bessel function of the second kind of order ν. This gives f Z (z) should be a valid probability density function, and we can verify it indeed integrates to 1 over its support z ∈ [0, ∞) using the following identity [4]: By symmetry about z = 0 we can integrate the following: as Γ(1/2) = √ π. Furthermore, for gamma-distributed steps, the average segment length is given by (a + 1)θ. Thus, in our system of a one-dimensional folded strip of total length L 0 , k and θ are related as k = L 0 2(a + 1)θ for a strip folded into 2k segments. Note that the total length of the walk is distributed as By a change of variables u = w t/2L 0 , which to leading order in u yields by consequence of Supplementary Equation (8), where c is an integration constant. In order to solve for u, we recall the definition of the Lambert W function, or product logarithm, which gives the inverse solution where W 0 (y) is the principal branch of the Lambert W function valid for real x and y, and positive y. Defining a new Making a final asymptotic approximation, W 0 (y) ≈ log(y), we thus have that With the condition t(n = 0) = 0, we obtain the final relation t(n,∆; w) = 2L 0 w 2 log 1 + SUPPLEMENTARY NOTE 4 Critical confinement. Our model of a folded one-dimensional strip as a random walk relates geometric incompatibility to the random walk stepping outside a confinement distance w. This critical distance w is dictated by the geometry of the imposed confinement, and the way in which the one-dimensional strip folds into stacks of one or more folded layers to accommodate its full length within the allowed space. Let m represent the number of spaced stacks, and p the average number of layers per stack, in our strip of length L 0 , confined to a rectangular container of width R and height L, L ≤ L 0 . To satisfy the constraint of total length L 0 at any compaction∆ = L/L 0 , we must have mp (L/m) 2 + R 2 = L 0 . The critical width w of facets which would fragment under further confinement is given by w = L 0 /mp = (L/m) 2 + R 2 . At low confinement, p ≈ 1, and thus our constraint gives m = (L 0 /R) 1 −∆ 2 , resulting in At high confinement, the collapse of stacks leads to a decrease in m that scales in proportion to L, in turn scaling the number of layers p ∼ 1/L. Specifically, we can define which is consistent with our result at low confinement taken to the limit of small L. Thus we use Supplementary Equation (10) throughout. Substituting Supplementary Equation (10) into Supplementary Equation (9), we arrive at an expression for t solely in terms of n and∆: t(n,∆) =c 1 (1 −∆ 2 ) log 1 +c 2 ñ ∆(1 +∆) , wherec 1 = 2L 0 /R 2 andc 2 = αR 2 /L 0 √ 2π.