A three-dimensional phase diagram of growth-induced surface instabilities

A variety of fascinating morphological patterns arise on surfaces of growing, developing or aging tissues, organs and microorganism colonies. These patterns can be classified into creases, wrinkles, folds, period-doubles, ridges and delaminated-buckles according to their distinctive topographical characteristics. One universal mechanism for the pattern formation has been long believed to be the mismatch strains between biological layers with different expanding or shrinking rates, which induce mechanical instabilities. However, a general model that accounts for the formation and evolution of these various surface-instability patterns still does not exist. Here, we take biological structures at their current states as thermodynamic systems, treat each instability pattern as a thermodynamic phase, and construct a unified phase diagram that can quantitatively predict various types of growth-induced surface instabilities. We further validate the phase diagram with our experiments on surface instabilities induced by mismatch strains as well as the reported data on growth-induced instabilities in various biological systems. The predicted wavelengths and amplitudes of various instability patterns match well with our experimental data. It is expected that the unified phase diagram will not only advance the understanding of biological morphogenesis, but also significantly facilitate the design of new materials and structures by rationally harnessing surface instabilities.

A variety of fascinating morphological patterns arise on surfaces of growing, developing or aging tissues, organs and microorganism colonies. These patterns can be classified into creases, wrinkles, folds, period-doubles, ridges and delaminated-buckles according to their distinctive topographical characteristics. One universal mechanism for the pattern formation has been long believed to be the mismatch strains between biological layers with different expanding or shrinking rates, which induce mechanical instabilities. However, a general model that accounts for the formation and evolution of these various surface-instability patterns still does not exist. Here, we take biological structures at their current states as thermodynamic systems, treat each instability pattern as a thermodynamic phase, and construct a unified phase diagram that can quantitatively predict various types of growth-induced surface instabilities. We further validate the phase diagram with our experiments on surface instabilities induced by mismatch strains as well as the reported data on growth-induced instabilities in various biological systems. The predicted wavelengths and amplitudes of various instability patterns match well with our experimental data. It is expected that the unified phase diagram will not only advance the understanding of biological morphogenesis, but also significantly facilitate the design of new materials and structures by rationally harnessing surface instabilities. N umerous intriguing morphologies and phenomena on surfaces of growing animals, plants and microorganism colonies have fascinated artists and scientists for decades 1,2 . Abundant examples (Fig. 1A) can be found in various types of living creatures across multiple size scales, such as wrinkles on skins of mammalians, plants and fruits [3][4][5][6][7][8] , undulations in developing biofilms [9][10][11] , grooves on the cerebral cortex [12][13][14][15] , mucosal villi and folds of airways, esophagi and guts [16][17][18][19][20][21][22] , buckled tumor surfaces 23,24 , epithelial cell delamination due to tissue crowding 25,26 , and crumpled membranes of blood cells 27 . Although these biological patterns may be results of complex genetic, biological and biochemical processes, recent studies have suggested that growthinduced mechanical forces regulate the formation and evolution of biological patterns 2,16,18,[28][29][30] . Biological structures usually consist of multiple layers with strikingly different biochemical compositions and mechanical properties; for example, epidermis on the dermis or hypodermis of mammalian skins [3][4][5] , the epidermis on the ground tissue of plant skins [6][7][8] , biofilms on culture gels [9][10][11] , the grey matter on the white matter of cerebral cortexes 12,13 , the mucosa on the muscle layer of airways, esophagi and guts [16][17][18][19][20][21][22] , outer proliferative cells on the inner necrotic core of a tumor 23 , epithelial cell monolayer on the underlying tissue 25,26 , membranes on the cytoskeleton of blood cells 27 . During growth, development or aging, different layers of biological structures usually have different expanding or shrinking rates, thus resulting in mismatch strains between the biological layers. The surface topographical patterns have long been believed to be results of mismatch-induced compressive strains in the skin layers which have higher growth rates or lower shrinkage rates than the underlying biological layers [14][15][16]18,31 . Once the mismatch compressive strain rises to critical values, the initially flat surface of the film becomes unstable and bifurcate into different types of corrugated patterns (Fig. 1B), including (i) wrinkle -the film undulates sinusoidally but remains locally smooth (e.g., the pumpkin skin in Figs 1Ai and 1Bi) 6 , (ii) creasethe surface of the film folds into dispersed regions of self-contacts with sharp tips (e.g., the cerebral cortex in Figs 1Aii and 1Bii) 12,32 , and (iii) delaminated-buckle -the film delaminates from the substrate to form buckled regions (e.g., the biofilm in Figs 1Aiii and 1Biii) 10 . As the mismatch strain further increases, the wrinkles may further bifurcate into more complicated surface patterns, including (iv) fold -some valleys of the wrinkle fold into selfcontacts with sharp tips (e.g., the dog skin in Figs 1Aiv and 1Biv) 33 , (v) period-double -the sinusoidal wrinkle transits into a pattern with twice of the wavelength (Fig. 1Bv), and (vi) ridge -the wrinkle drastically increases its amplitude but decreases its wavelength, forming a high-aspect-ratio pattern that ceases to follow sinusoidal shape (Fig. 1Bvi). These instability patterns with distinctive topographical characteristics have been only studied and identified separately in different biological systems under varied physical and biological conditions 6,10,16,18,25 . However, a general model that can quantitatively predict the formation and evolution of various types of surface-instability patterns still does not exist; primarily because existing theories such as linear stability analysis cannot systematically analyze all modes of instabilities 12 , and existing experiments did not systematically vary the mechanical properties of film-substrate systems.
Here, we take biological film-substrate structures at their current states as thermodynamic systems, and regard each mode of surfaceinstability pattern as a thermodynamic phase. By systematically varying mechanical properties of the structure including moduli, adhesion energy and mismatch strain of the film and substrate, we calculate the initiation and evolution of various modes of growthinduced surface instabilities. We then compare potential energies of different instability patterns, and construct a quantitative phase dia-gram that accounts for all instability patterns discussed above, by assuming the current pattern seeks the lowest potential energy among all possible configurations. To validate the phase diagram, we impose different mismatch strains in polymeric film-substrate structures with systematically varied rigidity and adhesion energy to induce various modes of instability patterns. The resultant patterns indeed follow the phase diagram quantitatively. We further find that the phase diagram agrees well with reported data on growthinduced surface instabilities from a number of previous studies. It is expected that the phase diagram will not only advance the understanding of biological morphogenesis, but also significantly facilitate the design of new structures with innovative surfaces or interfaces for disease therapy 22,24 , active cell culture 34 , biofouling management 35 , tunable superhydrophobicity 36 and flexible electronics 37,38 .

