Mechanics of fragmentation of crocodile skin and other thin films

Fragmentation of thin layers of materials is mediated by a network of cracks on its surface. It is commonly seen in dehydrated paintings or asphalt pavements and even in graphene or other two-dimensional materials, but is also observed in the characteristic polygonal pattern on a crocodile's head. Here, we build a simple mechanical model of a thin film and investigate the generation and development of fragmentation patterns as the material is exposed to various modes of deformation. We find that the characteristic size of fragmentation, defined by the mean diameter of polygons, is strictly governed by mechanical properties of the film material. Our result demonstrates that skin fragmentation on the head of crocodiles is dominated by that it features a small ratio between the fracture energy and Young's modulus, and the patterns agree well with experimental observations. Understanding this mechanics-driven process could be applied to improve the lifetime and reliability of thin film coatings by mimicking crocodile skin.

Fragmentation of thin layers of materials is mediated by a network of cracks on its surface. It is commonly seen in dehydrated paintings or asphalt pavements and even in graphene or other two-dimensional materials, but is also observed in the characteristic polygonal pattern on a crocodile's head. Here, we build a simple mechanical model of a thin film and investigate the generation and development of fragmentation patterns as the material is exposed to various modes of deformation. We find that the characteristic size of fragmentation, defined by the mean diameter of polygons, is strictly governed by mechanical properties of the film material. Our result demonstrates that skin fragmentation on the head of crocodiles is dominated by that it features a small ratio between the fracture energy and Young's modulus, and the patterns agree well with experimental observations. Understanding this mechanics-driven process could be applied to improve the lifetime and reliability of thin film coatings by mimicking crocodile skin. S kin fragmentation (also known as crocodile cracking), is a common way for thin films to release deformation energy that is caused by tension [1][2][3] . This process is irreversible and leaves complex pattern as many cracks interact and go through the material surface. The geometry and size of those fragmentation patterns varies dramatically for different materials including organic and inorganic systems such as fruit skin, paintings and ceramics 4,5 , as well as a nanofilm with atomic thickness 6 (Fig. 1). They are generally considered a threat to engineering systems. For instance, steel and wood structures in buildings and turbines with fragmented coating materials are largely exposed to environment and subjected to corrosion. It is critical to understand the mechanism behind such phenomena crossing multiple scales. Especially, how tensile forces initiate the fragmentation in thin films and what material characteristic stabilizes the pattern remain largely illusive. Fundamental answering of those questions can help the design of surface coating with improved lifetime and resistance.
The fragmented skin on a crocodile's head provides a suitable case to investigate the skin fragmentation. The skin acts like scales of snakes, lizards and fishes 7,8 to protect the crocodile, but instead of gene regulation, its pattern is induced by skeletal growth during the embryonic stage 9 . This contradicts with the conventional thinking, suggesting that some fragmentation may actually enhance the protection function. Moreover, it is intriguing, but unclear, what is the mechanical basis for crocodiles to generate this skin fragmentation with a characteristic size. By learning from the fragmentation of crocodile skin, we expect to gain knowledge in improving the design of engineered surface coatings.
In our current study, we combine theoretical analysis with numerical simulation that uses a simple thin film model for the crocodile skin. We note that prior studies have investigated the fragmentation in surface coatings using numerical methods, but most of them require pre-existing defects for crack initiation [10][11][12][13] . It is difficult to decouple the effect of random distributions of weak regions and the effect of material property on the characteristic of fragmentation pattern by using those models. Here, we use a simple elastic network model with a uniform spring stiffness and strength (detail in the Methods section), without additional defects, to investigate how the characteristic size of fragmentation is affected by mechanical properties of the film material and loading conditions. It is shown that by adjusting the spring characteristic we can reproduce the characteristic of fragmentation pattern as is found in different materials including crocodile skin.

