## Abstract

We present the mechanics of folding surface-layer wrinkles on a soft substrate, i.e. inter-touching of neighbouring wrinkle surfaces without forming a cusp. Upon laterally compressing a stiff layer attached on a finite-elastic substrate, certain material nonlinearities trigger a number of bifurcation processes to form multi-mode wrinkle clusters. Some of these clusters eventually develop into folded wrinkles. The first bifurcation of the multi-mode wrinkles is investigated by a perturbation analysis of the surface-layer buckling on a pre-stretched neo-Hookean substrate. The post-buckling equilibrium configurations of the wrinkles are then trailed experimentally and computationally until the wrinkles are folded. The folding process is observed at various stages of wrinkling, by sectioning 20–80 nm thick gold films deposited on a polydimethylsiloxane substrate at a stretch ratio of 2.1. Comparison between the experimental observation and the finite-element analysis shows that the Ogden model deformation of the substrate coupled with asymmetric bending of the film predicts the folding process closely. In contrast, if the bending stiffness of the film is symmetric or the substrate follows the neo-Hookean behaviour, then the wrinkles are hardly folded. The wrinkle folding is applicable to construction of long parallel nano/micro-channels and control of exposing functional surface areas.

## 1. Introduction

Wrinkles of a thin stiff layer attached on a soft substrate, ranging from bio-cellular to geological-crust wrinkles, have been observed widely in nature. The formation process of such wrinkles has recently been recognized as an attractive means of making two-dimensional self-organized patterns at multiple length scales for various applications, such as flexible electronics, optics with pitch-controllable gratings, and nanotechnology of adhesion, wetting and particle-sieving control (Stafford *et al*. 2004; Genzer & Groenewold 2006; Chan *et al*. 2008; Rahmawan *et al.* 2010). Most compression-induced wrinkles are formed by buckling of the surface layer. In such buckling processes, very rich patterns of wrinkle clusters can be produced by extensive material property combinations of the bi-material system (Cerda & Mahadevan 2003; Efimenko *et al*. 2005; Moon *et al*. 2007). Upon further compression, some of these clusters eventually develop into folded wrinkles, i.e. inter-touching of neighbouring wrinkle surfaces without forming a cusp. This wrinkle-folding process provides a mechanical means of concealing and re-exposing partial surface areas, and can be used to construct long parallel nano- and micro-channels. To date, the mechanics of such large-amplitude folded wrinkles and the underlying formation mechanisms are still not well understood, owing to the broad-scale nonlinear deformation characteristics entailed in the wrinkle formation process.

Historically, interest in compression-induced wrinkles was motivated by the buckling of the sandwich panels in aircraft structures (Gough *et al*. 1940; Hoff & Mautner 1945). The buckling was then treated in the context of elastic stability by Timoshenko & Gere (1961) and of incremental elasticity by Biot (1957, 1965). Approximate structural analysis for the buckling of an elastic layer was well summarized by Allen (1969). The simplified structural buckling analysis was extensively used for various applications (Röll 1976; Bowden *et al*. 1998; Stafford *et al*. 2004). A more complete analysis of elasticity for the buckling was provided by Shield *et al*. (1994), which presented the exact formulation of the first-order perturbation analysis for the wrinkle formation with incrementally small-strain deformation of an isotropic substrate. A general incremental bifurcation theory that incorporates the effects of boundary elasticity was presented by Steigmann & Ogden (1997). Lee *et al*. (2008) applied the bifurcation analysis for determining the critical load for the onset of wrinkling and the associated wavelength to substrates in which the elastic modulus is an arbitrary function of depth. Recently, a high-order perturbation analysis was carried out for stiff-layer wrinkles on a neo-Hookean (Rivlin 1948) elastic substrate; however, the analysis was limited to isotropic elastic perturbation with respect to the undeformed configuration of the substrate, and to the single-mode wrinkle deformation with vanishing interface shear tractions (Song *et al*. 2008).

While a variety of patterns have been observed under multi-axial compression of the surface layer (Bowden *et al*. 1998; Chen & Hutchinson 2004; Huang *et al*. 2005; Moon *et al*. 2007), an array of parallel straight wrinkles formed by uniaxial compression are regarded as generic wrinkles. Multi-scale clustering of small-amplitude generic wrinkles has been observed experimentally on a UV-cured polydimethylsiloxane (PDMS) substrate (Efimenko *et al*. 2005; Huck 2005). The UV treatment of PDMS produces a graded distribution of elastic stiffness near the surface and the multi-mode wrinkles of small amplitude may be attributed to graded stiffness effect. On the other hand, a thin film layer of distinct uniform property on a soft substrate also produces multimode wrinkles as the amplitude grows large. The multi-mode generic wrinkles of large amplitude are generated by multiple bifurcations of the equilibrium wrinkle configuration. Formation of such multi-mode generic wrinkles of large amplitude is needed to fold wrinkles. However, the formation of the multi-mode generic wrinkles is insufficient to fold wrinkles, and only some of the multi-mode generic wrinkles develop into folded wrinkles as the amplitude grows larger. As will be shown in this paper, certain material nonlinearities in large deformation kinematics are required to trigger the wrinkle folding.

