## Abstract

We present mechanics of surface creasing caused by lateral compression of a nonlinear neo-Hookean solid surface, with its elastic stiffness decaying exponentially with depth. Nonlinear bifurcation stability analysis reveals that neo-Hookean solid surfaces can develop instantaneous surface creasing under compressive strains greater than 0.272 but less than 0.456. It is found that instantaneous creasing is set off when the compressive strain is large enough, and the longest-admissible perturbation wavelength relative to the decay length of the elastic modulus is shorter than a critical value. A compressive strain smaller than 0.272 can only trigger bifurcation of a stable wrinkle that can prompt a setback crease upon further compression. The minimum compressive strain required to develop setback creasing is found to be 0.174. If the relative longest-admissible perturbation wavelength is long enough, then the wrinkle can fold before it creases, and the specimen can be compressed further beyond the Biot critical strain limit of 0.456. Various bifurcation branches on a plane of normalized longest-admissible wavelength versus compressive strain delineate different phases of corrugated surface configurations to form a *ruga* phase diagram. The phase diagram will be useful for understating surface crease, as well as for controlling *ruga* structures for various applications, such as designing stretchable electronics.

## 1. Introduction

Free-surface deformation instability of a soft solid, with its modulus varying in depth, often develops surface patterns, such as wrinkles, creases, ridges or folds, under in-plane compression. Here, a large-amplitude state of a wrinkle, crease, ridge or fold is collectively denoted by the Latin word *ruga*, and the study of motions associated with configurational changes among such states is called *ruga mechanics*. The motions include progressive amplitude growth of a *ruga* and the snap jump of the amplitude from one state to another. For example, it is commonly observed that a free surface of a homogeneous soft solid develops a crease with a snap-buckling instability, and a soft substrate with a stiff layer progressively grows stable wrinkles [1,2]. So far, much of *ruga* mechanics research has been limited to post-bifurcation analysis of relatively small-amplitude wrinkles as incipient of *ruga* structures. In this paper, we present post-bifurcation stability analysis of various *ruga* states up to very large amplitudes for a neo-Hookean solid [3,4], with its modulus decaying exponentially from the free surface [5,6].

Regarding small-amplitude stability analysis of a flat surface, Biot [7] showed, in his linear bifurcation analysis, that the flat surface of a homogeneous neo-Hookean solid becomes unstable at the Biot critical strain *ε*_{B} of 0.456 plane-strain compression. However, this critical strain can be hardly reached experimentally because the flat surface unstably creases at a strain lower than the Biot critical strain. Gent & Cho [1] observed in their experiments that a rubber-strip surface compressed by bending creased at a nominal strain in the range 0.28–0.42. Later, Hohlfeld [8] used finite-element analysis (FEA) to show that such a crease is typically formed through a snap-buckling process at the Biot critical strain, and that post-bifurcation unloading shrinks the crease depth until it disappears at a compressive strain of 0.35. Subsequently, Hong *et al.* [9] could trigger progressive crease growth by pre-pinching the surface in their FEA so as to follow the stable subcritical creasing branch of pitchfork bifurcation beyond the compressive strain of 0.35. Cao & Hutchinson [10] studied the imperfection sensitivity of creasing by analysing multiple-mode post-bifurcation stability, with a first-order displacement perturbation in the Koiter approach [11,12]. They found that the post-bifurcation wrinkle of an arbitrary wavelength is unstable, and the critical strain is sensitive to the imperfection amplitude. They also noted in their FEA that the local strain to set off creasing is always the Biot critical strain.

