Direct numerical simulations of three-dimensional surface instability patterns in thin film-compliant substrate structures

A comprehensive numerical study of three-dimensional surface instability patterns is presented. The formation of wrinkles is a consequence of deformation instability when a thin film, bonded to a compliant substrate, is subject to in-plane compressive loading. We apply a recently developed computational approach to directly simulate complex surface wrinkling from pre-instability to post-instability in a straightforward manner, covering the entire biaxial loading spectrum from pure uniaxial to pure equi-biaxial compression. The simulations use embedded imperfections with perturbed material properties at the film-substrate interface. This approach not only triggers the first bifurcation mode but also activates subsequent post-buckling states, thus capable of predicting the temporal evolution of wrinkle patterns in one simulation run. The state of biaxiality is found to influence the surface pattern significantly, and each bifurcation mode can be traced back to certain abrupt changes in the overall load–displacement response. Our systematic study reveals how the loading condition dictates the formation of various instability modes including one-dimensional (1D) sinusoidal wrinkles, herringbone, labyrinth, and checkerboard.

current paper presents our first comprehensive study on surface patterns encompassing the full range of inplane compressive loading. Instead of treating surface patterns discretely under specific assumptions, we aim to provide a full picture by applying our numerical approach to directly simulate any possible forms of wrinkling instabilities and their transition.
It is worth mentioning that the terms "direct numerical simulations" and "instability" used in this paper should not be mistaken as similar terminologies in the fluid mechanics literature. Here the use of direct numerical simulations is to distinguish between the common multi-step approaches in the solid mechanics literature and our embedded imperfection technique which only requires a single-step analysis. Furthermore, instability in solid structures under compression is frequently referred to as the onset of buckling, as opposed to its definition in fluid mechanics as the onset of development of turbulence. Nevertheless, apparent similarities exist and one may draw analogies between the present surface instability analysis in solid materials and the studies of receptivity 39 or non-modal disturbance growth 40 in fluid mechanics.

Brief overview of theories
Various theoretical formulations of surface instability are available in the literature. Here we include representative analytical solutions, some of which are used for numerical model verifications in the current study. Detailed discussion of the theories is given in the Supplementary Information of this paper. For surface wrinkling of a thin film on top of a compliant substrate, one may categorize available analytical techniques into generic groups of 1D 2,13-16 and 2D buckling formulations [17][18][19][20][21][22][23]41 . The 1D solutions, based on the plane strain assumption, focus on the classical sinusoidal wrinkles (also termed 1D wrinkles). The 2D analytical solutions allow for in-plane biaxial loading.
Consider the thin film-substrate structure shown in Fig. 1. Under uniaxial compression, 1D sinusoidal wrinkling commences once the critical point for instability is reached. Assuming that the substrate is semi-infinite, film and substrate are fully bonded at the interface, and both layers are linear-elastic and isotropic, the wavelength of the wrinkles at the onset of bifurcation (primary instability) follows where cr is the critical wavelength of the 1D mode, and t f , E f , and ν f are, respectively, thickness, Young's modulus, and Poisson's ratio of the film layer. The parameter E s = E s / 1 − υ 2 s , where ν s and E s are Poisson's ratio and Young's modulus, respectively, of the substrate. An alternative expression for the parameter E s is given 17 www.nature.com/scientificreports/ Note that it is equivalent to the critical stress in Eq. (2) divided by the plane-strain modulus of the film, For a film-substrate structure subjected to pure equi-biaxial compression, the square-checkerboard pattern has been shown to possess the lowest energy in the buckled state 17,19 . The wavelength of square-checkerboard was derived 17,18,20 as It was also postulated that the critical stress introduced in Eq. (2) not only applies to 1D wrinkles but also applies to any possible biaxial wrinkling mode. The critical strain for the square-checkerboard mode, (e cr ) CB , is apparently derived 17,18,20 by dividing σ cr by the biaxial modulus of the film, In addition, the amplitude of the surface wrinkles, A , can be written in a general form as where the parameter is a function of Poisson's ratio and thickness of the film layer, with = t f for the sinusoidal 1D mode and = t f · 8/[ 3 − υ f 1 + υ f ] for the square-checkerboard mode 17,20 , and e/e cr is the applied compressive strain normalized by the critical value at the onset of primary bifurcation. Controversies exist in the literature regarding the wrinkling patterns pertaining to specific loading states. There are other well recognized surface instability patterns such as herringbone (or zigzag) and labyrinth, which emerge after the primary modes when the compressive stress is well beyond the first critical point. There are also uncertainties regarding their evolution and necessary loading condition in the post-instability regime. A detailed discussion is given in the Supplementary Information. On the experimental side, although the most energetically favorable mode under equi-biaxial loading is square-checkerboard 17,19 , other checkerboard patterns (including the hexagonal and triangular modes) have been observed in actual experiments 17 . A true square-checkerboard pattern can only form under pre-defined conditions using special fabrication methods 43,44 . As a consequence, researchers have speculated the possible causes of the discrepancy. A preexisting curvature of the film surface has been considered as a possible cause to initiate hexagonal-mode wrinkles 17,22 . Nonlinear elastic property of the substrate 17 and unequal substrate elastic moduli in tension and compression have also been proposed as the origin for the hexagon-based checkerboard patterns 21 .
It is evident that, even with a flat film-substrate system with a simple elastic behavior, a unified theme for the evolution of wrinkle patterns is lacking. Furthermore, how transitions occur from one post-instability mode to another is not at all clear. The effects of loading biaxiality are also in need of investigating. In addition to demonstrating the modeling capabilities using the embedded imperfection approach, the present work seeks to address all these issues and offer an overarching scheme of surface wrinkling under biaxial loading.