Results
Here we use a simple network model to study the fragmentation behavior of the thin film as depicted schematically in Fig. 1g (details see Methods section). We apply biaxial tensile strain to this model to investigate how fragmentation generates and grows. The fragmentation pattern is identified as an assembly of polygons as shown in Fig. 2a. Notably, cracking only appears beyond a critical level of deformation (denoted by the strain e f , where strain is defined as the length increase divided by the initial length), as shown in Figs. 2b. The deformation in the material decreases for strain beyond e f by creating more fragments and decreasing their sizes. We measure the average fragment size as a function of strain as shown in Fig. 2c (detail in Methods section). It is shown that for strain beyond e f the fragmentation size keeps decreasing, until it reaches a secondary critical deformation state (denoted by the strain e c ) when the pattern becomes stabilized and strain increases exclusively, leading to further opening the already-cracked edges of polygons. This results in an asymptotic fragment size. Using fracture mechanics theory, we derive how such fragmentation follows deformation by relating the deformation energy in the material to the surface energy used to create cracks (details see Methods section). We find that the fragment size during stretching can be expressed by: where E is the material's Young's modulus, n is Poisson's ratio, 2c is the fracture energy per unit area of the material, e is the applied tensile strain caused by the mismatch between thin film and substrate, L ' is the length of the residue scale, which is the asymptotic size and a is the stiffening factor that relates to the nonlinear material property of skin (detail in the Methods section). From this equation, we find that the evolution of the fragment size is governed by two factors: the applied strain as well as the ratio between the fracture energy and the material stiffness. The comparison of the importance of those two factors is adjusted by the nonlinearity given by a. We summarize some plausible values for 2c/E for different skins (human hand, chicken and crocodile) measured in experiments as shown in Table 1. Those values are used to fit the simulation result as shown in Fig. 2c. We find that only small 2c/E value, as is observed in crocodile skin, gives a good interpretation for the simulation result. From this analysis we also obtain the empirical value of a as 2.4, which corresponds to a stiffening hyper-elastic material, as common in nature. Moreover, we find the best fit with unphysical negative L ' values for human hand and chicken, while physical positive L ' value for a crocodile. This result suggests that it is the mechanical force and releasing of deformation energy that dominates the fragmentation of crocodile skin and makes the pattern eventually reach the stable characteristic length of L ' (Fig. 1a). The negative L ' value, in contrast, suggests the surface geometry of those skins is not caused by fragmentation, as they can form corrugated pattern with more complex geometry but not asymptotic size (until the cellular size) (as shown in Fig. 1c). We change the mechanical characteristic of the network structure (changing from a 5 2.4 to a 5 1 and 7) by using different power exponent of the constituent springs and compare the characteristic nature of the fragmentation in those different systems, as shown in Fig. 3. It is seen that the fragmentation starts to appear in the linear elastic material (a 5 1) at a smaller strain (e f 5 0.05) than the hyperelastic material and the hyper-elastic material with strong strainstiffening characteristic (a 5 7) generate the fragmentation at a much larger strain (Fig. 3b). Moreover, for a 5 7, it is observed that the fragmentation starts with many disconnected small cracks as shown in Fig. 3a, suggesting that this strain-stiffening material is able to diffuse the deformation energy in the entire deformed region instead of concentrating at several crack tips 14 . Further increasing the deformation for this material quickly lead to the fragmentation pattern with many small fragments (L 5 0.2 cm as shown in Fig. 3b). We measure each fragment length of the three materials (a 5 1, 2.4 and 7) at the end of simulation of e 5 0.2 and obtain the probability distribution as shown in Fig. 3c. It is found that both a 5 1 and 7 lead to smaller fragments with more uniform (L max /L 5 1.9 and 1.8, for a 5 1 and 7, respectively, where L max is the maximum length and L is the average length of all the fragments) and symmetric distribution of the fragment length than a 5 2.4 which contains fragments much larger than the rest (L max /L 5 2.5). The result obtained for a 5 1 quantitatively agrees with the experimental observation for oxide coating, which is taken as a linear-elastic material, as L max /L 5 1.65 , 1.8 is obtained for this material under different applied strains 12 . These results agree with what we observe for crocodile skin and other thin films ( Fig. 1) as the fragments in crocodile skin ( Fig. 1a) are larger and have more irregular length than in linearelastic material (ceramics and surface coating in Figs. 1f and b, respectively) or hyper-elastic material with strong strain-stiffening  characteristic (human skin in Fig. 1c), suggesting again that the fragmentation pattern of crocodile skin is strictly governed by its mechanical properties.
We now change the strain rate and investigate how the deformation speed alters the fragmentation pattern. The stabilized patterns under multiple strain rates are as shown in Fig. 4a, where it is demonstrated that a faster strain rate induces more fragments with smaller size. The asymptotic size of polygons L ' , as systematically shown in Fig. 4b, more clearly illustrates that the fragmentation is strongly rate-dependent at small deformation rate while becomes insensitive of large deformation rate. By considering-as posed by classical viscoelasticity-the scaling E(t) , E(t 5 0)/(t/t 1 1) (t is a characteristic time) 15 , since e~ð _ e t ð Þdt~_ et where t is the time for applying deformation of constant strain rate, assuming the asymptotic size is proportional to the fragment size at largest strain we obtain that L ? !c Ee az1 e c is the applied strain for obtaining asymptotic fragment size). This relationship is used to fit the simulation results as shown in Fig. 3b. Indeed, this mechanism, combined with the anisotropic growth of the head skull, explains the anisotropic and irregular fragmentation pattern seen in crocodiles (Figs. 5h and i). This result also suggests that the fragmentation pattern varies for paints because of the different deforming history caused by environmental change, generating the unique fingerprint on each of them. It is also noted that for small loading rate, the fragment size approximates exponential scaling, but the size departs from the exponential at larger loading rate because the resolution of the network structure for generating the fragments is limited by the initial length of the springs in the model, which has the physical meaning of the size of the constituent particles such as cellular size of the skin. This factor can be involved by combining experimental measurements of cellular size to improve the biological insight of the model in future.
We systematically compare the simulation result against in situ observations of the fragmentation of crocodile head, as illustrated in Fig. 5. It is clearly seen that the emerging fragmentation pattern very well reproduce what is observed on the crocodile head. The connectivity of the polygons has an average value of 3.03 to 3.21 (with and without boundary), and this connectivity range yields a mean number of sides of polygons between 5.88 and 5.31. This is in excellent agreement with the experimental observation of 5.532 ref. [9], and indicates that polygons are nearly hexagonal and those closer to the boundary have more sides. Moreover, the model shows that introducing heterogeneous deformation-considering that the strain rate (or bone growth) in one direction is faster than that in the other, orthogonal direction-captures the characteristic laddering shape on the top of crocodile's head as illustrated in Figs. 5h and i. This finding demonstrates that the laddering scale pattern observed on crocodile head is caused by the fact that the head's longitudinal growth is faster than in other directions 9,16 . To demonstrate the significance of the  fragmentation of crocodile, we summarize the connectivity of the fragmentation pattern of different systems as shown in Fig. 6. It is shown that the fragmentation of crocodile head has similar geometry as is observed in cantaloupe, ceramics and paint and their patterns are all captured by our mechanics model. The significant difference for the skin crumpling of hand skin indicates that the corrugated pattern with complex geometry is governed by involving other mechanisms rather than tension-induced fragmentation.