Now, questions arise. What distribution of the substrate modulus can stabilize the post-bifurcation wrinkle? If it can be stabilized, can we sustain the stability up to an arbitrary level of compression? If the stable post-bifurcation wrinkle becomes unstable at a large compression, does the second post-bifurcation *ruga* become stable or unstable? If unstable, can the unstable bifurcation develop period-multiplication modes [13]? To answer these questions systematically, we use a neo-Hookean solid model with an exponentially decaying modulus (EDM) for the study of *ruga* development. Then, analysis of the EDM model provides the deformation characteristics of a half-space from those of a homogeneous medium to those with a stiff film [14–16] as a single-parameter family. The model was also used for different purposes by Biot [5] and Lee *et al.* [17]. Behaviour of the EDM model is analytically tractable, and the results are transformative for making interpretations of homogeneous material behaviour by taking the limit of the decay length to infinity. It is shown in this paper that the EDM model can capture not only the behaviour of *instantaneous creasing*, but also that of *setback creasing* delayed by progressive growth of wrinkles or folds. It is revealed in later sections that Koiter's approach with second-order perturbation analysis can determine whether the nonlinear bifurcation mode is a progressive wrinkle or an instantaneous crease. In turn, the Koiter analysis provides a criterion to delineate instantaneous creasing from setback creasing. However, once the wrinkle amplitude becomes large, nonlinear deformation processes in the post-bifurcation growth of wrinkles or folds are hardly tractable analytically. Therefore, FEA is used to analyse setback creasing.

Figure 1 shows the schematics of possible creasing processes on an EDM solid surface under isoperiodic compression. The creasing is found to be a snap-buckling process. The creasing process is, in general, sensitive to loading boundary conditions and/or imperfections. We thereby limit our analysis of bifurcation processes to those caused by uniaxial isoperiodic compression. In isoperiodic compression of EDM solids, the longest-admissible wavelength is typically set by the domain size in a numerical analysis or the specimen size in an experiment. In addition, the periodic boundary condition suppresses period multiplying bifurcations, because the wrinkle with the longest wavelength has the lowest energy at a given compressive strain. As a consequence of analysing isoperiodic compression of an EDM solid, we can construct a *primary ruga phase diagram* on a plane of the smallest admissible nominal wavenumber and the nominal compressive strain. In the phase diagram and as depicted in figure 1, the phase boundaries are constructed by analysing transitions from one phase to another, following isoperiodic transitions from a flat to a wrinkle (*C*_{12}), a flat to a crease (*C*_{13}), a wrinkle to a crease (*C*_{2c}), a wrinkle to a fold (*C*_{24}) and a fold to a crease (*C*_{4c}).

## 2. Perturbation analysis of creasing processes: instantaneous versus setback creasing

Using a linear perturbation analysis, Biot [6] examined bifurcation of a flat neo-Hookean EDM solid surface to a sinusoidal wrinkle under plane-strain uniaxial compression. He determined the bifurcation wavenumber as a function of the applied compressive strain. However, the linear perturbation analysis cannot determine the stability of post-bifurcation modes and their amplitudes. Here, we revisit the Biot problem in §2*a* and extend the study to nonlinear bifurcation analysis of the wrinkle by carrying out Koiter analysis with second-order perturbation solutions in §2*b*.

### (a) Linear perturbation analysis of surface wrinkles on a neo-Hookean exponentially decaying modulus solid