Numerical model description
The overall problem geometry and boundary conditions are schematically shown in Fig. 1. A thin-film layer, with thickness t f = 0.1 µm , is on top of a compliant substrate. Both the film and substrate materials are taken to be isotropic linear elastic in all simulations. In this study we consider two film materials separately: P3HT:PCBM (poly-3-hexylthiophene conjugated polymer and phenyl-C61-butyric acid methyl ester fullerene derivative) and PEDOT:PSS (poly-3,4-ethylenedioxythiophene and polystyrene sulfonate acid). These are common polymeric thin films used in organic optoelectronic devices. For the P3HT:PCBM film the elastic modulus is E f = 7300 MPa 9 ; and for the PEDOT:PSS film it is E f = 2000 MPa 45 . This large difference in stiffness helps establish the generality of the current numerical predictions. The Poisson's ratio is taken as ν f = 0.35 for both film materials. The substrate material is PDMS (polydimethylsiloxane), with elastic modulus of E s = 2.97 MPa 46 and Poisson's ratio of ν s = 0.495 (set slightly smaller than 0.5 to avoid potential convergence issues).
The simulation domain is defined in such a way to represent a unit cell of a periodic structure. This periodiccell approach is built upon our earlier 1D wrinkling study 37 and is now extended to full 3D simulations. One embedded imperfection is placed in the substrate adjacent to the film-substrate interface, as shown in Fig. 1. The pre-existing defect is a regular finite element carrying the elastic property of the film material instead of that of the substrate. To preserve symmetry, the imperfection is of the square shape and exactly at the center of the xz-plane (note that all the elements in the model are initially perfectly square shaped in the xz-plane). Following the same approach in our two-dimensional simulations 37 , the size of the model is scaled by the analytical critical wavelength of 1D wrinkles ( cr as defined in Eq. 1). The domain dimensions of W x = W z = 10 cr and www.nature.com/scientificreports/ overall depth D = 5 cr were chosen. The film surface is thus of square shape, and the substrate is sufficiently thick compared to film thickness. As a consequence of such scaling, the domain dimensions for the cases with P3HT:PCBM film and PEDOT:PSS film were not identical. Nevertheless, this approach is more advantageous since it provides better control over the total number of elements and leads to computationally efficient simulations with straightforward comparisons. It should be noted that the domain dimensions of W x = W z = 10 cr were intentionally chosen such that 10 cr ∼ = 7( cr ) Cb , meaning that the model size is also sufficiently large for the formation of 7 cycles of square-checkerboard waves (based on Eq. 4) in each of the x and z directions. In addition, the model size is presumably suitable for the post-buckling studies of herringbone and labyrinth patterns, following the rationale discussed in the Supplementary Information that the short wavelength of the surface wrinkling patterns remains invariant and identical to the wavelength of 1D mode in post-instability regimes. Therefore, the model size is appropriate for capturing a variety of wrinkling features, and is also suitable for verification studies by comparing the results with available analytical solutions. It is noted that, since the embedded imperfection is a regular finite element, the imperfection size will be affected as the mesh density is altered. From our previous studies, the simulated surface wrinkling features were found to be essentially independent of the imperfection size 37,38 and material property 36 as long as mesh convergence is achieved. It was also shown that the out-of-plane imperfection thickness of 0.5 t f , used also in the current study, has a negligible effect on the numerical prediction 38 . In addition, it was further illustrated 37,38 that the dependency on imperfection distribution can be avoided if one uses an appropriately sized periodic unit-cell model with only one imperfection at the center as shown in Fig. 1. Therefore, there is a high degree of generality of the present numerical approach.
The simulations were performed under displacement control with a full range of biaxial compression considered. The boundary conditions and applied displacement directions are schematically shown in Fig. 1b, c. The roller boundary condition is imposed on faces z = 0 and x = 0 with only tangential slide allowed. The node at the origin (a corner point at the bottom of the substrate) is entirely fixed in space. The faces z = W z and x = W x are constrained to remain perpendicular to the z and x axes, respectively, during deformation 47 ; and the bottomsubstrate and top-film faces are traction-free. The compressive displacement was applied incrementally with its magnitude varying as 0 ≤ u x ≤ u x (and similarly 0 ≤ u z ≤ u z ), where u x and u z are the maximum applied displacements in x and z directions, respectively. The displacement increments for the static analysis were kept sufficiently small to avoid the potential increment-size dependency of the solutions 37 . The applied nominal strain therefore varies as 0 ≤ e xx ≤ e xx (and 0 ≤ e zz ≤ e zz ) where e xx and e zz are the maximum applied compressive strains in x and z directions, respectively. The biaxiality ratio (BR) is then defined as Due to the square geometry in the xz-plane, the strain ratio can also be written in terms of the applied displacements ( u z and u x ) as seen in Eq. (7). Simulations for various biaxiality ratios were performed in that u x was kept constant and u z was altered within the range of 0 ≤ u z ≤ u x for different cases of BR. As a consequence, the biaxiality varies within the range of 0 ≤ BR ≤ 1, where BR = 0 pertains to uniaxial compression (applied in the x-direction) and BR = 1 corresponds to the perfectly equi-biaxial compression. Any other BR values represent non-equi-biaxial loading conditions.
The simulations were performed using the finite element software package ABAQUS (Version 2017, Dassault Systems Simulia Corp., Johnston, RI, USA). The 20-noded second-order continuum brick elements were used throughout the model. A uniform element distribution was used for the film layer (four layers of elements along the thickness of the film). A graded element distribution is applied for the substrate, with the element size increasing gradually from top (interface) to bottom. The detailed procedure of mesh generation and element definition can be found in the ABAQUS user's manual 48 . To stabilize the numerical solutions in the post-buckling state, viscous damping was included and the corresponding damping factor was calculated via the adaptive automatic stabilization scheme 48 for each time-increment via an iterative process until the converged solution is ensured. The iteration is controlled by the convergence history and the ratio of the energy dissipated by the viscous damping to the total strain energy (termed "accuracy tolerance") 48 . In this study the accuracy tolerance was specified to be less than 10% to prevent damping-dependent solutions. Note that incorporating damping is discretionary, however it facilitates more efficient computations for large-scale simulations. It is worth mentioning that viscous damping was not employed in our earlier studies featuring problems with smaller scopes [35][36][37][38] .
The simulations were conducted using the Message Passing Interface (MPI) parallelization technique, and were implemented at the Center for Advanced Research Computing (CARC) at University of New Mexico. A total of 50 computing nodes (400 Intel-Xenon cores) were employed for each simulation. Preliminary analyses have verified that the numerical solutions were independent of the parallel computational procedures.

