Surface wrinkling of a film coated to a graded substrate

We study the surface wrinkling of a stiff thin elastic film bonded to a compliant graded elastic substrate subject to compressive stress generated either by compression or growth of the bilayer. Our aim is to clarify the influence of the modulus gradient on the onset and surface pattern in this bilayers. Within the framework of finite elasticity, an exact bifurcation condition is obtained using the Stroh formulation and the surface impedance matrix method. Further analytical progress is made by focusing on the case of short wavelength limit for which the Wentzel-Kramers-Brillouin method can be used to resolve the eigenvalue problem of ordinary differential equations with variable coefficients. An explicit bifurcation condition is obtained from which asymptotic the critical buckling load and the critical wavelength are derived. In particular, we consider two distinct situations depending on the ratio $\beta$ of the shear modulus at the substrate surface to that at infinity. If $\beta$ is of $\mathcal{O}(1)$ or small, the parameters related to modulus gradient all appear in the high order terms and play an insignificant role in the bifurcation. In that case, it is the modulus ratio between the film and substrate surface that governs the onset of surface wrinkling. If, however, $\beta\gg1$, the modulus gradient affects the critical condition through leading-order terms. Through our analysis we unravel the influence of different material and geometric parameters, including the modulus gradient, on the bifurcation threshold and the associated wavelength which can be of importance in many biological and technological settings.


Introduction
The wrinkling of thin layers on substrates is a universal morphological phenomenon of great importance in engineering materials and biological tissues [54,35].A typical set-up to study wrinkling is to consider a thin hard film coated to a soft substrate.If the film shrinks, grows, is prestrained, or the entire system is loaded, the film may develop wrinkles [1].This type of wrinkling morphologies appear in biological systems and is known to establish some of their physiological functions such as folds in human brain [10,38,4] or in gastrointestinal tracts and airway walls [5,30,29].From the viewpoint of mechanics, wrinkled or buckled shapes in skins, fruits or vegetables can be viewed as a result of a bifurcation process related to growth or aging [48,94,25].In engineering, wrinkling in film/substrate bilayers can be used to design specific patterns that can be used to alter the optical property [83,56,52], to regulate surface adhesion or tension [19,93,51], to optimize the wetting characteristic [95,68], to assist in measuring the material properties of soft polymer networks [76,88,17,18], or to help design novel flexible sensors [96,87,50].Moreover, certain surface morphologies brought by structure buckling are expected to occur and will play a positive role in some special cases, such as the creation of gecko's palm surfaces with adhesive force [3].Therefore, it is of great significance to explore the buckling and post-buckling of film/substrate bilayers to understand and further to control multiple pattern formations.
The fundamental importance and potential applications of film/substrate bilayers have led to a flourish of work related to instability problem when, typically, both film and substrate are homogeneous materials.Early works were dedicated to determining the onset of surface instabilities induced by axial compression [28,72,63,7,77,14].Within the framework of nonlinear elasticity, [13] conducted both the linear and weakly non-linear analyses for a compressed neo-Hookean bilayer.The main finding is that there exists a critical value of the modulus ratio between two layers, across which a transition between supercritical and subcritical bifurcations occurs.Later, the influence of compressibility [61,15], curvature [43], pre-stretch [75,41], growth [53,55,9,64], multi-layering [22,84,97], viscosity [40,39], fibers [78], and anisotropy [42,65] were addressed theoretically, numerically and experimentally.In particular, [33] studied the buckling and postbuckling of a film-substrate bilayer when they have similar material properties, and [41] extended the analysis in [13] by considering the pre-stretch effect.Recently, [1,2] revisited this classical problem as well and further explored the growth induced surface instability for which higher-order asymptotic solutions for the critical buckling loads and the amplitude equation were obtained.
In many systems, such as human airway wall [29,62] and hydrogels with graded cross-linking density and surface modified polydimethylsiloxane (PDMS) [66,47,86], the elastic modulus or Poisson's ratio of the substrate may be position dependent due to the spatial variation in the microstructure, leading to a graded substrate.Previous studies have looked at the initiation of surface buckling for an elastic graded half-space or block [49,90,91,92,27,21].In particular [92] adopted an exponential modulus distribution and obtained an explicit bifurcation condition for a compressed block.Yet such choice for modulus functions lead to zero modulus at infinity in a halfspace.For film/substrate bilayers, [16] investigated surface instability of a stiff film on an elastic graded half-space.They derived an analytical solution of the critical compressive strain and the critical wavelength for power-law grading modulus and exponential grading modulus, respectively and the post-buckling solution was studied numerically.[44] established a semi-analytical finite element model to obtain the bifurcation threshold for an elastic graded cylinder covered by a stiff film under axial compression.[20] further developed an analytical model for the buckling of a stiff film on a graded half-space and performed a finite element analysis to validate the analytical model as well as to investigate the post-buckling behavior.
When the material has a graded structure, the governing equations include space-dependent variable coefficients which is particularly difficult to study analytically.For this reason, most studies have resorted to numerical methods to solve the bifurcation problem.The present work is interested in obtain analytical methods for the wrinkling problem in graded material in order to obtain a comprehensive picture on the role of geometrical and material parameters.Here, we use the fact that for thin stiff film resting on a soft graded substrate, the critical wavelength is small compared to the length of the film, indicating a large critical wavenumber [16,20].Hence, we can use the WKB Method (WKB stands for Wentzel-Kramers-Brillouin) [37] to constructing asymptotic solutions for ordinary differential equations with small parameters.It has been used by [31,6] to derive asymptotic results of buckling of spherical shells of arbitrary thickness.After that this asymptotic method was successfully applied to deal with the bifurcation analysis arising from the problems of cylindrical tubes [70,69,59,46], spherical shells [36,45,26], and bending of blocks [24,71].
In this paper, we work within the framework of finite elasticity and study how the modulus gradient affects the critical buckling condition of a bilayer system with a graded substrate.We obtain explicit form for the critical buckling stretch and the corresponding critical wavenumber for the case of a compressed block and the case of a growing bilayer as shown in Figure 1.In Section Figure 1: Schematic illustration of a film coated to a graded half-space under uniaxial compression or growth.The initial length of the bilayer is L and the thickness of the film is h 0 .In both loading cases, the sliding end conditions are used such that there is no shear stress components in the two ends.In addition, the length L is fixed in the growth scenario so the stress is induced by restricted growth.
2, we characterize the primary deformation prior to surface instability and present the incremental theory under in-plane compression and restricted growth.A linear bifurcation analysis is carried out using the Stroh formulation and the surface impedance matrix method in Section 3. We use WKB method to obtain an explicit bifurcation condition in Section 4. The asymptotic solutions are established in Section 5 for β small or of O(1) where β is the ratio of the shear modulus of the surface of the substrate to that at infinity.The reverse case where β is large is explored in Section 6.