For the perturbation analysis, the material coordinate (*X*_{1},*X*_{2}) and the distribution of the shear modulus
2.1ain the substrate are depicted in figure 1*a*, where *Q*_{0} is the shear modulus at the surface, and the coefficient *a* represents the reciprocal of the decay length of the modulus. In the deformation analysis, the plane-strain mapping from the stress-free reference configuration *X*_{i} to a deformed configuration *x*_{i} is denoted by
2.1bHere, *λ*_{(i)} are the in-plane principal stretch ratios of the fundamental state aligned in (*X*_{1},*X*_{2}) directions, and *u*_{i} the plane-strain perturbation displacement components. For a neo-Hookean material, the strain energy of a plane-strain deformation is defined as
2.1cfor variations of *u*_{i}(*X*_{1},*X*_{2}) and *p*(*X*_{1},*X*_{2}). Here, *F*_{ij} is the two-dimensional deformation gradient tensor defined as *F*_{ij}=∂*x*_{i}/∂*X*_{j}, *J*=*det*[*F*_{ij}] the Jacobean, and *p* the hydrostatic pressure that plays the role of a Lagrangian multiplier function to enforce incompressibility in variational analysis. Summation convention is applied to repeated indices throughout this paper. The principal stretch ratios of the fundamental state are set as *λ*_{(1)}=*λ* and *λ*_{(2)}=1/*λ* for incompressibility. For the perturbation analysis, the displacement *u*_{i}(*X*_{1},*X*_{2}) and the hydrostatic pressure *p*(*X*_{1},*X*_{2}) are expanded with respect to a perturbation parameter *ζ* up to the second order,
2.2aand
2.2bwhere and *p*^{(0)}(*X*_{1},*X*_{2}) are, respectively, the uniform displacement from the fundamental-state configuration and the hydrostatic pressure in the fundamental state. The parameter *ζ* is taken to be the amplitude of the first-order vertical displacement component at the top surface, . The superscripts in parentheses indicate the perturbation order. The first Piola–Kirchhoff stress, defined for an incompressible solid as
2.3is used with the independent pressure expression in (2.2b) to obtain the explicit form of the stress field up to the second order as given in (A 12).

Then, the equilibrium equations up to the second order are reduced from (A 12) to
2.4aand
2.4band the incompressibility condition up to the second order to
2.4cThe first-order displacement components and hydrostatic pressure in the perturbed configuration are determined by solving (2.4), as presented in appendix A,
2.5a
2.5b
and
2.5cwhere
2.6aand
2.6bwhere is the normalized wavenumber in the undeformed configuration. The exponents *α*_{3} and *α*_{4} have negative real parts. Thus, the constants *C*_{3} and *C*_{4} in (2.5) are set to zero to ensure that the perturbed displacement components vanish at .

Enforcing the traction-free boundary conditions at the top surface of the substrate up to the first order, we construct the following highly nonlinear eigenvalue problem:
2.7Then, the characteristic equation provides the normalized bifurcation wavenumber as a function of the stretch ratio for *λ*<1 as
2.8

Figure 2*a* shows the normalized bifurcation wavenumber in the undeformed configuration in terms of the applied nominal compressive strain, *ε*=1−*λ*. The solid curve in figure 2*a* indicates the formation of stable wrinkles for compressive strains smaller than 0.272, the transition strain at *M*. The dashed curve represents the formation of unstable wrinkles under compressive strains between *ε*_{M} (=0.272) and the Biot critical strain, *ε*_{B} (=0.456). The wavenumber diverges to infinity at the Biot critical strain, indicating that the homogeneous neo-Hookean half-space unstably bifurcates at the Biot critical strain for perturbations with an arbitrary wavelength. The chain curve shows a lower bound of the crease strain, which could be gradually reached when the bifurcated crease was unloaded in our FEA.

### (b) Second-order perturbation analysis for nonlinear stability of a wrinkle

In this section, nonlinear stability of the wrinkles along the solid and the dashed curves presented in figure 2*a* is examined by second-order perturbation analysis. A solution procedure of the second-order perturbation is described in appendix A, and only the final expressions of the second-order displacement components and the hydrostatic pressure are presented here for further discussions on energy variations,
2.9a
2.9b
and
2.9cwhere explicit expressions of *F*_{m} (*m*=1,2,3) are given in appendix A(c), and *B*_{m},*H*_{1}(*X*_{2}) and *H*_{2}(*X*_{2}) are given in appendix B.

Nonlinear stability of the wrinkle near the bifurcation point is then tested by the energy variation with respect to the displacement variations up to the second order. The energy variation, with *λ* and fixed at a bifurcation point on the bifurcation curves in figure 2*a*, can be obtained from (2.1), (2.2), (2.5) and (2.9) with respect to the displacement variation as
2.10Here, *φ* denotes the strain energy per period, which is the same as the potential energy per period for the variation, and the superscript ‘*b*’ indicates that the quantity is evaluated at the bifurcation point. In (2.10) as well as in the following (2.11) and (2.12), the chain rule operation implies variation *δφ*^{b} of the functional *φ*^{b}, caused by variation of the displacement field, . Because the first- and the second-order displacement fields are orthogonal for the energy variation and the other terms up to the *ζ*^{3} term in (2.10) vanish for equilibrium at a bifurcation point, the expression of the nonlinear bifurcation energy variation is reduced to
2.11The terms on the right-hand side of equation (2.11) are determined for a neo-Hookean EDM solid as
2.12aand
2.12bwhere *i*, *j*, *k*, *l*, *m*, *n*=1 or 2, and a linear operator means (*U*^{(m)}*W*^{(n)})=*U*^{(1)}*W*^{(2)}+*U*^{(2)}*W*^{(1)}.

Equation (2.12a) exhibits the explicit coupling effect of and for the energy variation, whereas (2.12b) entails the effects of implicit coupling between the first- and second-order fields.

Figure 2*b* shows the nonlinear bifurcation energy variation along the solid and the dashed bifurcation curves in figure 2*a*. As shown in figure 2*b*, the energy variation is positive at every bifurcation strain below the nominal strain *ε*_{M}=0.272. This indicates that the bifurcated wrinkle is stable near the bifurcation point within 0<*ε*<*ε*_{M}, for which . However, the energy variation becomes negative for the bifurcation points within *ε*_{M}<*ε*<*ε*_{B} and . The result implies that the bifurcated wrinkle is unstable near its bifurcation point in this range, and its growth includes simultaneous unstable growth of the second-order mode of the wrinkle. Once the wrinkle becomes simultaneously unstable with the second-order mode, the wrinkle becomes unstable, regardless of all higher-order modes, because all the modes higher than the second order will contribute to the energy variation only up to *O*(*ζ*^{5}). As we can observe in (2.5) and (2.9), the second-order mode escalates the local curvature and compressive strain at the coordinate origin, localizing the valley of the wrinkle. It is observed in the following sections that the nonlinear bifurcation instability up to the second-order mode implies creasing of the surface. Consequently, the M-point (*ε*_{M}=0.272; ) breaks up Biot's bifurcation curve for 0<*ε*<*ε*_{B} in figure 2*a* into a progressive-wrinkle curve for 0<*ε*<*ε*_{M} and an instantaneous-crease curve for *ε*_{M}<*ε*<*ε*_{B}.

## 3. Finite-element analysis of post-bifurcation wrinkle growth and setback creasing

In §2*b*, second-order perturbation analysis predicted that isoperiodic compression of a flat EDM solid surface sets off progressive wrinkling or instantaneous creasing. In this section, we carry out finite-element method (FEM) simulations to trace the bifurcation modes and corresponding post-bifurcation processes. Isoperiodic compression is simulated for various wavenumbers in the range .

The FEM simulations are conducted as a two-dimensional plane-strain problem in ABAQUS/standard using the stabilized Newton–Raphson method. Hybrid elements CPE4RH are adopted to simulate the incompressible neo-Hookean material. For the isoperiodic compression process, periodic boundary conditions are enforced on the two sides of a planar unit cell in the (*X*_{1}, *X*_{2})-plane, with width *L* and depth *H*=5*L*, as shown in figure 3*a*. To introduce surface imperfections, the simulations are first run by enforcing a sinusoidal displacement boundary condition at the top surface as
3.1where *L* is chosen as the imperfection wavelength in the undeformed configuration, and as the relative imperfection amplitude. Then, the deformed geometry is imported, and the imperfection command in ABAQUS is used to create a perturbed mesh that conforms to the surface undulation of the half-space. A frictionless self-contact interaction is applied at the top surface of the half-space to capture the crease and fold behaviour without self-penetration.

The classical Lagrange multiplier method of constraint enforcement was adopted using the hard contact option in ABAQUS. The minimal mesh size is about 1/2000 of *L*. The shear modulus distribution (*Q*=*Q*_{0}e^{aX2}) is built up in the ABAQUS program by assigning a fictitious temperature distribution *θ*=*θ*_{0}e^{aX2} with a linear relationship between *Q* and *θ*, and a vanishing thermal expansion coefficient.

As shown in the left inset of figure 3*a*, the isoperiodic compression process is simulated by compressing the domain with periodic boundary conditions and small incremental compression steps, *Δε*, less than 0.0005 or 0.0008. In the simulation, the normalized wavenumber is controlled by varying *a* for a fixed domain length *L* (=2*π*/*k*). Figure 3 shows five different *ruga* configurations, observed near various bifurcation points during isoperiodic compression of an EDM solid with three different , 0.471 and 0.100). The colour plot exhibits the level of the nominal (or Biot) strain . In each frame, a single period result is repeatedly displayed twice to show a configuration of two periods. Figure 3*b* shows the post-bifurcation configuration of an instantaneous crease near the bifurcation point for , whereas figure 3*c* presents the wrinkle near the Biot bifurcation point and figure 3*d* the setback crease at the wrinkle valley for . This is the normalized wavenumber with which a setback crease can be generated by the minimum applied strain, *ε*_{T}=0.174. Below this strain, a crease cannot be formed, regardless of whether it is an instantaneous or a setback crease. This critical strain contrasts *ε*_{M}=0.272, below which an instantaneous crease cannot be created. Figure 3*e* and *f* respectively display a folding configuration and setback creasing at the bottom of the fold channel for .