Large-amplitude wrinkles of a stiff surface layer are formed on a soft substrate generally pre-strained in finite deformation. The finite deformation of the substrate is highly nonlinear and induces elastic anisotropy for incremental deformation. The large-amplitude wrinkling often bends the surface layer beyond its elastic limit. In this paper, we address the effects of such nonlinear deformation on the wrinkle formation from its inception to final folding. Figure 1 shows a series of schematics for the wrinkling process of an initially un-stressed thin stiff layer on a pre-stretched soft substrate. In these schematics, a long strip of a soft polymeric (PDMS) substrate of width *W* (figure 1*a*) is pre-stretched in the width direction to a stretch ratio *λ*_{ps} (figure 1*b*) in plane-strain deformation. The pre-stretching is followed by deposition of a thin metallic (gold) film of thickness *h* on the surface (figure 1*c*). Then, the substrate is released to a stretch ratio of *λ*<*λ*_{ps}, inducing thin film buckling to form post-buckling wrinkles (figure 1*d*). In our experimental study reported in §3, we measure the wavelengths and shapes of the wrinkle patterns, at various levels of *λ*, with atomic force microscope (AFM) topography and scanning electron microscope (SEM) imaging of the wrinkle cross sections exposed by focused ion beam (FIB) sectioning. In the experiment, we employ the notion of nominal compressive strain *ε*^{C}_{L}=(*λ*_{ps}−*λ*)/*λ*_{ps} applied to the surface layer by variation of the substrate stretch from the pre-stretch *λ*_{ps} to a released state of stretch *λ*.

## 2. Bifurcation of a generic wrinkle on a neo-Hookean substrate

Buckling of a thin stiff layer on a pre-stretched soft substrate is analysed in this section. We show how eigenmodes of the substrate deformation are involved in buckling the surface layer under compression to form bifurcated wrinkle patterns with two different wavelengths, and how the bifurcated states evolve as the film is further compressed by shrinking the pre-stretched substrate. The soft substrate is typically made of a non-porous rubbery material such as PDMS which exhibits characteristics of incompressible finite deformation at room temperature (Yoo *et al*. 2006). The nature of such incompressible finite deformation is well epitomized by the neo-Hookean constitutive relation (Rivlin 1948), for example, for uniaxial tension of a typical PDMS in the range of stretch ratio from 0.7 to 1.5. The eigenmodes of perturbed substrate deformations in a uniformly stretched neo-Hookean nonlinear finite deformation are analysed first in the following.

### (a) Analysis of the substrate deformation

We consider a plane-strain deformation of the substrate for the half space of and in which vertical surface displacement is perturbed by on a pre-stretched configuration of lateral stretch ∂*x*_{1}/∂*X*_{1}=*λ* and vertical stretch ∂*x*_{2}/∂*X*_{2}=1/*λ*. Here, the wrinkle wave number *k*=2*π*/*L* and the perturbation amplitude parameter *δ*≪1 are employed with *L* the undeformed-configurational wavelength of the perturbation. Then, the generic form of the mapping of the configuration *x*_{2}(*X*_{1},*X*_{2}) is assumed to be
2.1
where *α*>0 is a parameter to be determined from the equilibrium condition. The incompressibility condition, up to *O*(*δ*^{2}), requires the mapping *x*_{1}(*X*_{1},*X*_{2}) as
2.2
Here, we assume these functional forms of perturbation since the finite pre-stretching of the substrate generates an incrementally anisotropic elastic state (Steigmann & Ogden 1997) and the anisotropic linear elastic field in a two-dimensional half space loaded by a sinusoidal traction can be expressed by superposition of two elastic fields that are the same as the assumed perturbation fields with two distinct *α* values (A. C. Pipkin 1978, unpublished data). Furthermore, as the hydrostatic pressure for incompressible neo-Hookean solids has to be determined independently from the deformation kinematics, the pressure distribution is assumed to be
2.3
with *p*_{0} and *p*_{1} to be determined from boundary conditions.