Discussion
In this study, we identify two critical strains (e f for initiating of fragmentation and e c for reaching the asymptotic fragment size), which are induced by skeletal growth, that govern the initial and final mechanics of skin fragmentation. Those two critical strains divide the deformation procedure into three regimes, which agrees with what is observed in the cracking of other thin film systems 17 , except that our results do not statistically depend on the initial random distribution of defects. The first critical strain explains why there is no fragments found before 45 days of embryonic crocodile, as is found by in situ observation 9 , likely because the accumulated strain remains below the threshold (e f ) during that period. The second critical strain explains why the pattern does not change (or redevelop anew) repeatedly over time, as no new polygons are created by the increasing strain above the other threshold (e c ). Our result also shows that the fragmentation depends on the ratio between the fracture energy and Young's modulus of the skin. The mechanical property of crocodile skin makes it unique to form the polygon pattern during development. Flexible skins of other animals, such as mammal and bird, are hard to reach a stable fragmentation pattern in the skin under normal state. However, there are extreme conditions 18 , for   example dehydration, oldness and illness [19][20][21][22] , which can alter the mechanics of the skin and make the stable fragmentation possible. From a different angle of view, the skin fragmentation of crocodile can be attributed to mechanical adaption 23 . The hard skin of crocodile with fragmentation can greatly increase the flexibility (with reduced bending stiffness because cracks set fragments free-by reduced entropic elasticity 24 -to move out of plane) and reduce the difficulty for movements at the end of jaw. A similar strategy is also observed to be adopted by armoured fish as its discrete scales for body protection and mobility 20,25 . The fragmentation on crocodile heads is unique for animal skins but shows a similar geometry as many inorganic thin films. Here, we investigated the mechanism hidden behind these phenomena by introducing a mechanical model and probing the entire fragmentation process. Our results demonstrate that the ratio between the fracture energy and Young's modulus together with deformation rate governs the characteristic size of the fragmentation pattern. By using this knowledge learnt from crocodile skin, we may design and produce surface coating with improved stability by using synthetic materials. For instance, we may use heterogeneous distribution of stiff and soft materials, which mimic the crocodile scales and their joints, for coating under the help of digital manufacturing technologies. By using the strategies one can make the coating stable from further fragmentation and hence have longer lifetime and provide more efficient protection.