Results
A three-dimensional phase diagram. While the development of instability patterns in biological structures may involve compli- cated biological processes, determining the instability patterns at current states can be solved as mechanics problems 2,16,18,[28][29][30] . To focus on essential physical features, we simplify the layered biological structures at the current states as a homogeneous film adhered on a homogeneous underlying substrate, both undergoing plane-strain deformation (Fig. 1C). To account for large deformation, both the film and the substrate are taken as incompressible neo-Hookean materials with shear modulus m f and m s , respectively. If the film and the substrate at the current state are detached from each other, they will have lengths L f and L s and thicknesses H f and H s , respectively (Fig. 1Ci). We define the mismatch strain between the film and the substrate at current state as e M 5 (L f 2 L s )/L f . Since film thickness H f is much smaller than all the other dimensions (i.e., L f , L s and H s ) in the system, it is the only relevant length scale for analyzing the instability patterns. We further define the adhesion energy between the film and the substrate, C, as the work required to detach the film from a unit area of the substrate in the stress-free state.
Within the time scale of determining instability patterns, we take the film-substrate structure as a thermodynamic system, and assume the current surface-instability pattern always seeks the lowest potential energy among all possible configurations (Fig. 1D), i.e., following the Maxwell stability criterion [39][40][41][42] . The potential energy per unit width of the film-substrate system under plane-strain deformation can be expressed as 39 where U f and U s are strain energies per unit width of the film and substrate, respectively, and D is the current delaminated length of the substrate measured in the stress-free state (Fig. 1Ci). This simplified model involves five physical parameters that determine the instability patterns: m f , m s , H f , C and e M . By dimensional argument, they can be normalized into three dimensionless parameters: modulus ratio m f /m s , normalized adhesion energy C/(m s H f ) and mismatch strain e M . The types of instability patterns will be solely determined by the three dimensionless parameters, and therefore governed by a threedimensional phase diagram. It should be noted that biological structures can take different paths to induce mismatch strains such as expansion of films or shrink of substrates (see e.g., Supplementary Figs S5 and S7); however, structures with the same set of m f /m s , C/(m s H f ) and e M should reach the same type of instability pattern at the current state, given the Maxwell stability criterion is followed. Next, we discuss the process to quantitatively construct the phase diagram. A plane-strain finite element model is developed to calculate the formation of instability patterns (Methods and SI). To induce mismatch strains in the model, we assume the detached stress-free substrate in Fig. 1Ci is pre-stretched by a ratio of L f /L s , adhered to the film (Fig. 1Cii), and then relaxed to length L (Fig. 1Ciii), during which all deformation occurs in plane-strain condition. The overall compressive strain in the film is defined as e 5 (L f 2 L)/L f (Fig. 1Ciii). As e increases to critical values, patterns of surface instabilities can initiate and transit into others (Fig. 1D). Force perturbations and mesh defects have been introduced into the model as fluctuations to facilitate the system to seek minimum-potential energy states (Fig. 1D). When the substrate is fully relaxed (i.e., L 5 L s and e 5 e M , shown as the black solid circle on Fig. 1D), the resultant pattern is the instability pattern of the film-substrate system with mismatch strain e M , which represents a point of one phase in the phase diagram (Fig. 2). The boundaries between regions of different phases give the phase boundaries on the phase diagram. We can also determine the phase boundaries by comparing the potential energies of different patterns with the same set of m f /m s , C/(m s H f ) and e M 39,43,44 , i.e., where P i and P j are the potential energies of two different patterns on film-substrate models with the same properties and dimensions (Fig. 1D). Following this method, we categorize all modes of surface instabilities patterns discussed above into a three-dimensional phase diagram with quantitatively determined phase boundaries (Fig. 2).
To understand the phase diagram, we first consider the scenario in which the adhesion between the film and the substrate is so strong that the film does not delaminate from the substrate (i.e., D 5 0). The instability patterns are thus governed only by m f /m s and e M , giving a two-dimensional phase diagram (i.e., C/(m s H f )R' on Figs 2 and 3A). When the mismatch strain e M is sufficiently low, the flat filmsubstrate structure has lower potential energy than any instability pattern. As the mismatch strain increases to critical values, the flat state will transit into either wrinkled or creased state, depending on the modulus ratio. When m f /m s , 1.3 (i.e., relatively compliant film), the film tends to fold against itself without deforming the substrate to  minimize the potential energy of the system. The phase boundary between the flat and creased states, which is calculated by setting P flat 5 P crease 12,43 , is a vertical line on Figs 2 and 3A, where e C crease is the critical mismatch strain, at which the structure transits from flat to creased state. It is noted that, for 0.5 , m f /m s , 1.3, the creases may further develop into folds under larger mismatch strains, i.e., e M . 0.45 ( Supplementary Fig. S1).
On the other hand, when m f /m s . 1.3 (i.e., relatively stiff film), the film tends to undulate together with the substrate to minimize the potential energy of the system. The phase boundary between the flat and creased states, which is calculated by setting P flat 5 P wrinkle (See SI and Supplementary Fig. S2 where e C wrinkle is the critical mismatch strain, at which the structure transits from flat to wrinkle state. Notably, the triple point between flat, creased and wrinkled states in a film-substrate structure with mismatch strain is at e M 5 0.35 and m f /m s 5 1.3. As the mismatch strain further increases, the wrinkled filmsubstrate structure can further bifurcate into more complicated pat- terns including fold, period-double and ridge, depending on the modulus ratio m f /m s . Qualitatively, the pattern of fold develops at a lower range of m f /m s (i.e., relatively compliant films) than period double and ridge, because the fold requires large deformation and self-contact of films. Quantitatively, the calculated phase boundaries between fold, period-double and ridge are approximately horizontal lines at m f /m s < 12 and m f /m s < 800, respectively, in the region of 0.3 , e M , 0.6 (Figs 2 and 3A). In addition, the calculated phase boundary between wrinkle and fold is in the region of 1.3 , m f /m s , 12 and 0.3 , e M , 0.35 (Figs 2 and 3A, Supplementary Fig. S3); between wrinkle and period double in the region of 12 , m f /m s # 800 and 0.3 , e M , 0.32 (Figs 2 and 3A, Supplementary Figs S4 and S5); and between wrinkle and ridge states in the region of 800 , m f /m s , 10 6 and e M < 0.33 (Figs 2 and 3A, Supplementary Figs S6 and S7). If the mismatch strain increases beyond 0.6, more complicated patterns will develop, such as period triple 45,46 , period quadruples 45 , and co-existence of fold, period double and ridge, which will be not be covered in the current paper. It should be noted, since the Maxwell stability criterion is followed here, the formation of these instability patterns is independent on the processes of inducing mismatch strains. Our finite-element calculations indeed give the same instability patterns, while following different pathways to induce the mismatch strains, such as substrate pre-stretching and relaxing, film swelling and substrate shrinking (Fig. 1C, Supplementary Figs S5 and S7).
Next, we discuss the scenario in which the adhesion between the film and substrate is relatively weak (i.e., relatively low C/(m s H f )) and the film debonds from the substrate to form delaminated-buckles (i.e., D . 0). To calculate the potential energy of a delaminatedbuckled pattern P buckle , we first prescribe various pairs of delaminated length D and delaminated-buckle wavelength l in the finiteelement model (See Supplementary Figs S8A, S9A and S10A). Minimization of the potential energy of the structure requires 47 By solving Eqs. 5 and 6, we can obtain a set of D, l and P buckle for the film-substrate structure, which is assumed to be in the delaminated-buckled state (Supplementary Figs S8B, S9B and S10B). To calculate the phase boundary between the delaminated-buckle and other states, we set P buckle 5 P ud for two models with the same properties and dimensions, where P ud represents the potential energy of an un-delaminated state. For fixed values of m f /m s and e M , it is evident that P buckle is a monotonically increasing function of C/(m s H f ) but P ud is a constant (Supplementary Figs S8C, S9C, and S10C). Therefore, the critical value of C/(m s H f ) on the phase boundary between the delaminated-buckled state and un-delaminated states is determined by the intersections of the two curves of P buckle and P ud (Supplementary Fig. S11).
Overall, the three-dimensional phase diagram for growth-induced surface instabilities in film-substrate structures can be understood as follows (Fig. 2). As the mismatch strains in film-substrate structures reach critical values, the initially flat structures can transit into various types of surface instability patterns. If the normalized adhesion energies between the films and substrates are sufficiently high, the preferable patterns are crease and fold for structures with relatively compliant films, but wrinkle, period double and ridge for structures with relatively stiff films. As the normalized adhesion energy decreases to critical values, the un-delaminated patterns transit into the delaminated-buckle patterns.
Wavelengths and amplitudes of growth-induced surface instabilities. We next study the wavelengths and amplitudes of instability patterns on the phase diagram (Fig. 2). For a film-substrate system with m f /m s , 1.3, a mismatch strain above 0.35 induces nucleation of scattered creases on the surface of the film; further increasing e M leads to a pattern of creases with a wavelength l crease (Supplementary Fig. S13). The wavelength of the crease pattern approximately follows a linear relation with e M as 48 which matches consistently with the current experimental data and the previous study on crease pattern in materials under compression 48 (Supplementary Fig. S13).
If the modulus ratio of the system m f /m s is above 1.3, the initially flat surface first bifurcates into the wrinkle phase. The wavelength of the wrinkles at initiation (i.e., just transited from flat state) can be calculated by linear stability analysis and expressed as a function of m f / m s (see SI and Supplementary Fig. S2D) 49 . For large modulus ratios (i.e., m f /m s . 10 3 ), the wavelength of wrinkles at initiation can be approximated as 50 Further increasing the mismatch strain will decrease the wavelengths of the wrinkles. We adopt a modified accordion model to calculate the wavelength of wrinkles with increased mismatch strain. Without loss of generality, let's consider a case that the mismatch strain is induced by the shrinkage of the substrate and the length of the film maintains constant during the process, as illustrate in Supplementary Fig. S14. At a critical mismatch strain, wrinkles with initial wavelength sets in the film-substrate structure. The modified accordion model assumes that the number of undulations in the wrinkles does not change as the substrate further shrinks to increase the mismatch strain above the critical mismatch strain. Therefore, the wavelength of the wrinkles varies according to the mismatch strains at wrinkle initiation and current state, as 37,38,45 where e C wrinkle is the critical mismatch strain for wrinkle initiation, given by Eq. 4. The modified accordion model further assumes that the number of undulations in the patterns still maintain the same when the wrinkles transit into folds, period doubles and ridges due to further shrinkage of the substrate. Therefore, we can calculate the wavelengths of folds, period doubles and ridges under various mismatch strains, respectively, as It is noted that the period-double has a wavelength twice of the corresponding wrinkles, as given in Eq. 11.
Once the wavelengths of wrinkles are obtained from Eq. 9, we can further calculate the amplitude of wrinkles by approximating the arclength of wrinkles at various mismatch strains to be equal to the wavelength of wrinkles at initiation (See Supplementary Fig. S14) 37,38 , i.e., where ()9 is differential operation for x.  Fig. S12) 39 . In order to avoid the film-substrate delamination, a very high adhesion energy (i.e., C . 10 3 Jm 22 and C/(m s H f ) . 10 3 ) is achieved by smearing a thin adhesive layer between the film and substrate 47,48,51 . Since the adhesive layer is much thinner than the film and its modulus approximates that of the substrate, the adhesive layer does not affect the instability patterns 39 . We first discuss the five modes of patterns observed in the filmsubstrate structures with strong adhesion that prevents delamination: (i) If the film is more compliant than the substrate, for example m f /m s 5 0.3 or 0.64, the structure maintains flat under relatively low mismatch strain. When e M reaches ,0.36 (for m f /m s 5 0.3 and 0.64), the initially flat surface suddenly forms discrete creases as indicated by arrows in Fig. 3B, which then evolve into periodically distributed creases with the rise of e M 12,43 . (ii) When the modulus ratio increases As the mismatch strains further increase to critical values, the wrinkles can bifurcate into folds, period doubles or ridges, depending on the modulus ratio. (iii) If the film is slightly stiffer than the substrate with m f /m s 5 1.86, 3.64 or 9.79, the wrinkled surface folds against itself to form creases at some valleys once e M reaches ,0.33 (Fig. 3C). With further increase of mismatch strain, all valleys sequentially collapse into folds and further penetrate into the substrate (e.g., e M 5 0.59 in Fig. 3C). (iv) When the modulus ratio is further increased to a higher range, with m f /m s 5 14.77, 67.24 or 130.74, the wrinkles transit into period-doubles at e M < 0.32, by growing the amplitude of one wrinkle at the expense of its neighbors (Fig. 3D). With further increasing e M , the crests of period doubles may contact each other to form channels in the valleys (see e M 5 0.48 in Fig. 3D). (v) If the film is much stiffer than the substrate, for example m f /m s 5 1.59 3 10 3 , 9.11 3 10 3 or 7.94 3 10 4 , the wrinkles bifurcate into ridge at e M < 0.33 (Figs. 3A and 3E). From the comparisons in Fig 3A, it can be seen the observed patterns of creases, wrinkles, folds, period doubles and ridges follow consistently with the prediction of the phase diagram. Next, we discuss the observed patterns in film-substrate structures with moderate adhesion energies, which allow delamination between films and substrates. To compare the experimental observations with the calculated phase diagram, we section the three-dimensional phase diagram at various values of normalized adhesion energy, i.e., C/(m s H f ) 5 0 .13, 0.28, 0.46, 0.81, 3.99 and 66.63 (Fig. 4). The phase boundaries between delaminated-buckle and other patterns are highlighted as red curves in these sections. From Fig. 4, it can be seen that the observed transitions of phases with the increase of e M indeed follow the calculated phase diagram, for various values of C/ (m s H f ) and m f /m s . In particular, the delaminated-buckle can coexist with other instability patterns, and the phase boundaries between delaminated-buckle and other patterns can consistently predict whether delamination occurs in the film-substrate structures.
Finally, the wavelengths and amplitudes of the instability patterns are also validated by our experimental results. As shown in Fig. 5A, the experimentally observed wavelengths of the wrinkles, folds, period-doubles and ridges at varied mismatch strains match consistently with the predictions from Eqs. 9-12. In addition, combining Eq. 13 and l C wrinkle in Eq. 8, we compute the amplitude of the wrinkle for various m f /m s and e M in Fig. 5B. The calculated amplitudes of wrinkles match well with the experimental data for low modulus ratios, and the discrepancies for high modulus ratios are within 10% (Fig. 5B).

Discussion
The three-dimensional phase diagram is not only validated by our experimental results, but also by reported data on surface instabilities in various biological and biomimetic film-substrate structures. Figure 6 summaries the growth or swelling induced surface instabilities from a number of biological and biomimetic film-substrate systems, including animal tissue growth 14,16,18,21,22 , epithelial monolayer growth 19,25,26 , blood cell growth 27 , fruit growth or shrinkage 3,6,8 , biofilm growth 9-11 , and swelling of biomimetic hydrogels and elastomers 24,34,52-57 . The reported or estimated film-substrate modulus   Table 1 (see details in SI).
From Fig. 6, it can be seen the reported surface instability patterns indeed follow the three-dimensional phase diagram. If the film and the substrate are well bonded (Fig. 6A), swelling gels constrained by a rigid substrate on the bottom with m f /m s , 10 23 develop creases on the surfaces 34,[52][53][54]56 . The growing tissues 14,16,18,19,21,22,25,26 , tumors 24 and blood cells 27 , plant skins 3,6,8 , and mammalian skins 3 , with m f /m s in the range of 10 0 -10 4 , generally develop wrinkles in the systems with small mismatch strains, i.e., e M 5 0.0520.25; the wrinkles may further bifurcate into folds in the growing tissues 14,16,18,22 with larger mismatch strains (e M . 0.3), or transit into period-double in the mucosal guts 21 and on severely drying fruits 8 . The transition from wrinkles to folds has also been validated in swelling hydrogel bilayers 24,55 . In addition, in a system consisted of a swelling elastomer film on an underlying elastomer substrate (m f /m s < 10 4 ), the transition from wrinkle to ridge has been observed 57 .
If the adhesion energy between the film and substrate is relatively low (Fig. 6B), multiple delamination patterns have been observed in growing biological systems. For example, growing biofilms may delaminate to form buckle regions to facilitate nutrient transportation [9][10][11] ; the epithelial cells may delaminate due to overcrowd or self-metabolism 25,26 ; the swelling hydrogels loosely bonded on substrates may delaminate to form delaminated-buckles 54 . Since the mismatch strains in these phenomena are generally less than 0.4, we stack the three-dimensional phase diagram in the region of e M 5 0.1 2 0.4 into a two-dimensional phase diagram (Fig. 6B), where the boundary between delaminated-buckle and other phases is represented by a grey band. From Fig. 6B, it can be seen the reported delaminated-buckle patterns indeed fall in the delaminated-buckle region predicted by the phase diagram.
It should be noted that the phase diagram presented in this paper (Fig. 2) is valid for a wide range of dimensional parameters of growing biological systems. From Fig. 6, Table 1 and Supplementary methods, we can see that the film thickness varies from nanometer to centimeter, the modulus from Pascal to Gigapascal, the adhesion from 10 22 J/m 2 to 10 3 J/m 2 and the mismatch strain e M from 0 to ,0.6. The cases with excessively large film thickness (.m), large modulus ratio (m f /m s . 10 5 ) and large mismatch strain (.0.6) have not been considered in the phase diagram.
In summary, we present a three-dimensional phase diagram that can quantitatively predict various modes of growth-induced surface instabilities in biological film-substrate structures. By combining theory, computation and experiment, we show that the initially flat biological layers can systematically transform into instability patterns of wrinkle, crease, fold, period-double, ridge, delaminated-buckle and their coexistences, depending on three non-dimensional parameters: mismatch compressive strain, film-substrate modulus ratio, and normalized adhesion energy. The three-dimensional phase diagram offers a unified model for understanding morphogenesis in biological film-substrate structures on a mechanical base. The method for constructing the phase diagram opens new venues to study the formation of more complicated patterns, for example, in multi-layer structures, structures with intrinsic surface curvatures 6,58 , and inhomogeneous and/or anisotropic structures 6,8,16,18,59,60 . To the end, the phase diagram can potentially guide the rational design of a variety of biomimetic topographical-structures for engineering applications as diverse as disease therapy 22,24 , active cell culture 34 , biofouling management 35 , tunable superhydrophobicity 36  Experimental procedure. To simply induce a mismatch strain between the film and substrate, we first uniaxially pre-stretched the substrate to a prescribed ratio of L f /L s (Fig. 1Cii). The Sylgard film is then carefully attached on the pre-stretched substrate by uniformly pressing it on the substrate with two rigid plates, and the bilayer then is baked at an oven with a controlled temperature for 10 min. Thereafter, the prestretched substrate is gradually relaxed to the original length with a strain rate around 1 3 10 23 s 21 . The surface morphology is captured by a camera (Cannon, USA) with a tilted angel 20u. To enhance the film-substrate adhesion, a very thin layer of uncured Ecoflex can also be smeared on the substrate prior to attaching the films 48,51 .
Measurement of the adhesion energy. The adhesion energy between the film and substrate is measured with the peeling test ( Supplementary Fig. S12) 62 . A strip of Sylgard film with width b is carefully attached on a prestretched substrate, and the bilayer then is baked at an oven with a controlled temperature for 10 min. A force F is then applied to the Sylgard film to peel off the substrate along an angle h at a low peeling rate 1 3 10 24 m?s 21 . The adhesion energy is calculated by C 5 (1 2 cos h 1 e d /2)F/b, where e d is the strain in the detached section of the Sylgard film.
Finite element calculations. The finite element calculations (Supplementary Figs S3-S11) are implemented by software ABAQUS 10.1. The film and the substrate are taken to be incompressible neo-Hookean materials with shear moduli m f and m s respectively. The thickness of the substrate is much larger than that of the film (.20 times). Three methods are used to induce mismatch compressive strains in the films: In the first method, a film-substrate laminate is prestretched by a ratio of L f /L s , and then the strain in the film is released, followed by subsequently relaxing the laminate to the original length L s (Supplementary Figs S3, S4, S5A, S6, S7A and S8-S11). The second method involves swelling of the film in the horizontal direction ( Supplementary Figs S5B and S7B), while the third method is shrinking the substrate