Because we use an imperfection on the surface profile as shown in (3.1), we inspect the evolution of the surface profile by taking the Fourier-series expansion of the deformed surface profile as a function of *X*_{1} at every step of the compression, as
3.2for *N* up to 8, where *y*_{0} is the profile-base constant. Figure 4*a*–*c* shows evolutions of the Fourier coefficients (*A*_{1}, *A*_{2}, *A*_{4}, *A*_{8}) as the domain is compressed with three different (, 1.30 and 2.51) held fixed, respectively. Figure 4*d* exhibits the growth rates of the lowest-mode coefficient, d*A*_{1}/d*εs*, as a function of *ε* for the three different compression processes. Neglecting the coefficients less than 10^{−5}, figure 4*a* shows that all the coefficients grow quickly near Biot's wrinkle bifurcation point (*ε*^{(i)}_{w}=0.157, ); however, the growth rates are finite, and all are stable. The growth of the coefficients indicates sharpening of the wrinkle valley. The growth rate distinctly defines the bifurcation strain *ε*^{(i)}_{w} and thus the Biot wrinkle bifurcation point, as shown in figure 4*d*. The bifurcation strain *ε*^{(i)}_{w} is relatively insensitive to imperfections, and coincides closely with the prediction of the linear perturbation analysis. However, growth rates of all the coefficients diverge at the crease point, regardless of an instantaneous or setback crease, and the instantaneous crease point is imperfection sensitive. On the other hand, a setback crease requires substantial amplitude growth of the wrinkle before it locally sets off unstable creasing, and thus the setback crease is relatively insensitive to initial long-wavelength imperfections. However, the setback crease is sensitive to imperfections with wavelengths shorter than the final local zone size of the crease in the wrinkle valley.