Theoretical model
We consider a film/substrate bilayer composed of incompressible hyperelastic materials under compression or growth.The shear modulus of the substrate is assumed to decrease or increase uniformly in the thickness direction (away from the film).We assume that the bilayer is in a state of plane-strain deformation and therefore restrict our attention to the planar problem as shown in Figure 1.The origin O is located at the left end of the film surface, and the X 1 -and X 2 -axes correspond to the length and depth directions, respectively.We denote the initial length of the bilayer by L and the initial thickness of the film by h 0 .The substrate is assumed to be a half-space.Therefore, in the reference configuration the body B 0 occupies a region [0, L] × ([0, h 0 ] ∪ [h 0 , ∞)), and any material point can be represented in the orthonormal basis {e 1 , e 2 } by its initial position X = X 1 e 1 + X 2 e 2 .
We assume that the film surface (X 2 = 0) is traction-free and that there is no displacement as X 2 → ∞.On the vertical sides, the bilayer is allowed to slide along but not penetrate the sides at X = 0 and L. We assume perfect bonding of the interface (X 2 = h 0 ) so that both traction and displacement are continuous across this interface.
When loads are applied through uniaxial either compression or growth the bilayer undergoes a finite deformation and the resulting configuration, denoted by B r , occupies the region [0, l] × ([0, h] ∪ [h, ∞)) for the compression case or within the region [0, L]×([0, h] ∪ [h, ∞)) for the growth case.The position of a point originally at X is now x(X) = x 1 e 1 + x 2 e 2 .The deformation gradient associated with the deformation B 0 → B r is given by As the applied external load reaches a critical value, surface wrinkling may emerge as a bifurcation and this state is referred to as B t .Since there are two layers in the model, we use an upper hat for any quantity belonging to the film and no hat otherwise.For example, Ŵ represents the strainenergy function for the film while the strain-energy function for substrate is written as W . Since the derivations of the governing equations are similar for both layers, we only show explicitly the computation for the substrate.The equations for film are obtained by changing variables.We consider two loading scenarios, namely, uniaxial compression and restricted growth.In both cases, we will first characterize analytically the base state B r and then derive the linearized incremental governing equations from B r → B t .

Uniaxial compression
In the compression case, the deformation gradient (2.1) for the base state B r reduces to where λ 1 and λ 2 are the principal stretches.The incompressibilily condition det F = 1 implies λ ≡ λ 1 = λ −1 2 and the length and height of this deformed body are: Denoting the strain-energy function of the substrate by W = W (F) or, with a slight abuse of notation, in terms of the principal stretches by W (λ 1 , λ 2 ), the Cauchy stress tensor σ in B r is given by where p is the Lagrange multiplier enforcing the incompressibility condition and I is the secondorder identity tensor.Neglecting body forces, the equilibrium equation for a static problem is: where 'div' is the divergence operator in the current configuration.From the traction-free boundary condition at y = 0 and the continuity condition on the interface y = h, we obtain where W 2 = ∂W/∂λ 2 .

Restricted growth
In this scenario, the length between the two rigid plates in Figure 1 is fixed to be L and the deformation is induced by growth.We refer to the theory of volumetric growth [35] that uses a multiplicative decomposition of (2.1): where A is an elastic deformation tensor and G is the growth tensor.We assume that both A and G are diagonal tensors such that where α i (i = 1, 2) represent the elastic principal stretches, and g i are growth factors.There is no volume change in the i-direction when g i = 1, and g i > 1, g i < 1 represent growth and shrinkage, respectively.In this work, we only consider volume increase, meaning that the growth factors are always greater than unity.Since the deformation in the x-direction is limited so that λ 1 = 1, we have α 1 = g −1 1 .Further, according to the elastic incompressibility condition det A = 1, we have α 1 α 2 = 1 and hence α 2 = g 1 and λ 2 = g 1 g 2 .Here we denote the strain-energy function by W = W (A) (or W = W (α 1 , α 2 )), and an explicit expression of the Cauchy stress (2.4) is given by [35]: According to the equilibrium equation (2.5), the traction-free boundary condition on y = 0, and the traction continuity condition on the interface y = h, we find p = ĝ1 Ŵ2 , p = g 1 W 2 , (2.10) where W 2 = ∂W/∂α 2 .