Methods
Details of the elastic network model. The model we use to simulate the fragmentation behavior of keratinized epidermal under biaxial tension is a simple elastic network model. The model is composed of a collection of beads and inter-bead connections ref. [14,26]. The initial coordinates of all the beads are randomly distributed. The topology of the inter-bead connection is designed by Delaunay triangulation algorithm 27 . This method enables generating triangle meshwork from randomly distributed point and maximizes the minimum angle of all the angles of the triangles, avoiding irregular triangles. For each inter-bead connection in the triangulation, we use a breakable nonlinear springs to represent the interaction. The tensile force for each spring is given by f 5 dk 0 (Dr/r 0 ) a where r 0 is the initial length of the spring, Dr is the deformation of the spring, k 0 is the stiffness of a unit length spring, a is the stiffening factor that relates to the nonlinear property of the spring, d is a cut off function with its value given by d~exp Dr with J~300 as the smooth factor and e b 5 0.05 as the bond breaking strain. It is noted that by using this force function, the stiffness of each spring stochastically varies with its initial length, ensuring the uniform mechanics of the material. Without loss of generality, all parameters in the model, including bead mass, spring stiffness and film length are set unitless (without any fitting). We apply biaxial tensile strain to deform the entire film with an overall constant strain rates. Yet the deformation is applied in a discrete from: for every a hundred steps, we add the constant strains to the film at the first step and equilibrate the film during the rest steps with the single layer of beads at the four edges fixed. This loading method mimics the tension field caused by skeletal growth and applied on the epidermal on the crocodile head as illustrated in Fig. 1g (considering that the bone is 1-2 orders of magnitude stiffer than the epidermal and it can be modeled as a rigid substrate to generate homogenous deformation to the epidermal in every single direction 28 ). This model enables us to systematically investigate how growth history and growth rate are coupled 29 to determine the evolution of fragmentation pattern.  Monitoring fragmentation and computing its size in the simulation. During the post-analysis of the simulation we monitor the fragmentation pattern as assemblies of polygons by identifying all rupture springs (Dr/r 0 $ 0.05) and highlight them by taking snapshots. To compute the size of those polygons, we simply assume the cracks are in form of hexagonal lattices and the total lattice area is given by A tot~3 ffiffi ffi 3 p L 2 n 8 where n is the number of scales and L/2 is length of the hexagonal edge (L is the scale size), while the area of the polygon edge with width r 0 (the average spring length at equilibrium) is A edge~3 Lr 0 n=2. Noting that the ratio between the area of the polygon edge and total surface area equals the ratio between the number of beads involved in broken bonds (N crack as a function of strain) and the total bead number (N). These relations enable us to obtain the fragment size according to We have measured the distribution of fragment size for cases a 5 1, 2.4 and 7 and summarized the result in Fig. 3c. It is shown that the average value and standard deviation of the fragment is 1. Fracture mechanics analysis for the fragmentation size. Here we derive how fragmentation size follows deformation. The energy balance during crack propagation imposes that the variation of the total potential energy dP, equal to the variation of the elastic strain energy dW minus the variation of the external work dY, must be equal to the opposite of the work spent in creating the new surface of fragments dS, i.e dP 5 dW 2 dY 5 22cdS, where 2c is the fracture energy per unit area. Assuming linearity for the constitutive law of the film implies the validity of the Clapeyron's theorem 30 , i.e. dW 5 dY/2. Accordingly, dP 5 2dW and the condition for fragmentation becomes dW 5 2cdS. Note that this result is in general valid also for nonlinear systems under imposed displacements; for this case, in fact, the external work is identically zero and thus we do not need to apply the Clapeyron's theorem to obtain it. By integration, we have DW 5 2cDS, where DS is the new crack surface area created as DS 5 3/2nLb (for hexagonal fragmentation geometry) where n is the number of scales, b is their thickness and L is their size. We now generalize the stressstrain relation as s 5 Ee a by considering the nonlinearity of materials and estimate the deformation energy as DW 5 2AbEe (a11) /(12n)/(a 1 1), where e is the accumulated mismatch biaxial strain between skin and skeleton, A is the film surface area, b is its thickness, E is the secant modulus that equals to Young's modulus of linear material or general materials at small deformation, n is Poisson's ratio and a is the stiffening factor that relates to the nonlinear material property of skin (a , 1 for elastic-plastic, a 5 1 for linear, a . 1 for hyper-elastic material). Accordingly, the fragmentation is expected when the following system holds: where e f is the skin fracture strain as shown in computational modelling and A ' is the residue area that still subjects to deformation after fragmentation that we assume to be proportional to the fragment perimeter where the residual deformations concentrate 31 . Accordingly, since for the hexagonal geometry L 2~8 A 3 ffiffi ffi 3 p n , we can define LL ?~8 A ? 3 ffiffi ffi 3 p n and the fragment size during stretching is given by Eq. (1).