## 4. Discussion

In figure 5, three bifurcation curves, *C*_{12}, *C*_{13} and *C*_{2c} define four basic domains of surface configurations on the -plane. The curves *C*_{12}, *C*_{13} and *C*_{2c} correspond to the progressive-wrinkle, instantaneous-crease and setback-crease bifurcation curves, respectively. The four domains on the -plane constitute the basic structure of the primary *ruga* phase diagram, on which implies the homogeneous half-space limit and the thin-film limit. The phase diagram displays possible equilibrium configuration phases of an EDM solid surface under isoperiodic compression. Phase (I) denotes flat, phase (II) wrinkle, phase (IIc) setback-crease and phase (III) instantaneous-crease configurations. Regarding the setback-crease curve, as the compressive strain isoperiodically increases beyond Biot's bifurcation point, nonlinear instability is pre-empted, and a localized setback crease develops at the valley of a wrinkle, preventing the wrinkle from growing further. Collection of the new setback-crease bifurcation points produces the setback-crease curve *C*_{2c} on the phase diagram, the dashed line in figure 5, which delineates phase (IIc) from phase (II).

Furthermore, we found another phase of fold configurations, phase (IV). When we tracked the isoperiodic evolution of wrinkles with low wavenumbers, i.e. less than 0.471, we found that there is a critical wavenumber below which the wrinkles fold at *C*_{24} and form channels, as shown in figure 3*e*. Once the fold configuration is further compressed, a setback crease develops at the root of the fold channel (figure 3*f*). This creasing at *C*_{4c} separates the phase of a fold-crease configuration (IVc), from that of a fold configuration (IV), as shown in the inset of figure 5. Beyond this point, a question still remains whether this crease would further propagate to collapse the channel, defining another phase on the phase diagram. It is worth noting that the phases (IV) and (IVc) sit past Biot's singular strain limit. An EDM solid can be isoperiodically compressed beyond the Biot singular strain limit with a low wavenumber because growth of a large wrinkle amplitude can keep the local compressive strain reduced below the Biot critical strain at the valley of the wrinkle.