Incremental theory
Next, we are interested in the possible existence of some states B t close to B r .We obtain them by superimposing an infinitesimal displacement field u(x) = u 1 (x 1 , x 2 )e 1 + u 2 (x 1 , x 2 )e 2 on B r .In this way, the position vector of a material particle relative to the common orthonormal basis is denoted by x, which satisfies the following relation x(X) = x(X) + u(x). (2.11) The deformation gradient arising from where repeated indices imply a summation from 1 to 2. The linearized incompressibility constraint then reads where 'tr' is the trace operator.
The first Piola-Kirchhoff stress tensors in B r and B t are denoted by π and π, respectively.It is convenient to introduce the incremental stress tensor χ by where J = det F and 'T' is the transpose.The incremental equations for the substrate in terms of χ are and, similarly, for the film: div χT = 0. (2.16) The incremental boundary condition and continuity condition are: (i) the traction-free boundary condition, χe 2 = 0, on x 2 = 0; (2.17) (ii) the continuity condition, (iii) the decay condition, (iv) the sliding-end condition, χe 1 = 0, u 1,2 = 0, on x 1 = 0 and l. (2.20) For an incremental deformation, the magnitude of η is small and we can therefore expand χ in η and retain the linear terms to obtain where the pressure p has been defined in (2.4), ṗ is the incremental counterpart, and the fourthorder tensor A defines instantaneous moduli whose components are given by , in the compression case, (2.22) , in the growth case. (2.23) The set of equations (2.15)-(2.20)forms a linear system for which we seek solutions.More precisely, we are interested in finding values of the parameter (either compression or growth), such that there exists a non-trivial solution to this system.These values are the critical values for the bifurcation.