Results and discussion
Model verification. Verification of the numerical models in conjunction with the mesh convergence analysis are presented first. Consider uniaxial compression (BR = 0) applied to the film-substrate system. A typical form of simulated 1D wrinkles using the converged mesh is shown in Fig. 2a. The mesh refinement scheme is based on the number of elements in the x and z directions per critical wavelength of the 1D mode ( cr in Eq. 1), and equal mesh size in both x and z is maintained to ensure a square element shape in the plane. This refinement practice provides a neat control over the total number of elements in a wrinkled structure. Figure 2b shows the simulated critical wavelength of the 1D mode as a function of the number of elements per unit wavelength. The theoretical values are also included as horizontal lines for comparison. Three different thin-film materials were studied: a single-layer P3HT:PCBM, a single-layer PEDOT:PSS, and a composite (bi-layer) film comprising of  (6), valid for a single-layer film, need to be revised for a bilayer film via the definition of effective composite film moduli under bending and axial deformations 2,49,50 . For brevity they are not listed here.) As can be seen from Fig. 2b, mesh-insensitive solutions were achieved when the mesh density reaches about eight elements per unit wave. Further increases in mesh density result in constant critical wavelengths matching the theoretical values. Figure 2c shows how the amplitude of the 1D wrinkles (normalized by the film thickness) evolves with the applied strain (normalized by the critical value at the onset of primary bifurcation), for the two cases of singlelayer film material. The mesh density used is the converged 10 elements per wavelength. The corresponding plot for the bi-layer composite film is shown in Fig. 2d. The theoretical response displayed in Fig. 2c, d are based on Eq. (6). It is evident that the numerical results are generally in agreement with the analytical solutions in all cases. From the discussion in Supplementary Information, the critical wavelength of 1D wrinkles is the shortest wavelength compared to all other primary and post-buckling modes under general biaxial loading. In the remainder of this paper, all simulation results were based on the sufficiently fine mesh of 10 elements per wavelength to ensure accuracy. www.nature.com/scientificreports/ As a part of the verification study, we now present the evolution of stresses using the case of P3HT:PCBM film under uniaxial loading (BR = 0). A column of elements along the film-thickness direction from top (free surface) to bottom of film (adjacent to interface) and into the top substrate element, near the lower-right corner of Fig. 1c, were selected. Figure 3a shows the history of the elemental stress component σ xx in all these elements (evaluated at the elements' centroids) as a function of the applied strain in the x-direction, e xx . It is clear that the stress developments in all the film elements considered are identical at the pre-instability stage, and bifurcation starts at the same point with the σ xx value in agreement with Eq. (2). Upon instability the top and bottom film elements selected are at the concave and convex sides of a wrinkle, respectively, so their stresses become more and less negative, respectively, as deformation progresses. Note that the stress in the compliant substrate immediately below the interface remains close to zero throughout the deformation. The same type of stress history is plotted in Fig. 3b, using four random elements at the top surface of the film. The locations of the chosen elements are highlighted in the inset of Fig. 3b. Again, bifurcation initiates at the same point consistent with the analytical critical stress value. The post-instability stress diverges depending on the element location.
In addition to uniaxial loading, the critical stresses associated with the primary bifurcation mode for various biaxiality ratios were obtained from the numerical simulations and compared with the theoretical value of Eq. (2). The variation of the magnitude of (σ xx ) cr with the BR value is shown in Fig. 3c. As can be seen, the simulated (σ xx ) cr is independent of the biaxiality ratio and is close to the theoretical prediction. It is worth mentioning that, for biaxial loading, (σ zz ) cr will change with BR so overall the critical stress state is actually a function of the load biaxiality. While the case of BR = 0.7 was chosen arbitrarily, this loading ratio does lead to a series of pattern changes from 1D mode (primary) to herringbone (secondary) and then to labyrinth (tertiary), as shown in Fig. 4a. The entire sequence from the pre-instability stage can be directly captured in one simulation run. From the model www.nature.com/scientificreports/ output, the reaction forces in either compressive direction may be plotted against the applied strain. An example is shown in Fig. 4b using the reaction force in z, F z , and the applied strain e xx . The curve spans a wide range of deformation history from pre-instability to well into the tertiary instability mode, with the representative topview surface patterns included as insets in the figure. Whenever the instability mode change occurs, the curve displays a sudden change in slope or load drop. Before the onset of instability, the film-substrate system is under nominally uniform compression and the surface remains generally flat. In Fig. 4b, we include a pre-instability image with a highly amplified displacement map, illustrating that the formation of wrinkles is in fact a gradual process starting from very small disturbances around the imperfection site 38 . The onset of primary instability is at the applied strain of approximately 2.6 × 10 -3 , and the dominant 1D wrinkle mode lasts through the strain of approximately 3.6 × 10 -3 when the secondary bifurcation (herringbone mode) commences. As the strain approaches 4.0 × 10 -3 , an abrupt load reduction occurs which signifies a significant change in pattern into the labyrinth mode. Subsequent aberrations of the curve correspond to further adjustments of wrinkle configurations which will change the labyrinth form as the applied strain continues to increase. The simulation approach adopted here enables continuous visualization of the evolving instability patterns and the correlation with the macroscopic mechanical response.
It should be noted that the bifurcation points identified in Fig. 4b depend on the material properties and loading condition. The specific case of P3HT:PCBM film with BR = 0.7 is chosen for illustration. Figure 4c shows the reaction force F z versus applied strain e xx for seven biaxiality ratios ranging from 0 to 1. The corresponding plot for F x is shown in Fig. 4d. The e xx value is used as a reference for both figures to represent the extent of applied strain. In Fig. 4c the case of BR = 0 is for uniaxial loading along the x-direction so F z remains zero throughout the history. From Fig. 4c, d, a linear elastic response is observed at the pre-instability stage in all cases. In Fig. 4d a higher BR value results in a higher initial slope due to the Poisson's ratio effect of biaxial loading. For each BR the first change in the load-strain slope corresponds to the onset of the primary instability mode. Subsequent irregularities along the curves are the result of mode transformation or configuration change, as described in the previous paragraph. (Note that the curve in Fig. 4b is a zoomed-in view of the curve of BR = 0.7 in Fig. 4c.) The curves in Fig. 4c, d indicate that different biaxiality ratios result in significantly different post-instability behavior. The following section is then devoted to a systematic analysis of the wrinkle pattern formation over the entire BR span.

Wrinkle patterns versus biaxiality ratio.
To assess the dependency of wrinkle patterns on the biaxiality ratio, extensive simulations were conducted with the pattern evolution closely monitored. Due to the nature of the resulting patterns as discussed below, the cases of 0.9 < BR ≤ 1.0 require more detailed examinations. Therefore, this range is to be presented separately below. The progression of deformation is characterized by the normalized applied strain in the x-direction, e xx /e cr , where e cr is the numerically obtained critical strain corresponding to the first bifurcation point. Figure 5 shows a collection of top-view (xz-plane) wrinkle patterns displayed by the P3HT:PCBM film, for biaxiality ratios between 0 and 0.9 up through an applied strain of e xx /e cr = 5 . The color contours represent the out-of-plane displacement ( u y ) with the red and blue colors being the highest (peak) and lowest (valley) positions, respectively. Note that the same quantitative color scheme is applied to all the surface patterns presented in this paper.
In Fig. 5 the left-most images are patterns right before instability, at e xx /e cr ≈ 0.98 , showing very small disturbances originating from the imperfection site. (A very high scaling factor of 1000 is assigned to these preinstability images to make them visible.) Shortly after the deformation reaches e xx /e cr = 1.0 , all cases of BR from 0 to 0.9 have displayed well defined 1D wrinkles. For high BR values such as 0.9, the wrinkles are perpendicular to the x-direction since e xx was set to be the higher compressive strain than e zz in the models. For uniaxial loading (BR = 0), the 1D mode persists throughout the entire history. For biaxial loading with a non-zero e zz component, the parallel wrinkles will start to curve so as to transform into the herringbone-like structure as deformation continues. This transformation occurs earlier as the biaxial loading becomes more prominent (increasing BR). In addition, another major post-buckling transformation from herringbone to labyrinth appears for high BR values, which also tends to occur earlier as BR increases.
The herringbone pattern observed in Fig. 5 deserves further discussion. As stated in Supplementary Information, theoretical uncertainties exist in the literature regarding the evolution of herringbone pattern and the necessary loading condition in the post-instability regime. From the current simulations, the pattern is essentially an outcome of lateral undulation (post-bifurcation) in the 1D waves. In addition, this mode was not observed under pure equi-biaxial condition and was captured only for non-equi-biaxial cases of 0 < BR ≤ 0.90. It was observed that the short wavelength (perpendicular to the local wrinkle lines) remains invariant, but the long wavelength (lateral undulation) and the jog angle depend strongly on the applied strain. Our results about the herringbone pattern in general conform to the theoretical considerations in references 19,21,23 . Figure 6 shows the wrinkle patterns at various stages of deformation, for BR values greater than 0.9. With the loading conditions now closer to equi-biaxial, it can be seen that the 1D mode observed previously in Fig. 5 is suppressed right from the onset of first bifurcation. (Immediately before instability a very small disturbance in the form of ripples originating from the imperfection site can be detected under a very high scaling factor.) There is a tendency for the primary instability mode to become checkerboard-like as BR approaches unity. Soon after the primary mode, there is a transient stage where the checkerboard form breaks down and more continuous wrinkles emerge leading to a labyrinth structure. This transformation happens within the strain range of 1 ≤ e xx /e cr ≤ 2 , with the labyrinth mode staying dominant thereafter. Note that the transient patterns seen in Fig. 6 may be interpreted as a hybrid mode of checkerboard-herringbone; the term "bistable modes" have been used in the literature to refer to such hybrid instability patterns 21 www.nature.com/scientificreports/   Fig. 7 for a clear comparison. Within the range of 0 ≤ BR ≤ 0.91 (Fig. 7a) the mode is 1D wrinkles. It is worthy of note that, for a biaxial loading state with the applied compressive strain in one direction (z in the current study) as high as above 90% of that in the other direction (x), the primary instability mode is still 1D wrinkles. Some 1D wrinkles become "branched" when BR is increased to the range of 0.92 ≤ BR ≤ 0.94, which is also associated with the lateral kinks displayed by most of the waves (Fig. 7b). This is a transitional state from the 1D wrinkle to checkerboard modes. Figure 8a, b shows the evolution of branching patterns for BR = 0.92 and 0.94, respectively. In both cases the 1D wrinkles at the two opposite sides of the model appear to emerge in a staggered manner. With further deformation the two sets of waves link up, resulting in a kinked and branched structure. As BR is increased to 0.95 (Fig. 7c), the continuous waves are entirely broken up by the increasingly dominant kinking and branching, and the surface topography tends toward discrete islands.
Beyond BR = 0.95, a periodic wave form in both the x-and z-directions becomes more distinct and the checkerboard pattern takes shape. Note that the true square-checkerboard is obtained only when BR = 1 (Fig. 7f). With a slight deviation from perfect equi-biaxiality, e.g., at BR ≈ 0.99, the wavelength along the z-direction is greater than that along the x-direction (Fig. 7e). It is worth mentioning that this non-square checkerboard pattern is   Fig. 7 can be summarized as follows: • The wrinkling patterns are highly sensitive to the loading biaxiality, and variations of the checkerboard pattern are observed for BR ≥ 0.95. All these simulation results are with an initially perfectly flat surface, influenced only by the load biaxiality. This aspect has not been predicted by any theoretical studies discussed in Supplementary Information. • With BR ≥ 0.95 the loading state is very close to the pure equi-biaxial compression. This may explain why the perfect square-checkerboard pattern has not been observed in earlier experimental studies 1,3,17 (in addition www.nature.com/scientificreports/ to other possible reasons discussed in the literature), since a nominal equi-biaxial state in actual experiments is likely to have slight deviations so non-square types of checkerboard patterns may be easily triggered. Furthermore, any non-uniformity in the film and substrate materials may also affect the local stress state.

Construction of instability phase diagrams.
The results presented in Figs. 5 and 6 provided an overview of how the instability pattern evolves as deformation progresses for the various biaxiality states. Since our simulations are able to predict temporal evolution of one bifurcation mode to another, the critical biaxial strains for each mode may be graphically presented so as to create a phase diagram. Figure 9a shows the variation of critical strains in the x-direction, (e xx ) cr , as a function of biaxiality ratio for our model system of P3HT:PCBM film on PDMS substrate. (Different materials will lead to different quantitative values but the qualitative features  www.nature.com/scientificreports/ remain the same, as confirmed by our analyses using the PEDOT:PSS thin film.) Note that, with any combination of BR and (e xx ) cr values in Fig. 9a, (e zz ) cr is also determined. The critical strains here were obtained from the load-strain curves in conjunction with the wrinkling patterns, as discussed in previous sections. The critical values for the primary, secondary, and tertiary bifurcations are denoted as (e cr ) 1 , (e cr ) 2 , and (e cr ) 3 , respectively. These curves serve as the boundaries between the domains whose characteristic instability patterns are also displayed. Two theoretical values for the first bifurcation are included for comparison, one for the square-checkerboard mode and the other for 1D wrinkles under the plane strain condition. It is evident from Fig. 9a that the critical strains decrease as the biaxiality ratio increases. These domain boundaries converge toward the equi-biaxial state (BR = 1). However, at BR = 1 the three critical values do not actually coincide, as shown in the inset. There is a very small primary instability region where the characteristic pattern is checkerboard. The critical strains, (e cr ) 2 and (e cr ) 3 , on the other hand, are indistinguishable and thus the square-checkerboard transforms to labyrinth as the compressive strain increases. In the vast majority of biaxial loading conditions, 1D wrinkling is the primary bifurcation mode. Herringbone exists between 1D wrinkles and labyrinth. The stable herringbone pattern, however, occurs only within a relatively narrow region in this diagram.
An alternative form of the phase diagram is shown in Fig. 9b, using e xx and e zz to characterize the biaxial state (note that both axes are normalized by (e cr ) 1 along the x-direction). The diagonal line represents the border of this diagram, since e xx ≥ e zz in all cases of this study without the loss of generality. (Flipping the loading directions will lead to exactly the same results but in the other half of the figure). A dashed line representing BR = 0.9 is added in the figure as a reference (not as a "phase boundary"), to aid in the reading of the map. The horizontal axis itself corresponds to the case of BR = 0. The phase boundaries are the three critical strains (e cr ) 1 , (e cr ) 2 and (e cr ) 3 normalized by (e cr ) 1 (corresponding to the respective BR values). As in Fig. 9a, the three main domains featured in Fig. 9b have the characteristic instability patterns of 1D wrinkles, herringbone and labyrinth. The checkerboard pattern, as well as the transitional hybrid mode of checkerboard-herringbone mentioned in the last section, only exist in a very small region when BR > 0.95 between the pre-instability and labyrinth domains. For all other BR values, the development of bifurcation modes follows the order of 1D wrinkles, herringbone and labyrinth as deformation progresses.

Concluding remarks
The embedded imperfection approach is successfully applied to three dimensions, and a comprehensive study on the surface instability patterns induced by biaxial compression is presented. The technique is easily implemented for structures consisting of thin films above a compliant substrate. Instead of treating the surface patterns discretely under specific assumptions, the generation of wrinkling morphologies and their transformations can be directly obtained from the simulations. The results provide insight and mechanistic rationale for uncertainties seen from past theoretical and experimental considerations. The state of biaxiality is found to influence the surface pattern significantly, and each bifurcation mode can be traced back to certain abrupt changes in the overall load-displacement response. The square-checkerboard pattern has proven to be a dominant bifurcation mode only under strict equi-biaxial loading within a very narrow range of strains; physical experiments would easily miss the condition due to any slight deviation. A lower biaxiality ratio (closer to uniaxial loading) expands the dominance of the 1D wrinkles. However, for a biaxiality ratio as high as 0.9, the primary bifurcation mode is still 1D wrinkles except that it transforms into herringbone and then labyrinth patterns quickly. The phase diagrams constructed from the simulation results provide an overview of wrinkling configurations over the entire span of biaxial loading. It is worth emphasizing that all the results presented in the current paper are based on simple isotropic linear elastic material behavior with an initially perfectly flat film layer. The 3D numerical model is well suited for future investigations involving more complex material, geometric and loading conditions.