## Abstract

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.

## Introduction

Crumpling is a complex, non-equilibrium process arising in diverse systems across a wide range of length scales, from the microscopic crumpling of graphene membranes^{1}, to the macroscopic folding of Earth’s viscoelastic crust^{2}. Crumpled structures are highly porous, providing function for applications such as high-performance batteries and supercapacitors by increasing the electrochemical surface area^{3,4}. Controlled crumpling has also been used to tune electronic, optical, and surface properties in graphene films^{1}. Further, understanding the mechanics of crumpling is essential as flexibility and shape conformation become integral considerations in the design of thin, wearable devices^{5,6}. Despite its ubiquity, a complete understanding of crumpling dynamics remains elusive due to the complexity of the disordered process. Nevertheless, statistical properties of crumpled geometries are highly reproducible in experiment^{7,8,9,10,11,12,13} and confirmed via simulation^{14,15}, which suggests that this complex process is strongly dictated by universal aspects of thin sheets such as topology and self-avoidance^{13}.

Similarly adopting a coarse-grained perspective, Gottesman et al.^{16} revealed an unexpected order to ridge network evolution in crumpled sheets. By performing a protocol of repeated compaction and unfolding, as in the schematic of Fig. 1a, they demonstrated that the intricate details of ridge networks in crumpled sheets could be reduced to a single collective quantity, the total crease length, which evolves robustly as a logarithm in the number of crumpling repetitions across varying degrees of compaction. Notably, the incremental damage added upon re-crumpling the sheet was found to be independent of the sheets’ crumpling history -the sequence of prior compactions performed to produce the current crease network. Rather, the added crease length is determined only by the current total crease length and the new compaction depth. While processes that evolve logarithmically in time are readily observed in a variety of disordered physical systems, including stress or strain relaxation of a compacted sheet^{17,18,19}, conductance relaxation of disordered electronic systems^{20}, and creep dynamics of granular suspensions^{21}, the emergence of a logarithmic model in the specific context of damage evolution in crumpled sheets is clearly distinct, and has had limited physical justification thus far.

In this work, we take a novel approach to characterize crumpling and offer explanation for the logarithmic model by drawing a correspondence between crumpling and fragmentation processes. Fragmentation models have a rich history of theoretical development^{22,23,24} as well as industrial applications^{22,25} and use in modeling collision and fracture phenomena^{26}. Here we concentrate on a theoretical, physically based rate equation for modeling time-dependent fragmentation detailed by Cheng and Redner^{27}, which provides a general framework for processes that may be treated as successive, homogeneous breakups instigated by non-local stresses. The model has been flexibly applied to describe polymer degradation^{28} and volcanic fragments expelled in an eruption^{29}, for example, though to the best of our knowledge this is the first application of such concepts to describe crumpling.

Our work is organized as follows: we derive a scaling solution to the rate equation presented in Cheng and Redner^{27} which decomposes into a time-invariant distribution of scaled facet area and a time-dependent evolution of mean facet size. We demonstrate that the derived area distribution effectively reproduces key statistical features of experimental crumpled patterns. Fragment distributions are a natural point of comparison between theory and experiment; however, in this work we go a step further to draw additional correspondence in the temporal evolution of the patterns. The temporal parameter that chronicles the evolution of mean facet size serves as an intrinsic clock measuring the maturity of the fragmentation process. We connect this to experimental parameters driving fragmentation forward, namely the number of crumpling iterations and compaction strength which characterize the experiments of Gottesman et al.^{16}. To do so, we construct a simple geometric model that likens crumpling to a random walk and is informed by the statistical properties of the derived area distribution. We derive an analytical relation for how geometric frustration occurring in a confined random walk instigates new damage and advances the temporal measure of fragmentation maturity. We demonstrate how this approach allows one to recover the logarithmic evolution of damage in ridge networks observed in Gottesman et al.^{16} and explain the history independence of damage formation, thereby furnishing a missing physical basis for this unexpected result.

The key idea behind our model is the extension of fragmentation theory to incorporate a feedback loop: as facets become smaller, they make the sheet more compliant and therefore lower the rate of subsequent fragmentation. This idea may extend to many physical systems where the accumulation of damage inhibits further damage from occurring. Our work therefore shows how fragmentation theory could be applied more generally, and suggests that the universal damage evolution seen in crumpling may have analogs in other physical systems.

## Results

The collection and processing of experimental crumpling data used to verify analytical results presented in this work is fully detailed in the Methods section. Crease patterns obtained from uniaxially compressed Mylar sheets as shown in Fig. 1b are carefully segmented into individual facets as in Fig. 1c. The samples collected vary in compaction ratio \(\tilde{{{\Delta }}}\), the ratio of final to initial height, and in the number of successive crumples of the same sheet, *n*. A total of 24 segmented crease patterns is analyzed spanning 7 different compaction ratios and including *n* = 1, 2, 3, and 24 crumpling iterations.