The first Piola–Kirchhoff stress is given by
2.4a
for the incompressible neo-Hookean deformation energy potential of
2.4b
where *λ*_{i} are the principal stretches, *Q*=*E*_{0}/3 with *E*_{0} the elastic modulus of the neo-Hookean solid at *λ*_{i}=1, *F*_{ij}=∂*x*_{i}/∂*X*_{j} the deformation gradient, and the Jacobian *J*=det[*F*_{ij}]=1 for incompressibility. Then, the deformation gradient obtained explicitly from (2.1) and (2.2) is inserted into (2.4a) to have
2.5a
2.5b
2.5c
2.5d
The zeroth order (*δ*=0) solution is then obtained from (2.5) and as *p*_{0}=*Q*/*λ*^{2} and
2.6
With the above stress representations, the equilibrium equations up to the first order are expressed as
2.7a
and
2.7b
Then, equation (2.7) implies an eigenvalue problem of *α* as
2.8
The positive eigenvalues, *α*=1 and 1/*λ*^{2}, of equation (2.8) for non-diverging solution and the corresponding eigenvector solutions provide *p*_{1}=0 and *p*_{1}=2*π*(1/*λ*^{4}−1)*λQ*, respectively. The two distinct eigenvalues, for *λ*≠1, stand for two distinct eigenstates of the perturbed deformation on the incrementally anisotropic elastic state generated by the pre-stretched finite deformation. In neo-Hookean solids the incremental modulus softens in the stretching direction while it stiffens in the shrinking direction. For the degenerate case of *λ*=1, the two eigenstates are rearranged by to those of decaying in −*X*_{2} direction with e^{kX2} and *X*_{2}e^{kX2}, respectively, which have been used in conventional analysis of isotropic linear elasticity (A. C. Pipkin 1978, unpublished data). The two eigensolutions provide the first-order functional forms of the first Piola–Kirchhoff stress components as
2.9a
2.9b
2.9c
2.9d
where *A* and *B* are coefficients to be determined by the boundary conditions at the interface between the surface layer and the substrate. The Cauchy stress components *σ*_{ij}=*J*^{−1}*F*_{ik}*T*_{kj}, then, provide the expression of the traction components at the interface as
2.10a
and
2.10b
Corresponding displacement components are obtained from (2.1), (2.2) and (2.8) as
2.11a
and
2.11b

### (b) Plate-theory analysis of the surface-layer deformation

The expressions for the tractions and the displacements at the interface in (2.10) and (2.11) have to be compatible with the deformation of the surface layer. Although the compatibility relationship can be expressed exactly up to the accuracy of the linear elastic description of the layer deformation (Shield *et al*. 1994), the simple theory of plate stretching and bending with single-axis curvature is accurate enough and more conveniently applicable to layer deformations of asymmetric bending to be discussed in the following sections. The compatibility relationships of the interface tractions and displacements with respect to the linearized layer deformation consistent with the first-order perturbation of the substrate deformation are expressed as (Shield *et al*. 1994):
2.12a
and
2.12b
where is the cross-sectional average strain of axial compression in the layer, and *E*_{L} the Young’s modulus and *ν*_{L} the Poisson’s ratio of the surface layer. Then, as (2.10) and (2.11) are inserted into (2.12), we will have a highly nonlinear eigenvalue problem of the normalized wave number for buckling as
2.13
with and , where the characteristic equation gives the buckling strain *ε*_{B}
2.14

Since the stretch ratio *λ* of the substrate and the stiffness ratio between the layer and the substrate are positive, all the coefficients of the polynomials in the numerator of (2.14) are also positive, and thus the numerator is a monotonically increasing fourth-order positive polynomial function of , for . The denominator of (2.14) makes *ε*_{B} singular at and , having a single positive minimum in the domain of . The minimum *ε*_{B} is the compressive strain at the onset of buckling, *ε*^{cr}_{B}, and the corresponding determines the onset wave number of the wrinkle . By making the derivative of (2.14) with respect to vanish, i.e.
2.15
we have the unique positive root of (2.15) as the critical buckling (minimum) wave number .

Figure 2*a* shows plots of (2.15) for the substrate stretch ratio *λ* at which the surface layer critically buckles with the wave number for various stiffness-ratio values. Equations (2.14) and (2.15) asymptotically converges, for a very high stiffness surface layer on a very soft substrate, to
2.16a
and
2.16b
Equations (2.16a) and (2.16b) show that the critical compressive strain for the surface layer to be critically buckled at the stretch ratio *λ* of the substrate scales with {(1+*λ*^{2})/2}^{2/3}, and the wavelength in the deformed configuration (*L**=*λL*) of the wrinkle with {(1+*λ*^{2})/2}^{−1/3}. As an example, the roots of equations (2.14) (solid line) and (2.15) (dashed line) are plotted in figure 2*b* for elastic buckling of a porous gold film with nominal Young’s modulus of 20 GPa and Poisson’s ratio of 0.3, attached on a PDMS substrate. Since we use PDMS in our experiment, the deformation characteristics of which is close to those of the neo-Hookean solid with *Q*=0.15 MPa approximately up to *λ*=1.5, we plot the wavelength *L** in the deformed configuration as a function of the nominal compressive strain *ε*^{C}_{L}=(*λ*_{ps}−*λ*)/*λ*_{ps} near the buckling point. The dashed line shows the roots of (2.15), which stands for the trace of the critical buckling points for different *λ*. The critical buckling wavelength looks constant in figure 2*b*; however, it is slightly increasing for decreasing *λ* (or increasing *ε*^{C}_{L}=(*λ*_{ps}−*λ*)/*λ*_{ps}). On the other hand, the solid curve indicates the roots of (2.14), meaning the trajectory of the buckling wavelength versus the nominal compressive strain during release of the substrate strain from the prestretched state of *λ*_{ps}=1.5. As shown by the solid curve, the buckling state bifurcates beyond the critical nominal compressive strain of 0.065 per cent. The wavelength of the primary buckling mode decreases as the nominal compressive strain increases, while the long wavelength of the secondary buckling mode grows rapidly as the compressive strain increases. As the wavelength of the secondary mode diverges so quickly that most of the laboratory observations in a single microscopic scale are limited to the primary buckling mode. However, an observation of the bifurcation phenomenon in computational simulation of a thin layer buckling on a neo-Hookean solid substrate was reported by Efimenko *et al*. (2005) and Huck (2005). Such a bifurcation was also observed in experiments by Shield *et al*. (1994). As shown in figure 2*b*, the wavelength of the primary wrinkle quickly reduces as soon as the layer buckles and strains a little further within a fraction of the critical compressive strain. Therefore, the measurement of the critical buckling wavelength may not be accurate, if we rely only on direct experimental observations in using (2.16b) to measure the stiffness of a surface layer. However, the prediction of the post-buckling wavelength by the first-order perturbation has limited accuracy beyond the buckling point. High-order perturbation analysis is expected to improve the accuracy of the evaluation of the post-buckling wavelength. In the following sections, evolution of the post-buckling large-amplitude wrinkles towards folding is investigated by the hybrid method of experiment and finite-element analysis.

## 3. Experiments on evolution of large-amplitude wrinkles

The post-buckling wrinkle patterns were generated experimentally in a thin gold film attached on a PDMS substrate by the following procedures. We prepared the PDMS by mixing dimethylsiloxane monomers with Sylgrad 184 silicone elastomer curing agent in a weight ratio of 15 : 1, followed by degassing in a vacuum chamber and curing at 80^{°}C over 1 h. The PDMS samples were prepared in the size of 40.0×10.0×3.0 mm. As depicted in figure 1*b*, the PDMS sample was pre-stretched to *λ*_{ps}=2.1, using a screw-driven tensile device which could apply the nominal compressive strain, *ε*^{C}_{L}=(*λ*_{ps}−*λ*)/*λ*_{ps}, more than 50 per cent in the resolution of 0.5 per cent. The gold film was then ion-coated on the pre-strained PDMS using a plasma-ion coater (IB-3, Eiko) with 6 mA current, as illustrated in figure 1*c*. The transmission electron microscope (TEM) image shown in figure 5*a* reveals that the gold film was porous in the nanometre scale and the grain size was in the range 10–20 nm. In our experiments, we used samples with various film thicknesses ranging from 20 to 80 nm. The gold film coated on the pre-stretched PDMS was compressed as the pre-stretched strain was relaxed, generating periodic wrinkles grooved perpendicular to the compression direction. The wrinkled morphology was observed and measured *in situ* with optical and AFM (Digital Instrument), while the pre-stretch of the PDMS was released, as shown in figure 1*d*. After the pre-stretch was completely released, the PDMS was re-stretched back to the original pre-stretch, during which the gold film returns to the strain-free state. While the gold film was compressed and the amplitude of the wrinkles grew, the cross-sectional shapes of the wrinkles were observed by sectioning the film with FIB (NOVA200, FEI) milling. The cross-sectional shapes were recorded by the high-resolution scanning electron microscope (HR-SEM) built in the FIB-system, at four different stages of the nominal compressive strains, 12, 35, 41 and 46 per cent. Before being sectioned by the FIB, the surface of the wrinkled thin film was pre-coated with carbon in the thickness range 200–500 nm, followed by a platinum coating of 500–2500 nm thickness to hold the shape of the wrinkle from relaxation and to prevent potential damages of the film from gallium ion irradiation.

### (a) Observations on evolution of multi-mode wrinkles

The optical-image sequence (*a*1)–(*a*6) in figure 3 shows evolving wrinkles of a 66 nm thick gold film during a compression–decompression cycle. As the pre-stretch of the PDMS was released a wrinkle pattern emerged as shown in (*a*2), implying that the film already buckled at a nominal compressive strain *ε*^{C}_{L} less than 0.5 per cent. Non-uniformity of the wrinkles in (*a*2) indicates that nucleation sites of buckling were randomly distributed and incoherent phases of the buckles were merging together at the small nominal compressive strain, presumably owing to fluctuations of the structural properties such as film thickness, roughness of the substrate and variations of the stiffness of the film and the substrate. Once the nominal compressive strain was increased to 7 per cent, all the pitches of the wrinkles were coherently aligned as shown in (*a*3). At this nominal compressive strain two bifurcated modes of wrinkles were observed; the primary-mode wrinkles of 2.5 μm wavelength were nested in the enveloping secondary-mode wrinkles of 30–40 μm wavelength. Such bifurcation was predicted by the analysis in §2, and also observed experimentally and in computational simulations by Huck (2005) and Efimenko *et al*. (2005). The aspect ratio of the PDMS substrate width in stretching direction to its lateral span was too large to produce plane-strain deformation, and the Poisson expansion of the 7 per cent compression cracked the film at some places (not shown in the frame). In (*a*4), the compressive strain was 40 per cent and the primary-mode wrinkles of a single-period evolved into paired wrinkles of dual period. As the film was further compressed to 50 per cent, secondary pairing of paired wrinkles occurred, activating the quadruple-wavelength mode as shown in (*a*5). After this point, the compression of the film was released by re-stretching the PDMS substrate. When the applied compression was completely removed, the film was not fully recovered to its original flatness. Instead, plastically deformed wrinkle patterns remained as shown in (*a*6). The plastically deformed wrinkle patterns could not be totally removed even after the film was further stretched to 10 per cent tensile strain.

The three frames of figure 3(*b*1), (*b*2) and (*b*3) show triple, double and quadruple wrinkle pairs observed in gold films of 23, 30 and 44 nm thicknesses, respectively, at nominal compressive strains in the range 40–50%. In the *in situ* AFM measurements of wrinkle topography for a 30 nm thick gold film, a single-period wavy pattern of 1.3 μm wavelength was distinctly observed at 5 per cent compressive strain, while the paired wrinkles of dual period emerged at 41 per cent compressive strain. Fast Fourier transformation of the AFM topography showed a clear spectrum split of the wavelength into 0.9 and 1.8 μm for the paired wrinkles. At this point, the total arc length of the wrinkle measured by the AFM topography was reduced substantially. The reduction of the AFM-traced arc length implied folding of the wrinkle which did not allow access of the AFM tip to trace the arc length at the folded region. This observation prompted us to section-wrinkled films with FIB at different stages of compression. Upon complete removal of the applied nominal compressive strain in the film, the residual deformation of the paired wavy pattern remained permanently owing to irreversible plastic deformation.

Figure 4*a* shows the cross-sectional shapes of a buckled gold film of 74 nm thickness cut along the loading axis by FIB at four different places for four different levels of the nominal compressive strain (12, 35, 41 and 46%). The exposed cross sections were imaged by SEM at an angle of 38^{°} tilt with respect to the normal of the cross-section, towards the surface normal of the undeformed PDMS substrate. Therefore, actual amplitudes of the wrinkles are larger than those shown in the figure by factor of 1/cos 38^{°}. As shown in (i) of figure 4*a*, the primary wrinkle at the nominal compressive strain of 12 per cent exhibits single-period waviness. However, the curvature at the valley of the wrinkle is considerably larger than that at the peak, insinuating that the film might have asymmetric bending stiffness. This observation brought us about imaging the cross-sectional configuration of the nano-grains in the gold film with a TEM, which is shown later in figure 5*a*. As the compressive strain reached 35 per cent, the curvature contrast was even enhanced as shown in (ii), and repeated pairs of the wrinkles begin to bunch, activating the paired wrinkle pattern of dual period. A close look reveals that the platinum coating layer applied prior to FIB sectioning was cracked open in every other valley of the wrinkle. Implication of this phenomenon will be treated in §4 in comparison with computational results. The cross-sectional configuration of the wrinkles formed at the nominal compressive strain of 41 per cent, shown in (iii), unveils that indeed the bunching pairs of wrinkles touched each other to create periodic folding of wrinkles, for which the AFM tip could not follow the whole cord of the film. As the nominal compressive strain was further increased up to 46 per cent in (iv), the energy-release capacity of the primary and secondary bending modes in the folded film with dual period seemed to be saturated, and the increment of the compressive strain was rather accommodated by bending the film with quadruple-wavelength period. In this stage, every other fold of the film was opened by the period-doubling mode of deformation.

Figure 4*b* shows the variation of the wrinkle wavelength in the 74 nm thick film with respect to the nominal compressive strain. The variation was monitored by *in situ* AFM measurements. The wavelength variation nearly followed the substrate shrinkage rate, *L**≃*L*(1−*ε*^{C}_{L}), implying that the primary wavelength of the wrinkle in the undeformed configuration might be locked by plastic bending. When the nominal compressive strain reached 35 per cent, the paired wrinkles with dual period were produced; however, this dual-period mode of deformation did not affect the variation trend of the primary period. In the same manner, when the quadruple-period mode was set off, the wavelength variations of the primary mode seemed to be hardly influenced, strongly suggesting that the primary wavelength was set by plastic bending at an early stage of wrinkle development.

### (b) Material characterizations

As discussed above, nonlinear deformation characteristics as well as non-uniformity and anisotropy of the substrate stiffness are critically related to clustering and folding of surface-layer wrinkles. In addition, asymmetric bending stiffness and/or asymmetric yield moment of the surface layer produce bending configurations with high contrast of curvatures at the valley and the peak of the wrinkle, which seem to promote folding of the wrinkle. The cause of such asymmetry in the gold film with nanometre scale thickness was found to be an asymmetric distribution of grain grooves on the top and bottom surfaces. The cross section of a gold film of nominally 50 nm thickness, imaged by an S-TEM and shown in figure 5*a*, was porous and had asymmetric distribution of grain grooves. It had more and deeper grain grooves on the gold/PDMS interface side, presumably owing to growth history of the nano-grains during the deposition process. When the film is bent to close the deep grain grooves the effective bending stiffness is relatively higher than the bending stiffness to open the grooves, since the groove surfaces contact each other as they are closed. In addition, the bending to open the grooves will readily induce plastic yielding.

A schematic of the bending moment–curvature relationship is depicted in figure 5*b* as a dashed curve. However, direct experimental measurement of the moment–curvature relationship is difficult for such a thin film. We employed a hybrid method of experiment and computation to evaluate the bending characteristics by matching the experimentally observed shapes of the wrinkles with computationally simulated ones. The effective Young’s modulus *E*_{L} of the porous gold film at the onset of buckling was determined as 20 GPa, and the Poisson ratio *ν*_{L} was adopted as 0.3, by matching the shapes of relatively small-amplitude wrinkles. On the other hand, the effective bending stiffness for large bending curvature was set as , where *ξ*^{±} is the asymmetry factor of bending, and the superscripts ± stand for bending in the directions of closing and opening the interface grain grooves, respectively. We identified *ξ*^{+}=1 and *ξ*^{−}=0.1 for the ion-deposited gold film, with finite-element simulations.

Figure 5*c* shows the experimental uniaxial stress–strain curve of the PDMS substrate, together with the best fits of two hyperelasticity models. The second-order Ogden (1972, 1998) model curve is fitted to the experimental data with a Poisson’s ratio of 0.48. For the incompressible neo-Hookean model, the experimental data points of nominal strain between 0 and 0.5 are used for the data fitting. Both models are used in the present study. For the neo-Hookean model, we represented the hyperelastic energy potential as expressed in (2.4b). The stress–strain curve of a PDMS uniaxial tension test closely matches the neo-Hookean curve up to 50 per cent stretching with *Q*=0.1487 MPa, as shown in figure 5*c*. However, the stress–strain curve substantially deviates from the neo-Hookean curve beyond the nominal tensile strain of 0.5, and we employed the Ogden (1972, 1998) model to match the deformation behaviour all the way to the experimentally adopted value of pre-stretch, *λ*_{ps}=2.1. The general form of the Ogden (1998) strain energy potential is given by
3.1
where *J*=*λ*_{1}*λ*_{2}*λ*_{3} is the volume ratio, are the deviatoric stretches, and *μ*_{i}, *α*_{i} and *D* are material parameters. By taking *N*=2 and fitting the Ogden potential to the uniaxial stress–strain curve over the whole strain range, we obtain *μ*_{1}=0.0001440 MPa, *α*_{1}=14.24, *D*=1.263 MPa^{−1}, *μ*_{2}=0.2674 MPa and *α*_{2}=2.743.

## 4. Computational simulations of wrinkle folding: finite-element analysis

While the bifurcation of wavelengths at the first buckling point and the post-buckling evolution of relatively small amplitude wrinkles could be understood with the perturbation analysis described in §2, we analyse folding of highly nonlinear large-amplitude wrinkles of gold thin films on a PDMS substrate with finite-element method in this section. The finite-element simulations are carried out as a two-dimensional plane strain problem in ABAQUS/Standard using the stabilized Newton–Raphson method. Figure 6 shows a schematic of the finite-element model in its initial configuration. The film is modelled using 2-node linear two-dimensional beam (B21) elements, and the substrate is meshed with 4-node bilinear generalized plane-strain quadrilateral (CPEG4) elements. In the simulations, the substrate is first pre-stretched to *λ*_{ps}=2.1, with respect to the original substrate width. In the mean time, the film undergoes an artificial thermal expansion in the lateral direction, in order to maintain stress-free state while deforming together with the substrate. Then, the film and the substrate are compressed together as done for the experiments. During the first 0.5 per cent compressive strain, the film is guided by displacement boundary conditions to deform it into a buckling mode. After the system is relaxed by removing the boundary conditions on the film, the system is further compressed to desired compressive strain of *ε*^{C}_{L}. During compressive loading, small damping forces are added to the global equilibrium equations for stabilizing the wrinkling process. When the desired compressive strain is reached, an additional relaxation step is performed by removing the damping forces to ensure static equilibrium.

The material properties given in §3*b* are used for the finite-element analysis. For the simulations with incompressible neo-Hookean model, the material incompressibility of the substrate is approximated by setting a Poisson’s ratio of 0.4999, since perfect incompressibility is not supported by CPEG4 elements.

### (a) Effects of substrate material nonlinearity

Beyond the first bifurcation point for neo-Hookean substrates, we searched for the second bifurcation possibility of the primary wrinkle when the amplitude becomes larger under further compression. However, in our limited computational search range of 1≤*λ*≤*λ*_{ps}=2.1 for , the single mode of the primary wrinkle persisted and the second bifurcation was not observed in the model of a symmetric bending-stiffness layer on the neo-Hookean substrate. On the other hand, the second bifurcation occurred at a nominal compressive strain *ε*^{C}_{L} between 12 and 35 per cent in the symmetric bending-stiffness layer on the second-order Ogden model substrate. However, when the film was further compressed beyond 41 per cent nominal strain, the bi-modal primary wrinkle created by the second transition returned back to the single mode primary wrinkle as shown in figure 7(Og4*).

Noticing that the substrate material nonlinearity could induce high-order bifurcations of the primary wrinkles, we compared the wrinkle shapes of the same asymmetric-stiffness layers on neo-Hookean and Ogden model substrates. Figure 7(nH1), (nH2) and (nH3) is the result for the neo-Hookean model, while (Og1), (Og2) and (Og3) for the Ogden model. These results were also compared with experiments. In close examination of the experimental figures, it was found that the nominal compressive strain imposed by the experimental apparatus was different from the local values evaluated by measuring the cord length of the wrinkle in the cross-sectional view. While the imposed nominal compressive strains were 12, 35 and 41 per cent for the cross-sectional images shown in (Ex1), (Ex2) and (Ex3), the local measurements gave 17, 42 and 47 per cent, respectively. We used the local nominal compressive strains for the simulations. Although the cross-sectional shapes of the single-mode primary wrinkles of the neo-Hookean and the Ogden models were close to each other at the 17 per cent compression as shown in (nH1) and (Og1), the stress-concentration patterns were very different. The cross-sectional configurations became substantially different as well at the 42 per cent compression. The second bifurcation transition was just starting in the neo-Hookean model shown in (nH1); on the other hand, (Og2) indicates that the transition in the Ogden model occurred at a lower compression level. Finally, the Ogden model in (Og3) shows folding of wrinkles at 47 per cent compression as observed in (Ex3), while the neo-Hookean model in (nH3) did not.

### (b) Effects of asymmetric bending: bending nonlinearity of the surface layer

As mentioned above, a layer of symmetric bending stiffness on a neo-Hookean substrate hardly produced the second bifurcation transition. Although the layer of symmetric bending stiffness attached on an Ogden substrate generated the second transition, it did not show any folding of wrinkles. It did not exhibit the third transition to the quadruple-wavelength mode observed in the experiment. Furthermore, the maximum bending curvature was observed at different locations in the layer. The symmetric bending stiffness of the layer always gave the maximum curvature at the peak of the wrinkle, while the maximum curvature was spotted at the valley of the wrinkles in the experiment. As mentioned earlier, these observations directed us to consider asymmetric bending stiffness of the gold films. When we tuned the asymmetric bending stiffness factor *ξ*^{−} down to 0.1 from unity, we could match the third transition to the quadruple mode at , as shown in figure 7(Og4^{**}). In this simulation, we allowed inter-penetration of folding surfaces of the wrinkles. With the asymmetric bending stiffness factor *ξ*^{−}=0.1, not only the third transition mode but also the maximum curvatures at the right locations at different compression stages could be produced. Nevertheless, still some subtle differences between the experiment and simulation were observed. In figure 7(Og2), every other valley was much deeper than those observed in the experiment (Ex2). However, close examination reveals that the platinum coating in every other valley was cracked in the experiment, and the SEM images of the valleys far away from the FIB sectioned region indicate that the valleys with coating cracks were deeper than the valleys without coating cracks. Presumably, the stress relaxation during FIB sectioning open the coating crack near the sectioned region. The simulations were carried out by searching for the minimum energy configurations, floating the primary wavelength with symmetric boundary conditions for the simulation zone size twice the primary wavelength, except for (Og4^{**}). The symmetry boundary conditions were imposed for the simulation region of (Og4^{**}) with quadruple length of the primary wavelength.

When a stiff layer is buckled on a compliant substrate, delamination of interface may occur even without a noticeable growth of wrinkle amplitude if the interface bonding is weak (Shield *et al*. 1994). Therefore, we investigated the stress state produced by a folded large-amplitude wrinkle. Figure 8 shows the simulated wrinkling profile compared with experiment (figure 8*a*) and the von Mises stress distribution (figure 8*b*) in (Og3) of figure 7. The von Mises stress distribution was limited to 4.54 MPa. Figure 8*c* shows that the compressive interface traction at the bottom of the closed valley was very large down to −9.46 MPa, while the tension was limited to 0.87 MPa. The maximum shear traction was 2.53 MPa. The results imply that the interface is well protected from tensile delamination in large amplitude wrinkles on a soft substrate such as PDMS.

## 5. Discussion and conclusions

As introduced in §2, the bifurcation of the post-buckling equilibrium wrinkle configuration on a pre-stretched neo-Hookean substrate involves four eigenstates of incremental substrate deformation. The incremental elastic anisotropy of the pre-stretched substrate determines the −*X*_{2}-direction decay rate of the eigenstate. The amplitude sharing between the eigenstates is decided by the interface boundary conditions compatible with the deformation of the surface layer. The interface compatibility condition provides the bifurcation eigenvalue problem for the wave number of the sinusoidal buckling. Similar to the solution procedure of the first bifurcation point provided by the first-order perturbation analysis in §2, high-order asymptotic analysis is believed to be possible to follow the post-bifurcation equilibrium paths for the growth of wrinkles. The asymptotic analysis is also considered to be possibly extended beyond the neo-Hookean model of the substrate and the symmetric-bending-stiffness model of the surface layer.

When we followed the post-buckling equilibrium configurations of the primary wrinkles in 20–80 nm thin gold films attached on a pre-stretched PDMS substrate, we experimentally observed the second and the third bifurcation points within the stretch ratio of 1≤*λ*≤*λ*_{ps}=2.1. The second and the third bifurcations generated periodic clustering or bunching of wrinkles. FIB sectioning revealed that the clustering led to folding of the wrinkles when the surface layer was highly compressed. While multi-scale clustering of small-amplitude wrinkles may be possible on a substrate with graded material properties, folding of the wrinkles requires multiple symmetry-breaking nonlinear material characteristics of both the substrate and the surface layer. Our computational tracing of the equilibrium wrinkle configurations revealed that the neo-Hookean substrate could not fold the wrinkles regardless of the elastic surface-layer nonlinearity. In contrast, the Ogden substrate could fold the wrinkles as observed in the experiments, if the surface-layer stiffness was asymmetric in bending. However, when the bending stiffness of the surface layer was symmetric, the Ogden model could not produce the folded wrinkles either.

The folding of the wrinkles is believed to be retarded or promoted by inelastic behaviour of the film and the substrate. Viscoelastic behaviour of the substrate (Huang 2005) may provide time-dependent folding behaviour of the wrinkles. Plastic bending and unbending processes of the surface layer (Kim & Kim 1988) can also generate different wrinkle patterns than predicted by the elastic analysis. Our *in situ* AFM measurement of the wrinkle wavelengths strongly suggests that the phase of the wrinkle curvature was locked in the undeformed configuration of the film. The wavelength of the primary wrinkle conformally followed the stretch of the substrate when the gold film was compressed to be folded. In addition, the film flatness was not recovered when it was fully relaxed; instead, plastically deformed wrinkle grooves remained. Once the wrinkle phase is locked in the undeformed configuration by the surface-layer plasticity, incremental configuration forces (Kim *et al*. 2010) caused by the external compression loading may unlock the phase by unbending the plastically deformed surface layer. Further studies on inelastic effects of the materials on wrinkle folding are in progress. The multiple folds of thin films with nanometre-scale thickness are expected to have extensive applications such as nano-fluidics in the parallel nano-channels of the folding wrinkles, self-cleaning adhesion control with concealing and re-exposing surfaces, and creation of nanometre-scale solid reaction surfaces for catalysts and battery electrodes.

## Acknowledgements

S.X. and K.S.K. were supported by the U.S. National Science Foundation (DMR-0520651) through MRSEC at Brown University, and by Korea Institute of Science and Technology (KIST); M.W.M. by KIST and in part by the Ministry of Knowledge Economy of Korea through the Fundamental R&D Program for Core Technology of Materials; J.Y.S. and K.H.O. were supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (R11-2005-065). We also thank M. Diab for valuable discussions and the late H. M. Rho for his support for the hybrid computational analysis.

- Received September 17, 2011.
- Accepted October 26, 2011.

- This journal is © 2011 The Royal Society