Stroh formulation and the surface impedance matrix method
To describe the initial wrinkling pattern, we look for a solution of the form where U (x 2 ), V (x 2 ) and P (x 2 ) are unknown functions to be determined and n is the wavenumber.Note that in the current problem the same n is used for both the substrate and the film.It then follows from the sliding end condition (2.20) that where k is the dimensionless wavenumber and l = λL.Remember that in the growth case the principal stretch λ is always unity.The components of the incremental stress tensor χ can be expressed as where T 12 (x 2 ) and T 22 (x 2 ) are unknown functions.
To simplify book-keeping, we scale all coordinates and variables of length dimension by L. In addition, all variables of the stress unit are normalized by µ h .Here L is the length of the bilayer and µ h is the shear modulus on the substrate surface (X 2 = h 0 ) yielding where μ is the uniform ground state shear modulus for the film and γ measures the modulus ratio between the film and the substrate surface.To further simplify notations, we drop the asterisk from now on.We obtain relationship between U (y) and V (y) from the linearized incompressibility condition (2.13): where a prime denotes a derivative with respect to y.Similarly, P (y) can be solved from the expression of χ 22 : Next, we define the displacement vector u(y) = [U (y), V (y)] T , the traction vector T(y) = [T 12 (y), T 22 (y)] T , and write their combination as a displacement-traction vector Based on the equations (2.15), (2.21), (3.5) and (3.6) we recast the boundary-value problem given by (2.15)-(2.20)as a differential equation for the displacement-traction vector: This way of posing the problem is known as the Stroh f ormulation of the implied incremental problem [79,32], where Q(y) is the Stroh matrix possessing the following block formulation: and Q i (i = 1, 2, 3, 4) are 2×2 sub-blocks such that Q 2 and Q 3 are real symmetric, and In particular, the matrices Q 1 , Q 2 and Q 3 are given by Note that in obtaining the expression of Q 3 we have used the fact that A 2112 = A 1221 .We apply the impedance matrix method [11,12,73,74] to solve the incremental problem.To this end, we first assume a relation between the displacement vector u(y) and traction vector T(y): with K(y) an unknown matrix.It can be derived from the decay displacement boundary condition u(∞) = 0 that K(∞) = 0. Now substituting equations (3.7) and (3.11) into (3.8) and eliminating the dependence of T(y), we obtain a Riccati equation for K: For the film with a free lateral boundary, instead of equation (3.11), we use as an ansatz a different relation between û(y) and T(y), which is written as where Z(y) is known as the surface impedance matrix.It follows from the traction-free boundary condition T(0) = 0 that Z(0) = 0.
For the substrate, we obtain the following differential matrix Riccati equation by applying a similar derivation procedure: The expressions of Qi (i = 1, 2, 3, 4) are determined from (3.10) by suitable variable substitutions.We can solve numerically the equation (3.12) subject to the boundary condition K(∞) = 0.In this way, the traction on the interface is given by T(h) = K −1 (h)u(h).Meanwhile, a similar solution procedure for equation (3.14) can be applied and as a result the traction on the interface is written as Z(h)û(h).The continuity condition gives (Z(h which is the bifurcation condition.This bifurcation condition (3.15) is valid for both problems considered in this study.Once the strain-energy function and the shear modulus distribution function are specified, one can obtain out the onset of surface wrinkling by solving (3.15) numerically.

Preliminary for asymptotic analysis
In this subsection, we return to the original incremental equations and transform it into an equivalent fourth-order ordinary differential equation suitable for asymptotic analysis.For that purpose, we substitute (2.21), (3.1) into the governing system (2.13), (2.15) and eliminate all unknowns except for V (y).Ultimately, we acquire a dimensionless fourth-order ODE for V (y) by utilizing the dimensionless variables in (3.4), which takes the following form Moreover, the traction components are given by and the boundary conditions and continuity conditions (2.17)-(2.19)are transformed into T12 = 0, T22 = 0, on y = 0, This formulation is valid for both the compression and the growth cases.

WKB reduction
We consider the short wavelength approximation, i.e. k → ∞, to carry out an asymptotic analysis on the critical condition of surface wrinkling, following [46].
First, to make progress we must choose a material model.Here, we will use the neo-Hookean model but other choices can be used within the same setting.Note that an asymptotic formula for the critical buckling stretch has been derived by [13] when the substrate is homogeneous and the incompressible neo-Hookean model is employed.This will be used as a benchmark for validating the asymptotic solutions in our particular situation.In the plane-strain case, the strain-energy functions for the film and substrate are given by where μ is the homogeneous ground state shear modulus of the film, µ s (X 2 ) is the modulus function, and I 1 = tr FF T if the external load is given by axial compression or I 1 = tr AA T in the growth case.Before proceeding further, we introduce the dimensionless modulus where µ h has been defined in (3.4).

Uniaxial compression
For the neo-Hookean model (4.1), the fourth-order ODE (3.16) simplifies to 3) The boundary conditions and continuity conditions in (3.19) become The above equation (4.3) is for the substrate; its counterpart for the film can be found by replacing V by V and by setting all derivative terms of µ to zero as the film is composed of a homogeneous material.Accordingly, the fourth-order ODE for the film is given by Solving directly the associated characteristic equation yields four eigenvalues s 1,2 = ∓kπ/λ and s 3,4 = ∓kπλ.Correspondingly, the general solution to equation (4.5) can be written as where C i (i = 1, 2, 3, 4) are constants.Following the WKB method, we seek for a solution of (4.3) by using the following ansatz [37] V Substituting the solution (4.7) into (4.3)gives The next step is to look for an asymptotic solution to the above equation of the form where S m (m = 0, 1, 2 . . . ) are unknown functions of y.Substituting (4.9) into (4.8) and equating all coefficients of like powers of 1/k to zero provide infinitely many equations.In particular, the leading-order equation is a quartic algebraic one of S 0 from which four independent solutions can be obtained.Typically, a two-term approximation containing S 0 and S 1 in (4.9) already provides sufficient accuracy [46].The four independent solutions for S(y) are: Then, the general solution for V (y) can be written as where C i+4 (i = 1, 2, 3, 4) are arbitrary constants.Inserting (4.6) and (4.11) into (4.4),we obtain a vector equation where T and the coefficient matrix M is given by Note that the parameters âi , bi , a i , b i , q i , and f i are defined as The existence of a non-trivial solution to (4.12) requests the vanishing of the determinant of M. Finally, the bifurcation condition reads det M = 0, (4.15) which can be regarded as a function of the principal stretch λ, the wavenumber k, the modulus ratio γ, the thickness h 0 and the parameters involved in the modulus function µ(y).After a careful observation of the coefficients in (4.14), we find that f 2 and f 4 are both exponentially large so the terms whose factor is f 2 f 4 will dominate the value of determinant.As a result, the bifurcation condition (4.15) can be rewritten as where E.S.T denotes exponential small terms which will be neglected in further analysis and the algebraic equation Ξ = 0 provides a way to obtain an analytical formula for the critical bifurcation threshold.

Restricted growth
When the compressive stress is induced by growth, the bifurcation parameter changes is the growth factor.The strain-energy functions for the film and the substrate are provided by (4.1).In this case, the dimensionless fourth-order ODE (3.16) is where g 1 is the growth factor for the substrate in the horizontal direction.Correspondingly, the incremental boundary conditions and continuity conditions in (3.19) can be written as In our previous study [57], it was found that a gradient in the growth function has a negligible effect on the critical buckling threshold and a constant growth can provide a good approximation on interpreting complicated morphological formations in bilayer systems.Thus, we assume that both layers share the same constant growth factor such that ĝ1 = ĝ2 = g 1 = g 2 = g.Next, we aim to solve asymptotically the eigenvalue problem with variable coefficients arising from (4.17) and (4.18).The fourth-order ODE for the homogeneous film can be obtained by removing all derivative terms of the shear modulus in (4.17).Solving directly the associated characteristic equation gives a general solution of V (y) of the form where C i (i = 1, 2, 3, 4) are constants and the eigenvalues r i (i = 1, 2, 3, 4) are given by r 1,2 = ∓kπ and r 3,4 = ∓kπ/g 2 .
To solve (4.17), we use again a WKB-type method with ansatz given by and the integrand function can be expanded in k as where R m (m = 0, 1, 2 . . . ) are unknown functions of y.On substituting equations (4.20) and (4.21) into (4.17)we also figure out four solutions of R 0 from the leading-order equation.Similar to the solution procedure illustrated in Section 4.1, we are able to obtain the following four independent solutions The general form of V (y) is where C i+4 i = 1, 2, 3, 4 are constants.We further substitute (4.19) and (4.23) into (4.18) to obtain a vector equation of the same form as (4.12) but with all S (i) and s i in (4.13) replaced by R (i) and r i , and entries given by As in the previous case, the WKB analysis generates an approximate bifurcation condition Ξ = 0 from det M = 0, which will be used later on.

The shear modulus various exponentially
The problem cannot be solved for an arbitrary modulus function.Hence, we consider a modulus with an exponential decay of the form where β = µ ∞ /µ h ⩽ 1 and α = ζL is the normalized decay rate.We also impose µ 1 ̸ = 0 to avoid zero modulus [20].Using (4.26), we are able to express the bifurcation condition (4.16) as Ξ(λ, k, α, β, γ, h 0 ).For the growth scenario, the exponential modulus function (4.26) is used as well but we substitute λ = 1/g 2 to obtain a condition in terms of of g of the form Ξ = Ξ(g, k, α, β, γ, h 0 ).A particularly interesting situation arises when the substrate is more compliant that the thin film with γ large and h 0 is small considered in [16] and [20].Our study unifies and generalises both work by considering the two cases β ∼ O(1) and large β.

Asymptotic solution when β is small or of O(1)
We focus on β ≪ 1 or β ∼ O(1) and derive some asymptotic estimates for the critical buckling load as well as for the critical wavenumber.The compression and growth cases will be explored separately.Before carrying out the asymptotic analysis, we first present five bifurcation curves according to the exact bifurcation condition (3.15) in Figure 2 when the modulus ratio γ = 100 and the initial film thickness h 0 = 0.015.The purpose is to identify some valid parameter domains where the short wavelength assumption can be satisfied.In Figure 2 we see that each curve has a maximum, highlighted by a blue point; its vertical coordinate gives the critical stretch λ cr triggering surface wrinkling while the horizontal coordinate gives the critical wavenumber k cr .When α = 5 or α = 10 the critical wavenumber k cr is small, which does not satisfy the short wavelength assumption in the WKB expansion (4.7).Hence, we consider the case where α ∼ O(1).We also note that the specific value of β when β ∼ O(1) has a weak effect on the critical state when the modulus ratio γ is large.

Uniaxial compression
Next, we explore the influence of the film initial thickness h 0 on the critical state.Given γ = 100, α = 1, β = 0.01, Figure 3 shows nine bifurcation curves stemming from (3.15) for different values of h 0 .The dashed line is the envelope of the maximums of bifurcation curves.We observe that the critical wavenumber k cr is a decreasing function of h 0 .To satisfy the short wavelength assumption, we confine our attention to small h 0 as found in many engineering devices or biological tissues.It is worth noting that λ cr is independent of h 0 when the substrate is homogeneous [13, 61, 1, 85], but in our problem λ cr is no longer a constant as shown in Figure 3.In addition, we assume that the film is much stiffer than the substrate, which implies that γ is large.It is also known that ε = 1 − λ cr is small as the film is much stiffer than the substrate.In conclusion, we will place ourselves in the situation where the following four parameters are small: 1/k, h 0 , 1/γ and ε.In Figure 4, we compare the bifurcation curves based on the approximate bifurcation condition (4.16) (dashed line) and the exact one (3.15)(solid line) if 10 ⩽ γ ⩽ 500 and h 0 = 0.015, α = 1, β = 0.01.It can be seen that the asymptotic result is in good agreement with the exact (numerical) one for all cases, even when k cr = 3.Similarly, we checked that the asymptotic bifurcation condition (4.16) remains valid when α ∼ O(1), 1/k, h 0 , ε, and 1/γ are all small parameters.Therefore, we can now start from this approximation (4.16) to obtain explicit solution for λ cr and k cr .) First, we expand equation (4.16) in terms of the small parameters ε and h 0 to obtain where t i (i = 0, 1, 2, . . . ) are constant coefficients related to α and β, whose lengthy expressions are not written out for brevity.The leading-order balance and principle of least degeneracy [61] are used to derive the order relations among ε and kh 0 and γ.Keeping all possible leading-order terms, (5.1) can be written as where t 0 ∼ O(1) and I.H.O represents infinitesimal of higher order.Analyzing and discussing all possible balances following [46,61], we find that when ε ∼ O(γ − 2 3 ) and kh 0 ∼ O(γ − 1 3 ), the leading terms in (5.1) are of O(γ  ) and given by t 6 kh 0 γ, t 10 k 4 h 4 0 γ 2 and t 22 k 2 h 2 0 εγ 2 .Then, we extract all possible second-order terms for equation (5.1) and write where the terms t 8 k 2 h 2 0 γ, t 11 k 5 h 5 0 γ 2 and t 23 k 3 h 3 0 εγ 2 are all of order O(γ 3 ).Using the assumption that γ is a large parameter, we find that O(γ 1 3 ) ≫ t 0 and the second-order must be equal to O(γ 1 3 ).According to the least degeneracy principle, in order to keep the t 7 h 0 γ term, we can determine 3 ).Through the above analysis, the order relations are clearly given by ε 3 ).To derive an asymptotic solution for λ cr and k cr , in addition to the equation (5.1) we also need to compute ∂λ/∂k = 0, or equivalently Making use of the order relations we use where a 0 is a constant of O(1).We then seek an asymptotic solution for λ cr and k cr in the form of where m i and l i (i = 0, 1, 2 . . . ) are coefficients to be determined.
Substituting (5.6) into (4.16) and ( 5.4) provides the four-term asymptotic solutions ) As a check to our method, we recover the homogeneous case by taking α = 0 or β = 1 to obtain where an upper circle stands for the solution corresponding to a homogeneous substrate.Using (3.2) the estimates (5.10) gives + O(γ −5/3 ), (5.11) which is an exact match to the results found in [13,1,85] up to O(γ −5/3 ).A remarkable property of the solution (5.7) and (5.8), is that all parameters relevant to gradients (given by α and β), only appear in higher order terms.In other words, it is the modulus ratio between the film and the substrate surface that dominates the critical bifurcation condition, and we find that the modulus gradient has a minor influence on λ cr and k cr .To further unravel how the modulus gradient affects the critical state, we define ), (5.12) The sign of δ 1 implies whether a bilayer with graded substrate is more stable or more unstable than the corresponding homogeneous substrate.Thus, we summarize the results of the signs of δ 1 and δ 2 as follows Note that the prerequisite that (5.12) and (5.13) hold is β ∼ O(1).Therefore we may allow β > 1 but a large β is not granted.It can be seen that either α < 0, β < 1 or α > 0, β > 1 indicates a positive δ 1 but a negative δ 2 , suggesting that an increasing modulus distribution will delay surface wrinkling and lead to a larger wavenumber compared to its homogeneous counterpart.Contrariwise, the case α > 0, β < 1, with an exponentially decayed modulus, always gives rise to a more unstable structure and a lower wavenumber.These conclusions further confirm previous results [16].We now check the validity our asymptotic solutions given by (5.7) and (5.8). Figure 5 show how λ cr and k cr depend on γ.We see that the error between asymptotic and exact results is negligible for large γ.Further, we observe that the asymptotic solution remains valid even for small γ and at γ ≈ 2 the largest error is less than 0.6%.We also see in Figure 5(b) that (5.8) looses its validity γ ⪅ 5.These solutions also demonstrate that a larger mismatch of the stiffness between the film and substrate always tends to destabilize the structure with a higher λ cr but a lower k cr which is consistent with previous results.
Next we investigate the effect of h 0 in Figure 6.Again the agreement between exact and asymptotic solutions is excellent, and we observe that unless the substrate is homogeneous [61], the critical stretch is affected by h 0 .Indeed, λ cr in (5.9) is depends on h 0 to second and higher order terms for a graded substrate but with factors α and β − 1 that vanish identically for homogeneous substrates.The effect of h 0 is weak as seen in Figure 6(a).This results seems unexpected at first sight as the bilayer is more unstable with a thicker film.However, since k cr is proportional to 1/h 0 , a thicker film develops less wrinkles (see Figure 6(b)).We emphasize that the order relation γ ∼ O(h −3 2 0 ) implies that the film is much stiffer that the substrate for extremely thin films.