Throughout this work, we will refer to fragmentation in the context of crumpling as the successive partitioning of a thin sheet into smaller, flat facets separated by ridges or creases. To facilitate the construction of a model for this process, we begin with the general theory of fragmentation kinetics outlined in Cheng and Redner^{27}. Let *x* represent facet area and *c*(*x*, *t*) the concentration of facets of area *x* at time *t*; then the linear integro-differential equation describing the evolution of *c*(*x*, *t*) is given by

where the effective time *t* measures the progress or maturity of the fragmentation process, *r*(*x*) is the overall rate at which a facet of area *x* fragments, and *f*(*x*∣*y*) is the conditional probability that *x* is produced from the breakup of *y*, with *y* ≥ *x*. Inferred from this formulation are the assumptions that fragmentation occurs via a homogeneously applied external force, and independently of a facet’s shape.

### Breakup rates

In order to assess the correspondence between crumpling and a fragmentation process as described by Eq. (1), two relationships must be specified: the overall breakup rate *r*(*x*), and conditional breakup probability *f*(*x*∣*y*), which characterize fragmentation at the scale of an individual facet. Two principles help shape our formulations of the two: firstly, a common choice of *r*(*x*) consistent with physical breakup processes is the homogeneous kernel *r*(*x*) = *x*^{λ}^{27}. Furthermore, the conditional probability *f*(*x*∣*y*) must satisfy area conservation:

We use the collection of facets within each sheet as representative samples from which breakup rates may be determined. Figure 2a–c shows a typical example over three crumpling repetitions and traces the progressive fragmentation of selected facets. From such sequences, we estimate *r*(*x*) by determining the fraction of facets which fragment between two successive crumples as a function of their area *x*. The rates are computed separately for each sheet to ensure the same change in *t* elapses for all facets considered at a time. Without loss of generality, the values of *x* in all results are scaled so that 10 cm × 10 cm, the size of one sheet, corresponds to unit area. A breakup rate of the form *r*(*x*) = *x*^{λ} appears consistent with experimental breakup data, as shown in Fig. 2d. Results for samples at other compaction ratios \(\tilde{{{\Delta }}}\) are provided in Supplementary Fig. 6. We note one limitation of this analysis: sheets crumpled at a low compaction ratio may have too few facets for a robust sample size from which to infer a strong statistical trend; in the opposing extreme, sheets at high compaction likely undergo a cascade of multiple fragmentation events in a single crumpling iteration, and thus obscure the statistics of single breakup events. The power law relationship is motivated both by its consistency predominantly at low compaction, as well as the simplicity it affords later in our model.

To deduce *f*(*x*∣*y*), it is helpful to first examine the distribution *ρ*(*x*/*y*) of the area fraction *x*/*y* that a child facet occupies relative to its parent. That is, if *x* is the area of a facet at crumpling iteration *n*, and *y* the area of its enclosing facet at iteration *n*−1, then *ρ*(*x*/*y*)d(*x*/*y*) is the probability that a facet breaks to produce a fragment that is between *x*/*y* and *x*/*y* + d(*x*/*y*) of its initial area, for a small differential element d(*x*/*y*). To account for minor misalignment between successive scans, a child facet is identified if at least half of its area lies within the contour of the candidate parent facet. The area fractions display a power law distribution, as shown in Fig. 2e, and suggest a fit to a probability density function of the form

supported on *x*/*y* ∈ [0, 1]. This formulation introduces the assumption that fragmentation is a scale invariant process; while this is consistent with the present data, we note that a physical lower limit on facet area exists, and would expect deviation from scale invariant behavior as facet areas become comparable to the sheet thickness. Nevertheless, we observe clear indication of a power law relationship within our data, as shown in Fig. 2e. Extended results are provided in Supplementary Fig. 7; as previously noted, samples at high compaction undergo a succession of fragmentation events between crumples, and thus their distributions begin to depart from the power law dependence toward the more mature facet distributions of repeatedly crumpled sheets we present later, in our discussion of ridge length statistics. Taking *f*(*x*∣*y*) proportional to *ρ*(*x*/*y*) and obtaining the appropriate normalization which satisfies Eq. (2), we arrive at our final forms for the breakup rates:

It will prove useful to express the free parameter *β* as

With this definition, we demonstrate in a following subsection that the new free parameter *a* corresponds to the shape parameter for the distribution of crease length.

### Scaling solution