The most notable point on the primary *ruga* phase diagram is the M-point in figure 5. As discussed above, the M-point splits the Biot bifurcation curve into two parts, the unstable instantaneous crease curve and the stable progressive-wrinkle curve. This M-point splitting has an important implication on creasing of neo-Hookean EDM solid surfaces, in general. It defines the range of applied compressive strain 0.272<*ε*<0.456 that can instantaneously crease the surface. The branch curve *C*_{13} in figure 5 implies that instantaneous creasing sets off when *k*/*a* and *ε* have values on the curve. Once the surface creases, the unloading process of the crease is irreversible. The crease configuration recovers the flat state only when the compressive strain is reduced down to the lower bound crease strain shown as a chain curve in figure 2*a*. By taking the limit *a*→0, or , we conclude that the lower bound of the crease strain converges to 0.35 for homogeneous neo-Hookean solids, which is consistent with the findings of Hong *et al*. [9] and Hohlfeld & Mahadevan [18].

Now, we consider what causes the nonlinear bifurcation to be stable on *C*_{12} and unstable on *C*_{13}. Our energy calculation shows that the longer wavelength configuration has a lower energy level among admissible bifurcation modes with respect to the constraint of the boundary condition, typically the specimen size. Therefore, the longest-wavelength mode would be activated first, and the amplitude of the wrinkle would grow. Once the amplitude begins to grow, the wrinkle valley will encounter two competing effects. One is compressive strain relaxation owing to elongation of the cord length of the surface, and the other the local strain concentration caused by the new geometry of the surface. To follow these competing effects, we plot the local strain *ε*_{l} evolution against the applied compressive strain *ε* in figure 6 for isoperiodic compressions with different .

In figure 6, the 45^{°} straight line indicates stable isoperiodic compression keeping the surface flat. In all cases of different , the lines initially turn upwards as they are bifurcated, meaning that every local strain at the valley of the wrinkle is amplified as soon as the surface buckles. The bifurcated curves have finite slopes soon after the bifurcations for ; the steeper the slopes, the larger the is. However, the post-bifurcation slopes immediately diverge to infinity for implying instability. For stable bifurcations of , the strain concentration rates are reduced as they are further compressed. At this stage, the local strain relaxation rate caused by elongation of the surface cord length takes over the strain concentration rate owing to the geometry of the surface. Then, growth of the local strain at the valley of the stable wrinkle accelerates until the local strain reaches close to *ε*_{B}. As Hohlfeld & Mahadevan [18] and Cao & Hutchinson [10] noted in their FEM observations, the plane-strain crease-tip strain is believed to be the Biot critical strain *ε*_{B}. Regarding the curve for in figure 6, the local strain relaxation mechanism prevails, and the specimen can be compressed with its nominal strain much beyond the Biot critical strain limit. In addition, as described earlier, isoperiodic compression makes the progressive wrinkle eventually fold for . Once it folds, as shown in figure 3*e* and marked as (e) in figure 6, the strain concentration rate changes substantially, and creasing is delayed.

Returning to figure 5, the setback-crease curve *C*_{2c} has two peculiar points. One is the T-point (*ε*=0.174, ) at which the applied compressive strain *ε*_{T} takes the absolute minimum value 0.174 for setback creasing. While spontaneous setback creasing is possible in neo-Hookean EDM solids for *ε*_{T}≤*ε*<*ε*_{M}, setback creasing is only possible with initial imperfections in homogeneous neo-Hookean solids. The other is the R-point (*ε*=0.593, ); compression beyond this point folds the wrinkle. The folding is a characteristic of a thin stiff film attached on a soft substrate. This result shows that a compressive strain of more than 0.60 is required to fold a neo-Hookean EDM solid. Once the fold is further compressed, a setback crease sets off at the root of the fold. When we rescale the solid surface at the root of the fold, the renormalized effective compressive strain and the effective admissible minimum wavenumber satisfy the condition to sit on the instantaneous-crease curve *C*_{13} to set off the final setback creasing.

Moreover, it is found that, unlike a soft substrate with a stiff film, a neo-Hookean EDM solid does not activate period multiplication during its wrinkle growth because the critical wavenumber monotonically increases with the compressive strain. Creasing in a neo-Hookean EDM solid hardly develops period multiplication as well, although there is a setback-creasing branch along which the critical wavenumber decreases for increasing compressive strain. It is because setback creasing requires growth of a wavenumber of a preset wrinkle to trigger the unstable snap jump of creasing.

## 5. Conclusions

In this paper, we have explored creasing processes of a neo-Hookean EDM solid surface under isoperiodic plane-strain compression, with results leading to the construction of a primary *ruga* phase diagram. In constructing the phase diagram, our nonlinear bifurcation analysis revealed a critical point, the M-point (*ε*=0.272, ), on the phase diagram. The M-point splits Biot's bifurcation curve [6] such that a flat surface wrinkles stably under a compressive strain in the range 0<*ε*<*ε*_{M}=0.272, whereas it creases unstably for *ε*_{M}<*ε*<*ε*_{B}=0.456. The strain at which the flat surface starts to wrinkle or crease depends on the longest-admissible perturbation wavelength normalized by the decay length of the exponentially graded stiffness.

Moreover, our finite-element post-bifurcation analysis uncovered that the local compressive strain at the valley of the wrinkle is always amplified initially once the flat surface buckles. However, as the compression progresses further, growth of the large wrinkle amplitude tends to relax the local strain due to elongation of the surface cord length. For instantaneous creasing with *ε*_{M}<*ε*<*ε*_{B} and , the initial local strain amplification mechanism dominates to escalate the strain concentration to the level of strain localization of creasing. On the other hand, for stable wrinkle bifurcation with 0<*ε*<*ε*_{M} and , the local strain relaxation mechanism leads to progressive growth of the wrinkle amplitude until further compression sets off setback creasing. In tracing the setback creasing on the *ruga* phase diagram, we also found another critical point, the R-point (*ε*=0.593, , beyond which the wrinkle folds before it creases. In addition, another limit point, the T-point (*ε*=0.174, defines the minimum compressive strain required to trigger setback creasing. The FEM analysis also revealed that regardless of its type, instantaneous or setback, creasing is always accomplished when the local crease-tip compressive strain reaches the Biot strain, *ε*=0.456.

## Acknowledgements

This work was supported by the Korea Institute of Machinery and Materials, the Korea Institute of Science and Technology and the US National Science Foundation (DMR-0520651) through MRSEC at Brown University. The authors acknowledge that Eunice Kim in Classics introduced the Latin word *ruga* for this research.

## Appendix A. Perturbation analysis

**(a) The zeroth-order solution**

The zeroth-order solution is the affine deformation of *F*_{11}=*λ*, *F*_{22}=1/*λ* and *F*_{12}=*F*_{21}=0 with zero traction along the boundary *X*_{2}=0, for which the only non-vanishing stress components are
A1aand
A1b

**(b) The first-order solution**

The first-order incompressibility condition given by equation (2.4c) allows us to relate the first-order displacement components through a function *ψ*^{(1)}(*X*_{1},*X*_{2}) as
A2aand
A2bInserting equations (A 2) into the first-order expansion in (2.4a) and (2.4b), and eliminating the hydrostatic pressure by cross-differentiation, the system of partial differential equations becomes the following fourth-order homogeneous partial differential equation:
A3The general solution of (A 3) is obtained analytically as
A4where the coefficients *C*_{i} are determined only from the boundary conditions. Using equations (A 4) and (A 2), we can now determine the first-order displacement components as given in (2.5a) and (2.5b). The hydrostatic pressure is determined by the tractions in (2.3) with (2.5a) and (2.5b) and the traction-free boundary condition as
A5

**(c) The second-order solution**

Similar to the technique used to solve for the first-order displacement components, the second-order incompressibility condition is used to relate the second-order displacement components as
A6aand
A6bwhere
A6cand
A6dFollowing the same procedure as in the first-order solution, the system of partial differential equations becomes the following fourth-order inhomogeneous partial differential equation:
A7Then, (A7) admits the following solution:
A8where the exponents *β*_{m} (*m*=1,…,4) are obtained from the exponents *α*_{1} and *α*_{2} by replacing by in equation (2.6a). In (A 8), the coefficients *B*_{3} and *B*_{4} are set to zero to ensure the solution decays to zero at infinity.

By inspection of the right-hand side of the differential equation in equation (A 7), we can deduce that the particular solution *ψ*^{(2)}_{p}(*X*_{1},*X*_{2}) has the following form:
A9aThe coefficients *F*_{i} (*i*=1,2,3) are determined by inserting (A 9*a*) into (A 7) followed by matching the coefficients of the different eigenfunctions to obtain
A9band
A9cwhere
A9dwith *υ*_{1}=*λ*^{4}+1, and .

We then insert equation (A 8) into equations (A6a) and (A6b) to determine the second-order displacement components as (2.9a) and (2.9b) and the second-order hydrostatic pressure as (2.9c) given in the main text.

Enforcing the traction-free boundary conditions at the top surface of the half-space with the second-order Piola–Kirchhoff stress, we obtain the following two equations:
A10aand
A10bwhich constitute a system of linear equations for *B*_{1} and *B*_{2} given in appendix B. Normalizing the amplitude of the first-order displacement component to unity in (2.5b), and using (2.7), the constants *C*_{1} and *C*_{2} are determined as
A11aand
A11bThe second-order expansion of the first Piola–Kirchhoff stress is then given as
A12a
A12b
A12c
and
A12d

## Appendix B. *H*_{1}(*X*_{2}), *H*_{2}(*X*_{2}), *B*_{1} and *B*_{2}

B1a B1b B2a B2b and B2cwhere B3a3bwith B4a B4b B4c B4d B4e B4f B4g B4h and B4i

- Received December 28, 2012.
- Accepted June 4, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.