Moderately stiff film on soft graded substrate: γ ∼ O(h −1
0 ) From the leading-order balance analysis of (5.1) we discover an implied relation between γ and kh 0 , i.e. kh 0 ∼ O(γ − 1 3 ).We now assume that the film is not as stiff and have, instead, γ = O(h −1 0 ).In this case, we have k = O(γ   are now where m i and l i (i = 0, 1, 2, . . . ) are constants depending on h 0 , α and β and can be solved from equations (4.16) and (5.4).Following the same procedure as in the previous case, we obtain the asymptotic solutions: ) These solutions reduce to (5.9) and (5.10), when either α = 0 or β = 1.Further, the leading order term of (5.16) and the first two terms of (5.17) are identical to those in (5.7) and (5.8).
Nevertheless, an interesting difference is that the lowest order where the information of modulus gradient occurs is O(γ − 4 3 ) in (5.16) while the corresponding order is O(γ −1 ) in (5.7).Hence, the influence of modulus gradient on the critical stretch λ cr becomes negligible as the modulus mismatch between the film and substrate drops.For the critical wavenumber, the first two terms of (5.17) are consistent with those of (5.8), indicating that the buckled morphology is nearly impervious as order between γ and h 0 varies from γ ∼ O(h 5.1.3.Slightly stiffer film on soft graded substrate: γ ∼ O(h 3 ), and we obtain These results are independent of the modulus gradient parameters α, β and the thickness h 0 .Specifically, equations (5.18) and (5.19) coincide exactly with (5.9) and (5.10), respectively, up to the truncated order.This implies that for a stiff film, the graded substrate makes little difference compared to its homogeneous counterpart.
Combining the results for γ ∼ O(h ) and γ ∼ O(h −1 0 ), we conclude that for a film/substrate bilayer, where the film thickness is fixed and the shear modulus of the substrate either decays exponentially from the surface or is comparable to the modulus at infinity, the influence of the modulus gradient tends to shift to higher orders as the modulus ratio γ decreases, playing a more and more marginal role in regulating surface wrinkling as well as in shaping the final surface pattern.

Restricted growth
For a stiff layer coated to a soft substrate, it is known that the critical growth factor g cr leading to wrinkling is close to unity [1].Therefore, ε = g cr − 1 is small.Following the fundamental assumption in Section 5.1, we assume that α is of O(1) and ε, 1/k, 1/γ, h 0 are small.Figure 7: Comparisons of the bifurcation curves between the approximate (dashed line) and exact (solid line) bifurcation conditions when α = 1 and β = 0.01.The left subfigure displays the dependence of g on k by specifying h 0 = 0.015 while the right one exhibits the counterpart for fixed γ = 100.The blue dots correspond to the critical growth factor g cr associated with the critical wavenumber k cr .
To evaluate the accuracy of the approximate bifurcation condition, we compare in Figure 7 the the bifurcation curves the approximate (dashed line) and exact (solid line) bifurcation conditions.The dashed curves are an excellent fit for the solids curves, which indicates the validity of the explicit bifurcation as our new starting point to pursue an asymptotic solution.In Figure 7 we see that the bifurcation curves have a U −shape.The critical growth factor g cr and the critical wavenumber k cr correspond to the vertical and horizontal coordinates of the local minimum for each curve.

Very stiff film on soft graded substrate: γ ∼ O(h
We first expand the bifurcation condition Ξ = 0 in ε and h 0 .Referring to the order analysis in Section 5.1.1 and by using the leading-order balance and principle of least degeneracy [61], we arrive at the same order relations as before, namely, ε ∼ O(γ ).
Therefore, we look for an asymptotic solution of the form where m i and l i (i = 0, 1, 2, • • • ) are constants to be determined.Following the same procedure given in Section 5.1.1,we derive the following formulas ) For a homogeneous substrate, equations (5.21) and (5.22) simplify to + O(γ −5/3 ). (5.24) We emphasize that the leading order terms of (5.23) and (5.24) 2 are the same as those in [1], where there is no growth in the substrate and only film growth is considered.This reinforces our conclusion in [57] that growth gradient has a little effect on the critical buckling load and the associated buckling pattern.Further, it can be seen that the thickness h 0 and the parameters related to modulus gradient α, β are included at higher order terms.So the dominant factor that affects the critical buckling state is still γ.The same conclusion has already been observed in the compression case.Figure 8 illustrates the critical growth factor g cr and the critical wavenumber k cr as functions of the modulus ratio γ when α = 1, β = 0.01, and h 0 = 0.015.The dashed lines correspond to the asymptotic solution while the solid lines represent the exact solution.We note the surprisingly good agreement between the exact and asymptotic solutions, even when γ ≈ 5.It implies that the general order relation ε ∼ O(γ − 3 2 ) captures well yje surface instability in hyperelastic bilayers.We plot the relations of g cr , k cr and h 0 in Figure 9 to further validate our asymptotic solutions.We found from the second order term in (5.21) that the film thickness h 0 starts to affect the critical growth factor g cr , as confirmed in Figure 9(a).Note that h 0 does not impact the onset of surface wrinkling for a homogeneous substrate through (5.22), it will change the wavenumber, as shown in 9(b).We conclude that the asymptotic solutions are excellent approximations to describe the onset of surface instability and the corresponding pattern.
Hence, the effect of a modulus gradient becomes weaker as the modulus ratio γ declines.
5.2.3.Slightly stiff film on soft graded substrate: γ ∼ O(h 3 ) corresponds to a slightly stiff film coated to a soft graded substrate, for which, we have ) which does not depend on α and β.Meanwhile, equations (5.27) and (5.28) are identical to (5.23) and (5.24) 1 up to the retained order, respectively, implying that a homogeneous substrate can indeed replace a graded one when estimating the critical buckling threshold and the associated buckling pattern.In addition, we give the performance of the asymptotic solutions (5.21), (5.25) and (5.27) when γ is slightly greater than unity in Table 1.Again, the relative errors are small even when γ ≈ 2. The solution of γ = O(h 0 ) yields the best prediction.We conclude that in the confined growth scenario when both α and β are of O(1) and both h 0 and 1/γ are small, the influence of the modulus gradient on pattern formation in a stiff film boned to a compliant graded substrate is negligible.For smaller γ, the modulus gradient becomes irrelevant.
Table 1: Comparisons between the exact and asymptotic solutions of g cr as γ is not large.All other parameters are specified by h 0 = 0.015, α = 1 and β = 0.01.We denote ∆ the relative error between the exact and asymptotic solutions.
γ Exact g cr γ = O(h When the modulus ratio β is large we use the WKB approximation (4.9).Following the method given in Section 4.1, we obtain Both S 1 and S ′ 1 are negative under the assumption β ≫ 1 and α > 0 which implies that S 1 attains its local maximum at y = h 0 /λ with extreme value given by 3) It can be estimated from (6.1) and ( 6.3) that a boundary layer exists after which S 1 decreases.However, near the surface y = h 0 /λ, it is of O(β).The first term in (4.10) (kS 0 ) is of O(γ 3 ) and the second term (S 1 ) is of O(β).Bearing in mind that a valid asymptotic expansion requires that the latter terms are strictly smaller than the previous one, we consider a limit case β = O(γ 1 3 ), which means that the first and second terms are of the same order O(γ 3 ) in the boundary layer, and the first term is much larger outside the boundary layer.Hence, we write where b 0 is a constant of O(1).We seek an asymptotic solution for λ cr and k cr the same as (5.6).By applying the regular solution procedure, it is found that the coefficients l 0 and m 0 satisfy the following algebraic equations 2π 6 a 3 0 l 6 0 − 24π 4 a 0 m 0 l 4 0 + 12π 3 l 3 0 + 18π 2 αb 0 l 2 0 + 3πα 2 b 2 0 l 0 − 3α 3 b 3 0 = 0, (6.5) 16π 6 a 3 0 l 6 0 − 144π 4 a 0 m 0 l 4 0 + 60π where a 0 has been defined in (5.5).Unfortunately, equations (6.5) and (6.6) give rise to a polynomial equation of degree six in l 0 , which has no general solution.To overcome this difficulty, we resort to curve fitting using a specific function and the regression analysis yields l 0 = 0.170341 ln αβ γ 1/3 + 0.552186 + 0.564331 ) We see that the critical wavenumber is affected by the modulus ratio γ only through a logarithm term.A similar dependence was also observed in growing layers [23,46].The leading-order solutions of the critical bifurcation thresholds read λ cr = 1 − m 0 /γ 2/3 and g cr = 1 + m 0 /γ 2/3 while the critical wavenumber is k cr = l 0 γ −1/3 .We then illustrate the critical stretch and the critical wavenumber as functions of β in Figure 10 when γ = 500, h 0 =0.01, α = 1.For comparison, both the leading-order asymptotic solution (dashed line) and the exact one (solid line) are plotted.It is emphasized that an obvious gap occurs as the vertical axial ranges from 0.983 to 0.992.In fact, the maximum error for the critical strain 1 − λ cr is around 5% in Figure 10(a).On the other hand, we find that the leading-order coefficients are only valid when β is less than 40, which is the termination of β in Figure 10.Furthermore, the critical stretch λ cr is a decreasing function of β, indicating that a greater β tends to delay the instability.From Figure 10(b) we see that the critical wavenumber is an increasing function of β.So a larger β will produce more wrinkles, which is consistent with the results in [16].