With *r*(*x*) and *f*(*x*∣*y*) specified, we pursue an analytical solution to the fragmentation rate equation, Eq. (1). Specifically, we seek a scaling solution independent of initial conditions, a property that allows us to solve analytically and proves compatible with the chosen form of homogeneous breakup kernels^{27}. We thereby test a scaling ansatz *c*(*x*, *t*) = *ϕ*(*ξ*)/*s*(*t*)^{2} as proposed in Cheng and Redner^{27}, where *ξ* = *x*/*s*(*t*), and the mean area, *s*(*t*), carries all explicit dependence on *t*. The scaling function *ϕ*(*ξ*) satisfies \(\mathop{\int}\nolimits_{0}^{\infty }\phi (\xi ) {\mathrm{d}}\xi =1\) and \(\mathop{\int}\nolimits_{0}^{\infty }\xi \phi (\xi ) {\mathrm{d}}\xi =1\) such that \(\mathop{\int}\nolimits_{0}^{\infty }c(x,t){\mathrm{d}}x=1/s(t)\) gives the average number of fragments, and \(\mathop{\int}\nolimits_{0}^{\infty }xc(x,t){\mathrm{d}}x=1\) is the total area, conserved by construction. We note that *ϕ*(*ξ*) is a valid probability density function and represents the distribution of the scaled facet area *ξ*. The rate equation may be solved following the procedure in Cheng and Redner^{27} as detailed in Supplementary Note 1; by this approach we arrive at a solution *c*(*x*, *t*) = *ϕ*(*ξ*)/*s*(*t*)^{2}, valid at large *t*, with

where \(G(a,\lambda )={{\Gamma }}\left(\right.\frac{a+2}{2\lambda }\left)\right./{{\Gamma }}\left(\right.\frac{a}{2\lambda }\left)\right.\), and Γ(*z*) is the gamma function. We motivate a fixed choice of the breakup rate parameter *λ* = 1/2 both by its consistency with breakup statistics at low compaction, which more accurately reflect single breakup events as discussed earlier, as well as the simplification it provides to obtain an analytically tractable model. We thereby obtain the final forms

### Ridge length statistics

To facilitate comparison with *ϕ*(*ξ*) in Eq. (7a), the area of individual facets is scaled by the mean area for that sheet and plotted as a histogram using logarithmically spaced bins. Figure 3 shows the mean curvature, hand segmentation, and scaled area distributions for a typical example from our dataset at four different crumpling iterations *n*. By our preliminary observations from Fig. 2e, we notice sample-to-sample variation in the parameter *β* (correspondingly *a*), which suggests *a* is a function of *t*; however, we expect weak dependence on *t* such that \({\mathrm{lim}\,}_{t\to \infty }{\mathrm{d}}a/{\mathrm{d}}t=0\), to uphold the assumptions made in solving Eq. (1). Indeed, an individual fit of *a* to each distribution of facet areas reveals a dependence of the form

with a universal parameter *τ*, as shown in Fig. 4a. Thus, Fig. 3c additionally shows a best fit curve to Eq. (7a) with *τ* ≈ 24.041 as a universal fitting parameter across all samples, and individual *a* and *t* for each sample computed by solving Eqs. (7b) and (8) self-consistently. The complete set of segmented crease patterns and fitted area distributions for all data samples is provided in Supplementary Figs. 1 and 2.

The close correspondence between Eq. (7a) and experimental data supports the hypothesis that successive partitioning of the sheet’s surface into facets during crumpling evolves according to the fragmentation process described by Eq. (1). We can study the further implications of this statistical description on attributes such as the distribution of crease length, which has been explored in previous studies^{7,8,9,10,14,15,30}. Let *X* be the random variable representing the area of a single facet. Following from Eq. (7a), *X* is distributed as

with \(\theta =\sqrt{s/a(a+1)}=1/t\) by consequence of Eq. (7b), and mean area *s*. Let *Y* be a random variable representing the edge length of a facet in the ridge network. If *Y* scales as \(\sqrt{X}\), the consequent distribution of *Y* is a gamma distribution,

with *θ* the scale and *a* the shape parameter, alluded to in our discussion of breakup rates, and with mean edge length *a**θ*. The distributions of facet area and edge length provided by Eqs. (9) and (10) allow us to formulate an expression for the typical total crease length as a function of *t*, in tandem with the evolution of mean area *s*(*t*). First, we briefly restate the key empirical result of Gottesman et al.^{16} to which we will compare our model. The total crease length *ℓ* was found to vary according to a logarithm of the number of crumpling and unfolding repetitions *n*:

with \(\tilde{{{\Delta }}}\) the compaction ratio, and *c*_{1} and *c*_{2} fitting parameters. A striking property of this model is its implication that the rate at which new damage accumulates, as measured by added crease length per crumpling iteration *δ**ℓ*_{empir.} ≡ ∂*ℓ*/∂*n*, is independent of the details of the sheet’s preparation:

We observe from Eq. (12) that the added crease length *δ**ℓ*_{empir.} is uniquely determined by a sheet’s instantaneous state \((\ell ,\tilde{{{\Delta }}})\); moreover, the model is independent of the details of the crease network, such as the spatial homogeneity of damage across the sheet. The fitting parameters *c*_{1} and *c*_{2} are universal to all values of *n* and \(\tilde{{{\Delta }}}\). The facet segmentation of each crease pattern provides a second measurable quantity, *d*, equal to the sum of all interior perimeters of facets; i.e., the total length of all edges shared between two facets. We expect *d* and *ℓ* to be proportional, with differences arising due to incomplete scarring around facet perimeters as regions of the sheet restore elastically, particularly for mild compression. We find that \((1-\tilde{{{\Delta }}})d\) accomplishes the desired proportionality, and define \({\ell }_{\text{meas.}}\equiv (1-\tilde{{{\Delta }}}){d}_{\text{meas.}}\) to be the measured total crease length obtained from our segmented data. Next, working with the moments of our derived facet area and edge length distributions, we can estimate *d* analytically as the average length of an edge, *a**θ*, times the average number of edges. The latter may be expressed as the average number of facets, or the total sheet area divided by the typical facet area *s*, times the number of edges per facet *n*_{e}, halved to account for shared edges, which yields

Thus, we obtain that

Here, the superscript (*t*) denotes the explicit dependence of *ℓ*_{model} on *t*; in a following subsection, we develop a connection between *t* and *n* that allows *ℓ*_{model} to be expressed in terms of *n* and \(\tilde{{{\Delta }}}\), mirroring Eq. (11). A fit of Eq. (14) to the measured length *ℓ*_{meas.} reveals a value of *n*_{e} ≈ 4.2 as shown in Fig. 4b, which suggests an average of 4–5 sides per facet. Finally, Fig. 4c demonstrates the agreement between \({\ell }_{\,\text{model}\,}^{(t)}\) of Eq. (14), and *ℓ*_{empir.} of Eq. (11).

### Numerical evidence for the insensitivity to initial preparation

Now that the connection between the statistical model of facet area and total crease length has been presented, we briefly note on the insight that may be gained by additionally solving Eq. (1) numerically. A numerical integration scheme is implemented using second-order composite trapezoid rule for discretization in *x*, and second-order implicit multi-step discretization in *t*. The sample numerical result in Fig. 5 reveals a rapid convergence to the steady state analytical solution given by Eqs. (7a) and (7b), and thereby relative insensitivity to the initial state. To demonstrate the significance of this behavior, we reiterate the observed history independence of total crease length. As discussed in Gottesman et al.^{16}, sheets with different loading histories—one hand-crumpled and another deliberately folded along straight lines—yet nearly equal total crease lengths exhibited the same subsequent accumulation of damage when subjected to the protocol of Fig. 1a. Such sheets had clearly distinct initial facet area distributions: The facet areas of the deliberately folded sheet were sharply peaked near two different values, while those of the hand-crumpled sheet were broadly distributed. Thus, signatures of initial preparation appear to be quickly eclipsed by the strong attractor of the crumpled state, echoed in the rapid convergence to steady state seen numerically.

Thus far, we have established that an estimate of total crease length constructed from moments of the derived facet area and ridge length distributions shows consistency with the logarithmic scaling of Eq. (11). In the following section, we propose a simple mechanism for how the geometric incompatibility of a folded sheet and its confinement leads to further fragmentation, driving *t* forward. This argument establishes the evolution of *t* in accordance with *n*, and thus supplies the missing link to a physically based model that corroborates experimental findings.

### One-dimensional model

To offer an explanation for the observed logarithmic scaling, we develop a simple one-dimensional model that proposes how additional fragments may form when a crumpled sheet is re-crumpled, relying on the statistical descriptions of facet area and segment length formulated in the previous sections. Our goal can be summarized by the following two questions: (1) Given its current state and prescribed confinement, with what probability does a sheet undergo further fragmentation? (2) How does this probability relate to the continuous variables in the fragmentation model of Eq. (1)? First, we appeal to the axial symmetry of our confinement to simplify our view of crumpling to a 1D strip of length *L*_{0}, as shown in Fig. 6. The strip is characterized by a sequence of folds in alternating directions which divide the strip into random segments. The lengths *r* of the segments, which are equal to the cross-sections of the intercepted facets, are distributed according to the derived gamma distribution of Eq. (10), weighted by the horizontal facet width, which increases the likelihood of a facet’s occurrence within a randomly selected vertical strip. For facets of ~1:1 aspect ratio, the distribution of segment length is thereby

with the average segment length given by (*a* + 1)*θ*. A comparison of Eq. (15) with experimental data is provided for a strongly compacted sample in Fig. 6a, b, with extended results for all samples presented in Supplementary Fig. 5.

As a preliminary step, we derive the final displacement of the strip when folded at each break, in the absence of confinement. This problem can be mapped to the displacement of a walker performing a one-dimensional random walk with gamma-distributed steps. To enforce the concept of folding, the walker’s steps occur in alternating directions. The distribution *f*_{Z}(*z*) of position *Z* after 2*k* steps accurate for all *k* is derived in full in Supplementary Note 2. However, the salient trends may be likewise observed by applying the central limit theorem and considering the position *Z* valid for large *k*, or small step size, which gives

and describes a normal distribution of zero mean and variance *L*_{0}*θ*.

If a confinement is now introduced at the locations ∣*z*∣ = *w*, we next ask with what likelihood the walker steps beyond this confinement. One approach to approximate this probability is to integrate Eq. (16) for all ∣*z*∣ > *w*, producing a two-sided survival function of Eq. (16). Although this is not equivalent to our initial question, as intermediate steps may also have landed past ∣*z*∣ > *w*, it proves an acceptable estimate as the last step has the greatest variance. A more accurate calculation would be to evaluate the likelihood that a given walk escapes the confinement at any step; however, looking at the last step is useful for its simplicity in analytical form, and still captures the anticipated behavior. A comparison to the more accurate formulation is made numerically and provided in Supplementary Fig. 9. Once again, we pursue here the simpler form of the survival function valid for large *k*, and refer to Supplementary Note 2 for the exact derivation valid at all *k*. The survival function of Eq. (16), *S*_{Z}(*w*; *θ*) = *P*(∣*Z*∣ > *w*; *w*≥0), for a threshold confinement *w*, is given by

where erf(*z*) is the error function. In order for walkers at ∣*z*∣ > *w* to be restored within the limits of confinement, one or more of their steps must fragment, thereby increasing the number of steps taken and decreasing the overall average, which drives the evolution of fragmentation. This articulates our key claim: considering our original, cylindrically shaped sheets as a statistical ensemble of one-dimensional random walks, we suggest that the progression of fragmentation measured by a change *d**t*, over a single crumpling iteration *d**n*, should be proportional to the fraction of walks in the ensemble which leave the confinement imposed at ∣*z*∣ = *w*: d*t*/d*n* ~ *S*_{Z}(*w*; *θ*). Equivalently, this is the likelihood that a single random walk leaves the critical confinement. We note that this resulting fragmentation rate describes an average fragmentation likelihood given only a confinement *w* and current temporal parameter *t* = 1/*θ* describing the maturity of the fragmentation process thus far; it does not enforce direct correlations between successive crumpling iterations, whereby new folds should occur preferentially along previous ones. Instead, the decrease in fragmentation rate with *n* is encoded through the decreasing mean facet area with *t*. Moreover, while stronger correlation is expected between walks representing nearby transects of the sheet, here we consider the statistical behavior of the sheet as a whole, and account for the increased fragmentation likelihood for facets with larger horizontal extent through the weighting introduced in Eq. (15). At present, Eq. (17) gives the likelihood that new creases will form; however, it does not yet describe how much new damage is created, for which two additional factors should be considered: (1) When the sheet is strongly confined in closely packed layers, the layers tend to collectively fragment, as alluded to by Sultan and Boudaoud^{15} and Gottesman et al.^{16}, thus contributing a factor *p* ~ 1/*L* such that halving the final height doubles the number of additional ridges. (2) In the opposite limit of low compaction, facets are not in close proximity and need not behave cooperatively; thus, new damage scales linearly with the amount of compression *L*_{0}−*L*, as argued in Gottesman et al.^{16}. With these additional considerations, we propose that the evolution of the fragmentation process with crumpling iteration behaves as

where *α* is a fitted constant of proportionality. We indicate the explicit dependence on *t* here, as *t* and *θ* are inversely related. The critical width *w* is determined by the geometry of the imposed confinement, as illustrated in Fig. 6d; a complete derivation is provided in Supplementary Note 4:

where *R* is the radius of the container. By consequence of Eq. (14) we can directly relate Eqs. (18) and (12) as

and obtain a fit to the proportionality constant *α*. By performing an asymptotic approximation in the limit of large *t*, detailed in Supplementary Note 3, Eq. (18) may be analytically integrated to provide a scaling relation \(t(n,\tilde{{{\Delta }}})\) which bears similarity to \(\ell (n,\tilde{{{\Delta }}})\) of Eq. (11):

where \({\tilde{c}}_{1}=2{L}_{0}/{R}^{2}\), \({\tilde{c}}_{2}=\alpha {R}^{2}/{L}_{0}\sqrt{2\pi }\), and *L*_{0} and *R* are the sheet length (equivalently the confining container height) and container radius, respectively. Taken together, Eqs. (14) and (21) thereby provide a theoretically motivated expression \({\ell }_{\,\text{model}\,}^{(n)}\equiv \ell (t(n,\tilde{{{\Delta }}}),\tilde{{{\Delta }}})\) based on properties of fragmentation kinetics and a simple mechanism for re-fragmentation formulated as a random walk. Figure 7 compares the agreement of the empirical relations *δ**ℓ*_{empir.} and *ℓ*_{empir.}, as well as the derived models *δ**ℓ*_{model} and \({\ell }_{\,\text{model}\,}^{(n)}\), with the measured quantities \(\delta {\ell }_{\text{meas.}}\equiv {\ell }_{\text{meas.}\,}^{(n)}-{\ell }_{\,\text{meas.}\,}^{(n-1)}\) and *ℓ*_{meas.} for various *n*. Collectively, the results of Figs. 4 and 7 demonstrate clear consistency of the fragmentation model with the anticipated logarithmic growth.

## Discussion

By pursuing a correspondence between the crumpling of a thin sheet and a general fragmentation process, we have derived a physically based framework for the evolution of statistical properties of intricate crumpled patterns. Equipped with theoretical models in close agreement with experimental data, we have proposed a simple model of one-dimensional folding in which further fragmentation ensues due to a geometric incompatibility between the sequence of folds and the imposed confinement, likened to a random walk exceeding a critical allowed displacement. The predicted accrual of damage, quantified by added crease length, shows strong consistency with the logarithmic model of Gottesman et al.^{16}, and thereby supplies a possible physical basis for the puzzling origin of logarithmic scaling in repeated crumpling experiments. Furthermore, our model explains the history independence of the logarithmic scaling, since the area distribution of the crumpled state is such a strong attractor in the fragmentation process.

The consistency of crumpling with fragmentation theory hints at the possibility of universal behavior uniting more diverse fragmenting systems. For example, the activation of defects in the fragmentation of ceramics can locally slow down subsequent fracture, and may bear semblance to the slowing of damage accumulation as a re-crumpled sheet exploits its existing folds^{31}. Thus, studies of crumpled systems might offer a new lens through which to interpret other complex processes. An immediate extension of this work would be a validation of the results on sheets of varied thicknesses and material parameters, as well as those prepared according to different compaction protocols. One simplifying assumption of our analysis is that fragmentation of facets is a scale invariant process over the range of areas considered; however, this assumption starts to break down particularly for large crumpling iterations *n*. The work of Lechenault et al.^{32} offers a compelling approach for identifying this limit by considering the energetic competition between bending of facets and rigid folding along existing creases, with energy cost of the former proportional to the sheet’s bending rigidity, and the latter proportional to crease stiffness. The energy balance of these competing deformations provides a characteristic length scale which varies in proportion to the sheet thickness. This improvement to the current work would strongly benefit from further studies over a range of material parameters. Length scales of folds in crumpled systems have also been studied in the context of thermally driven dynamics, and it may thus be useful to draw possible connections to statistical mechanical models of crumpling^{33,34,35}. Moreover, it may be of value to explore slight generalizations of proposed functional forms introduced in this study, such as the breakup rates; this could allow variations across other experimental results to be explained, such as those arising between low and high compaction regimes^{12,15}, thereby providing a unifying framework for such observations.

Additionally, deeper understanding of crumpling dynamics can assist data-driven approaches to predicting damage network formation. Though machine learning methods are capable of unveiling hidden structure in complex, disordered systems^{36,37}, prior work has demonstrated the importance of preserving physical properties in making faithful predictions: for example, preserving vertex angle constraints in synthetic fold patterns to assist the task of ridge network reconstruction in crumpled sheets^{38}. In addition to encoding physical rules implicitly through data, future machine learning approaches may explicitly enforce constraints such as facet area and crease length statistics in predicting ridge network evolution. Strategies which couple detailed spatial data with coarse-grained theoretical insight could thus enable more comprehensive predictions of crumpling dynamics in future studies.

## Methods

The data analyzed in this work was collected for the study of Gottesman et al.^{16}; here, we briefly summarize the experimental protocol for reference. In total, 10 cm × 10 cm Mylar sheets are rolled into a 3-cm diameter cylindrical container and compressed uniaxially to a specified compaction ratio \(\tilde{{{\Delta }}}=L/{L}_{0}\), the ratio of final to initial height, with *L*_{0} = 10 cm (Fig. 1a). The resulting ridge network inscribed on each sheet is extracted by carefully unfolding and scanning the sheet using a custom laser profilometer, which produces a height map of the sheet. A two-dimensional map of mean curvature is determined from the spatial gradients of the height profile; sharp peaks in curvature mark the signature of a ridge (Fig. 1b). Successive re-crumpling and scanning of a single sheet is performed *n* times up to *n* = 24. Individual facets, characterized as contiguous regions of near-zero curvature, are delineated as shown in Fig. 1c. Due to noise and artifacts in data collection, not all facets are completely enclosed by a contour of ridges; breaks along a ridge, or smoothing out and softening of ridges, occur inevitably during re-crumpling and unfolding. Automated methods of crease detection and facet labeling were initially tested to perform the segmentation; however, these methods proved sensitive to noise and thus were prone to over-fragmenting the sheets. Each sample presented and analyzed in this work was digitally labeled by hand. Additional details of automated segmentation are provided in the Supplementary Methods, and resulting segmentations and facet area distributions are shown in Supplementary Figs. 3 and 4. With manual segmentation, care was taken to identify not only the dominant lines of each pattern as seen in the examples, but also the less pronounced softer scarring. The segmentation was performed for sheets after iterations *n* = 1, 2, and 3 at seven different compaction ratios: \(\tilde{{{\Delta }}}=0.63,0.45,0.36,0.27,0.18,0.09\), and 0.045. Each series of successive crumples was compared across all iterations *n* for consistency, to ensure that labeled facets from earlier iterations persist in later ones. Samples with \(\tilde{{{\Delta }}}=0.63,0.45\), and 0.27 were also labeled after *n* = 24 crumples, for a total of 24 samples overall. We acknowledge that samples at *n* = 24 are more prone to missing detail as older scarring is obscured by newer ridges, but are nonetheless valuable to the study. The results of manual segmentation and corresponding facet area distributions are provided in Supplementary Figs. 1 and 2.

## Data availability

No new experimental data was produced for this study; all data analyzed was previously collected and reported in Gottesman et al.^{16}. The subset of data from Gottesman et al.^{16} used in this article is provided in post-processed form with our analysis codes.

## Code availability

The data processing and analysis codes are available on GitHub at https://github.com/jandrejevic12/fragmentation_model^{39}.

## Change history

### 15 April 2021

A Correction to this paper has been published: https://doi.org/10.1038/s41467-021-22791-z

## References

- 1.
Zang, J. et al. Multifunctionality and control of the crumpling and unfolding of large-area graphene.

*Nat. Mater.***12**, 321–325 (2013). - 2.
Beloussov, V. The origin of folding in the Earth’s crust.

*J. Geophys. Res.***66**, 2241–2254 (1961). - 3.
Song, J., Yu, Z., Gordin, M. L. & Wang, D. Advanced sulfur cathode enabled by highly crumpled nitrogen-doped graphene sheets for high-energy-density lithium–sulfur batteries.

*Nano Lett.***16**, 864–870 (2016). - 4.
Wen, Z. et al. Crumpled nitrogen-doped graphene nanosheets with ultrahigh pore volume for high-performance supercapacitor.

*Adv. Mater.***24**, 5610–5616 (2012). - 5.
Kaltenbrunner, M. et al. An ultra-lightweight design for imperceptible plastic electronics.

*Nature***499**, 458–463 (2013). - 6.
White, M. S. et al. Ultrathin, highly flexible and stretchable pleds.

*Nat. Photon.***7**, 811–816 (2013). - 7.
Deboeuf, S., Katzav, E., Boudaoud, A., Bonn, D. & Adda-Beida, M. Compaction of thin sheets: crumpling and folding. in

*21st French Mechanics Congress, Bordeaux, France*https://api.semanticscholar.org/CorpusID:56431939 (2013). - 8.
Deboeuf, S., Katzav, E., Boudaoud, A., Bonn, D. & Adda-Bedia, M. Comparative study of crumpling and folding of thin sheets.

*Phys. Rev. Lett.***110**, 104301 (2013). - 9.
Andresen, C. A., Hansen, A. & Schmittbuhl, J. Ridge network in crumpled paper.

*Phys. Rev. E***76**, 026108 (2007). - 10.
Blair, D. L. & Kudrolli, A. Geometry of crumpled paper.

*Phys. Rev. Lett.***94**, 166107 (2005). - 11.
Balankin, A. S. et al. Intrinsically anomalous roughness of randomly crumpled thin sheets.

*Phys. Rev. E***74**, 061602 (2006). - 12.
Balankin, A. S. & Huerta, O. S. Entropic rigidity of a crumpling network in a randomly folded thin sheet.

*Phys. Rev. E***77**, 051124 (2008). - 13.
Balankin, A. S. et al. Fractal features of a crumpling network in randomly folded thin matter and mechanics of sheet crushing.

*Phys. Rev. E***87**, 052806 (2013). - 14.
Vliegenthart, G. & Gompper, G. Forced crumpling of self-avoiding elastic sheets.

*Nat. Mater.***5**, 216–221 (2006). - 15.
Sultan, E. & Boudaoud, A. Statistics of crumpled paper.

*Phys. Rev. Lett.***96**, 136103 (2006). - 16.
Gottesman, O., Andrejevic, J., Rycroft, C. H. & Rubinstein, S. M. A state variable for crumpled thin sheets.

*Commun. Phys.***1**, 70 (2018). - 17.
Lahini, Y., Gottesman, O., Amir, A. & Rubinstein, S. M. Nonmonotonic aging and memory retention in disordered mechanical systems.

*Phys. Rev. Lett.***118**, 085501 (2017). - 18.
Balankin, A. S., Huerta, O. S., Méndez, F. H. & Ortiz, J. P. Slow dynamics of stress and strain relaxation in randomly crumpled elasto-plastic sheets.

*Phys. Rev. E***84**, 021118 (2011). - 19.
Matan, K., Williams, R. B., Witten, T. A. & Nagel, S. R. Crumpling a thin sheet.

*Phys. Rev. Lett.***88**, 076101 (2002). - 20.
Eisenbach, A. et al. Glassy dynamics in disordered electronic systems reveal striking thermal memory effects.

*Phys. Rev. Lett.***117**, 116601 (2016). - 21.
Bérut, A., Pouliquen, O. & Forterre, Y. Brownian granular flows down heaps.

*Phys. Rev. Lett.***123**, 248005 (2019). - 22.
Rosin, P. Laws governing the fineness of powdered coal.

*J. Inst. Fuel***7**, 29–36 (1933). - 23.
Kolmogoroff, A. Über das logarithmisch normale verteilungsgesetz der dimensionen der teilchen bei zerstückelung (translated as ‘the logarithmically normal law of distribution of dimensions of particles when broken into small parts’).

*CR (Doklady) Acad. Sci. URSS (NS)*,**31**, 91–101 (1941). - 24.
Grady, D.

*Physics of Shock and Impact*, Vol. 1. 2053–2563 (IOP Publishing, 2017). - 25.
Brown, W. K. & Wohletz, K. H. Derivation of the Weibull distribution based on physical principles and its connection to the Rosin–Rammler and lognormal distributions.

*J. Appl. Phys.***78**, 2758–2763 (1995). - 26.
Brown, W. K. A theory of sequential fragmentation and its astronomical applications.

*J. Astrophys. Astron.***10**, 89–112 (1989). - 27.
Cheng, Z. & Redner, S. Kinetics of fragmentation.

*J. Phys. A Math. General***23**, 1233 (1990). - 28.
Wang, M., Smith, J. & McCoy, B. J. Continuous kinetics for thermal degradation of polymer in solution.

*AIChE J.***41**, 1521–1533 (1995). - 29.
Kaminski, E. & Jaupart, C. The size distribution of pyroclasts and the fragmentation sequence in explosive volcanic eruptions.

*J. Geophys. Res. Solid Earth***103**, 29759–29779 (1998). - 30.
Balankin, A. S. & Flores-Cano, L. Edwards’s statistical mechanics of crumpling networks in crushed self-avoiding sheets with finite bending rigidity.

*Phys. Rev. E***91**, 032109 (2015). - 31.
Levy, S. & Molinari, J.-F. Dynamic fragmentation of ceramics, signature of defects and scaling of fragment sizes.

*J. Mech. Phys. Solids***58**, 12–26 (2010). - 32.
Lechenault, F., Thiria, B. & Adda-Bedia, M. Mechanical response of a creased sheet.

*Phys. Rev. Lett.***112**, 244301 (2014). - 33.
Yllanes, D., Nelson, D. & Bowick, M. Folding pathways to crumpling in thermalized elastic frames.

*Phys. Rev. E***100**, 042112 (2019). - 34.
Bowick, M. J. & Travesset, A. The statistical mechanics of membranes.

*Phys. Rep.***344**, 255–308 (2001). - 35.
David, N., Tsvi, P. & Steven, W.

*Statistical Mechanics Of Membranes And Surfaces*(World Scientific, 2004). - 36.
Cubuk, E. D. et al. Identifying structural flow defects in disordered solids using machine-learning methods.

*Phys. Rev. Lett.***114**, 108001 (2015). - 37.
Carrasquilla, J. & Melko, R. G. Machine learning phases of matter.

*Nat. Phys.***13**, 431–434 (2017). - 38.
Hoffmann, J. et al. Machine learning in a data-limited regime: augmenting experiments with synthetic data uncovers order in crumpled sheets.

*Sci. Adv.***5**, eaau6792 (2019). - 39.
Andrejevic, J., Lee, L. M., Rubinstein, S. M. & Rycroft, C. H. A model for the fragmentation kinetics of crumpled thin sheets. Release 1.0.1 https://doi.org/10.5281/zenodo.4411510 (2020).

## Acknowledgements

This research was partially supported by the National Science Foundation through the Harvard University Materials Research Science and Engineering Center under grant nos. DMR-1420570 and DMR-2011754. J.A. acknowledges support from the National Science Foundation under grant no. DGE-1745303. C.H.R. was partially supported by the Applied Mathematics Program of the U.S. DOE Office of Science Advanced Scientific Computing Research under Contract No. DE-AC02-05CH11231.

## Author information

### Affiliations

### Contributions

S.M.R. and C.H.R. conceived the study and supervised the project. L.M.L. performed new experiments. J.A. performed the image segmentation. J.A. and C.H.R. carried out analytical derivations and numerical validation. J.A. and C.H.R. wrote the paper with input from all authors.

### Corresponding author

## Ethics declarations

### Competing interests

The authors declare no competing interests.

## Additional information

**Peer review information** *Nature Communications* thanks Jean-François Molinari, and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Publisher’s note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Supplementary information

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

## About this article

### Cite this article

Andrejevic, J., Lee, L.M., Rubinstein, S.M. *et al.* A model for the fragmentation kinetics of crumpled thin sheets.
*Nat Commun* **12, **1470 (2021). https://doi.org/10.1038/s41467-021-21625-2

Received:

Accepted:

Published:

## Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.