Discussions and conclusions
We investigated the surface wrinkling of a stiff film resting on a soft graded substrate engendered by in-plane compression or restricted growth.A theoretical model was used to identify the onset of surface wrinkling within the framework of nonlinear elasticity and by use of the incremental theory without specifying a material model and a modulus function of substrate.An exact bifurcation condition was presented by means of the Stroh formulation and the impedance matrix method [58,85].In the short wavelength limit, the WKB technique was applied to deal with the eigenvalue problem of ODEs with variable coefficients, and as a result an explicit bifurcation condition was obtained.By considering an exponential modulus function, the effectiveness of the approximate bifurcation condition was validated based on the exact one, and the valid domains where the WKB expansion is effective were identified for various parameters.
We determined the order relations and showed that the modulus ratio γ is extremely significant in regulating the onset of surface wrinkling in bilayer systems as found before [57,58] and that the critical relation ε ∼ O(γ − 2 3 ) is consistent with most existing studies on pattern formation in film/substrate structures [13,61,15,46,45,85].However, if surface wrinkling is replaced by a Euler-type buckling with mode number 1 or 2, the aspect ratio of the whole structure becomes dominant [60].
When β is not large, we found that the asymptotic solutions for the critical loads and the critical wavenumber give an excellent approximation and that these solutions do not strongly depend on the modulus gradient as it only appears in higher-order terms.However, as gamma decreases, and unlike an homogeneous bilayer, the film thickness h 0 starts to alter the critical load and unexpectedly a thicker film tends to destabilize the structure.Furthermore, surface wrinkling is more likely if the shear modulus of the substrate decays away from the surface.
To compare the axial compression case to the confined growth case, we define the critical thickness h cr at the bifurcation point: The small difference at O(γ −2/3 ) is due to the fact that the structure is allowed to grow in both the horizontal and the vertical directions.If we only consider horizontal growth, the corresponding formulas would be identical up to O(γ −4/3 ).
In the case of large β and assuming β ∼ O(γ 1/3 ), we found the leading-order solution from a regression analysis.In this case, the parameters related to modulus gradient emerge in the leadingorder terms, showing that the modulus gradient can not be neglected when the shear modulus of the substrate increases from its surface.To study the role of grading parameter α that controls how fast the shear modulus of the halfspace decays or grows from its surface, we show in Figure 11(a) the relation between λ cr and α and in Figure 11(b) the relation between k cr and α.We consider two scenarios, namely, exponentially decaying modulus (β = 0.01) and the reverse one (β = 100).When β = 0.01, we observe from 11 that λ cr gradually ascends but k cr experiences a sudden drop from 5 to 1 as α ≈ 5, and the λ cr curve is non-smooth as well.With increasing α, the shear modulus of substrate will decay faster and the system will behave like triple layer structure.If a compliant substrate is coated with two stiff layers a mode transition has been observed as the intermediate layer becomes softening [97].This may explain the jump behavior at α ≈ 5. Indeed, a greater α is equivalent to reducing either the thickness or the modulus of the middle layer.In doing so, the tendency of k cr -α curve in Figure 11(b) (for β = 0.01) is consistent with the conclusion in [22].Moreover, the grading parameter only generates a minor correction to λ cr but creates a large change in the critical mode.When α > 5, the critical wavenumber remains one which indicates that α does not affect the pattern.
When β = 100, with an exponentially growing modulus in the substrate, the results are different.Then both λ cr and k cr are highly dependent on α.In particular, the critical stretch λ cr is a monotonically decreasing function of α while the critical wavenumber k cr is a monotonically increasing function of α.As α goes to infinity λ cr will approach a limit value similar to the one obtained for instability of a soft wire embedded into a half-space.Indeed, for α = 5000, we find λ cr = 0.5724, similar to the Biot value 0.543 at which surface instability in a compressed half-space takes place [8].
Finally, we summarize the effect of the modulus gradient.If the grading parameter α ∼ O(1), our analytical solutions show a relatively small role of the grading for an exponentially decayed modulus in the substrate as the modulus gradient only affects high-order terms.If, however, the modulus increases away from the interface, then the gradient plays a leading role through a firstorder logarithm dependence.For α large, corresponding to a change of modulus in a very thin region, the grading effect has little influence on the critical state but has a very pronounced effect in shaping the pattern for an exponentially decaying modulus.Conversely, for a growing modulus, the modulus gradient has a strong effect on both the critical stretch and the associated pattern.As α → ∞ the bilayer degenerates to a half-space with a soft line inclusion.Then, the classical Biot instability is recovered but with a finite wavenumber.The results presented here give a systematic way to consider wrinkling instabilities in bilayered grading material and a natural complement to computational studies.

Figure 2 :
Figure 2: Bifurcation curves based on the exact bifurcation condition (3.15) for γ = 100, h 0 = 0.015.The blue dots determine the critical principal stretch λ cr associated with the critical wavenumber k cr .
Dependence of k cr on γ.
Dependence of k cr on h 0 .
Dependence of g cr on γ.
Dependence of k cr on γ.
Dependence of k cr on h 0 .
Dependence of k cr on β.
Dependence of k cr on α.

Figure 11 :
Figure11: Influence of the grading parameter α on the critical stretch and the wavenumber with γ = 100 and h 0 = 0.015.We are concerned with the behavior as α ⩾ 0.
Dependence of λ cr on h 0 .