Revisiting the wrinkling of elastic bilayers I: linear analysis

Wrinkling is a universal instability occurring in a wide variety of engineering and biological materials. It has been studied extensively for many different systems but a full description is still lacking. Here, we provide a systematic analysis of the wrinkling of a thin hyperelastic film over a substrate in plane strain using stream functions. For comparison, we assume that wrinkling is generated either by the isotropic growth of the film or by the lateral compression of the entire system. We perform an exhaustive linear analysis of the wrinkling problem for all stiffness ratios and under a variety of additional boundary and material effects. Namely, we consider the effect of added pressure, surface tension, an upper substrate and fibres. We obtain analytical estimates of the instability in the two asymptotic regimes of long and short wavelengths. This article is part of the theme issue ‘Rivlin's legacy in continuum mechanics and applied mathematics’.


Introduction
In the development of rational mechanics after WWII, there was much emphasis on foundations, abstraction and formalism. Yet, by its very nature and applicability to the world around us, mechanics is not an abstract concept. In this era, Ronald S Rivlin stood out as a pragmatic researcher, as he was interested in identifying new phenomena, in conducting experiments and in developing general methods to tackle scientific problems. He was at heart a problem-solver and his enduring legacy is to have given us the tools to understand the world around us. Two of his contributions are particularly noteworthy in this regard.  First, he systematically derived exact solutions for simple geometries (cuboids, spheres, cylinders [1]), and second, together with Green, Shield and Pipkin [2,3], he developed a general method, the small-on-large theory, to obtain both vibration modes and stability regimes for mechanical systems in large deformations.
The present contribution follows in Rivlin's footsteps as it uses both ideas to their full potential by presenting a general theory for the stability of a bilayer system in large deformations.

(a) Multi-layered elastic materials
Multi-layered elastic materials have been widely studied in engineering literature with applications ranging from construction materials [4] to three-dimensional printing. The past few decades have seen a resurgence of interest in the development of a mathematical description of these structures, moving from the classical linear theory of elasticity to a more accurate nonlinear hyperelastic description that is valid for the large deformations relevant in many biological contexts. Additionally, where traditional analyses often made use of thin-film approximations [5] or extreme ratios of stiffness between layers, newer works have begun to explore the intermediate parameter space.
Owing to the complex nonlinearities present in this problem, it is often too difficult to compute the solution of boundary-value problems explicitly in this setting. As such, we must often make use of indirect methods such as perturbation expansions to identify bifurcation points in the system as parameter values are changed, as well as the stability of any non-trivial solutions that may emerge. Such methods have been used to great success to compute the critical uniaxial compression required to cause buckling of an elastic half-space coated in a thin, stiffer elastic film [6]. Following works have considered variations of the physical setting such as pre-stretching the substrate [7], further compressing the buckled bilayer to induce a second, periodic-doubling bifurcation [8], the limiting behaviour of the system as the stiffness ratio of the layers tends to unity [9], the effect of adding reinforced fibres to the substrate [10] and the replacement of compression with growth as a mechanism to induce buckling [11].
This final modification is of particular relevance in the study of biological materials, as it allows us to model the formation of complex organs such as the brain in utero [12]. The physical structure of mammalian brains consists of distinct layers of cells with similar, but different mechanical properties and thicknesses. In particular, we can divide the brain into the outer layer of grey matter (the cortex), which primarily consists of neuron cell bodies, and the inner white matter (the subcortex), which primarily consists of axons and their insulating myelin sheaths. The process that leads to the evolution of the characteristic convolutions (gyri) of the brain has only recently seen convincing explanations after several decades of scientific research [13]. There are several competing mechanical theories [14] that seek to describe the evolution of cortical folds, but all current evidence points towards differences in the growth rates in the layers as the driving force behind cortical folding [11,[15][16][17]. In particular, it was recently demonstrated that the wrinkling instability also correctly captures variations of thickness between gyri and sulci [18].
Experimental verification of these theories is currently limited due to the difficulties involved in acquiring and mechanically testing brain matter in utero, but advances in non-invasive imaging techniques may provide the data needed to better validate their predictions [19,20]. Recently, it has been demonstrated that it is possible to capture the mechanical response of brain tissue in an elasticity-based framework [21]. Furthermore, preliminary numerical simulations have been able to demonstrate-at least phenomenologically-brain morphogenesis in this framework (see [22,23]).

(b) Mathematical studies of compression-induced instabilities
One of the first mathematical studies of wrinkling in elastic solids came from Biot in his seminal paper of 1963 [24]. In this work, he considers an incompressible elastic half plane under uniform compression of two different types, each of which induces a surface instability-at some critical

The model (a) General formulation
The basis of our computations is a three-dimensional formulation similar to those presented in [11,36]. We can substantially simplify the problem by only considering two-dimensional deformations, which is achieved by assuming that the material is in plane strain-that there is both no displacement in the transverse dimension and no dependence of the other components of the displacement on the spatial coordinate in that dimension. We consider the following model, illustrated in figure 1: let the region B s represent the initial unstressed infinite elastic substrate and B f be an elastic film bonded to its upper surface. Together, these form the domain B = B f ∪ B s . Let μ s and μ f represent the shear moduli of their respective layers, β = μ f /μ s be their ratio and X be a coordinate system across the two layers in the reference configuration. Let us henceforth fix our  Figure 1. Geometry of the domain. The system is composed of a bilayer with an infinitely deep layer of width 2L bonded by a film of thickness 1. Considering only plane strain, the problem is reduced to the deformation of a two-dimensional system under either compression or growth causing wrinkling. The boundary conditions are: continuity of traction and displacement between the layers, sliding vertical boundaries with no detachment, traction free upper layer and no displacement or detachment at infinity. (Online version in colour.) the new material coordinates of the deformed configuration are given by x(X) with deformation gradient We consider the two extreme cases that we label growth and compression.

(b) Growth
The fundamental assumption which allows us to incorporate material growth into the framework of elasticity comes from the theory of morphoelasticity, as described in [37]. We assume that any residual stresses within the material in the absence of applied loads are the result of growth on a local level and hence that we can decompose the deformation gradient multiplicatively as where A is an elastic deformation tensor and G is a growth tensor describing the local effects of growth. The application of the growth tensor alone to the reference configuration may not produce a physically realisable body, but the following application of the elastic tensor introduces stresses that enforce the boundary conditions and remove unphysical phenomena such as selfintersection. If we assume that our material is incompressible, we impose the constraint det A = 1 almost everywhere-the only local changes in volume of the material come from the growth process. For a hyperelastic material, with elastic strain-energy density function W, we can define an augmented energy density functional for the composed deformation by (see [38]) Here, p is a Lagrange multiplier that imposes the incompressibility constraint. In the particular case of a neo-Hookean material that we use here, the strain-energy density is given by

(c) Growth and compression
We can also consider an additional lateral compression in addition to or in place of growth in our system. As with growth, we can specify an initial diagonal stretch tensor A 0 to prescribe the external stretches that are applied to the bilayer. Since no new material is generated in this process, we must have det A 0 = 1. Our multiplicative decomposition is now Since A 0 represents an elastic process, our energy density functional (2.3) is unchanged. Indeed, we only separate it from A for notational convenience. We will study initial stretch tensors A 0 given by and a growth tensor G satisfying det G = J. In our case, we fix G = gI, where We then have J = g 2 . In our study, we will specialize to the compression-only case by taking J=1.
We can now specify our mathematical problem. We look for deformations x that are local minimizers of the total elastic energy of our system, subject to the elastic incompressibility constraint det A = 1. More explicitly, given G and A 0 , our variational problem is where Λ is a suitable Lagrange multiplier space that allows us to impose the pointwise constraint on A. The Euler-Lagrange equation for this system yields a necessary condition on minimizers of the energy.

(d) Stream functions
Taking advantage of the two-dimensional nature of the problem, we can make use of a technical tool to automatically satisfy the elastic incompressibility constraint. Essentially, we can find a stream function for the deformation, which is named after a similar construction used for the twodimensional Stokes flow. The difference here is that the domain of the stream function is a mixed coordinate space-it is a function of coordinates in both the reference and deformed configurations. The idea was first proposed in this setting by Rooney & Carroll [39] and used in [11,40]. Let x(X) = (x(X, Y), y(X, Y)) be any two-dimensional deformation for which where J is constant in each subdomain of B. More general growth conditions can be incorporated into this formulation, but for simplicity we will only study the constant case. Any such F can be decomposed multiplicatively as in (2.2). Away from some pathological cases, we can use an implicit function theorem based argument to define a function Ψ on the mixed coordinates (x, Y) such that From these representations, we can also compute the partial derivatives found in F = ∂x ∂X to rewrite the deformation gradient as Explicitly computing the determinant of F, we find hence the determinant constraint (2.10) is automatically satisfied exactly. To translate our energy functional into this stream function formulation, we make the change of integration variables (2.14) Since det A = 1 by construction, the Lagrange multiplier term in (2.9) disappears and leaves us with the minimization problem where Φ is the set of admissible stream functions.
To obtain the Euler-Lagrange equation of (2.15) and its boundary conditions explicitly, we must compute the first variation of its integral functionalĨ. For notational simplicity, we rewritẽ I asĨ (2.16) The Euler-Lagrange equations for the system are then given by The physical constraints we impose on the system at the boundaries are illustrated in figure 1.
We impose that all displacements in the substrate vanish at infinity, that the material may slide along-but not penetrate-the vertical sides of the domain and that the upper surface of the film is traction free. We define separate stream functions Ψ f and Ψ s for each layer of the system and we seek to simultaneously solve for the energy-minimizing stream function of each layer. The problems for each layer are coupled by the introduction of boundary conditions at the layer interfaces that impose continuity of traction and displacement between layers. From the physical condition that the two layers cannot detach from one another, we obtain, at Y = 0: On top of film (Y = 1), we obtain the traction-free conditions through the same process: Finally, we impose the decay conditions The two fourth-order PDEs for Ψ f and Ψ s in (2.17) and the eight boundary conditions given by (2.18)-(2.21) form the full Euler-Lagrange system. It should be noted that the explicit form of these Euler-Lagrange equations and their boundary conditions are lengthy with significant nonlinearity, making their direct solution impossible through analytic means.

(f) Perturbation
Despite the difficulties that a complete characterization of solutions to this problem presents, it is easy to see that the homogeneous growth solution given by with corresponding stream functions is a solution of (2.17). Consider a perturbation of the form Ψ = Ψ (0) + Ψ (1) , where is a small positive parameter. To linear order in , the Euler-Lagrange equations for the system read with boundary conditions given explicitly in appendix A. Assuming a periodic decomposition of the form Ψ (1) (x, Y) = sin(kx)h (1) (Y) for some k > 0, we arrive at the ODEs  27) in the case λ = γ and otherwise. Substituting these expressions into our boundary conditions, we obtain a homogeneous system of six linear equations in the six unknown coefficients c : where M is a 6 × 6 matrix. This system will only have non-trivial solutions if det M = 0, thus giving us a solvability condition for our system.

Linear analysis
We now focus our attention on two specific cases: a bilayer that is compressed unilaterally but experiences no growth and a bilayer that is under no compression but has a growing upper layer. The determinant of M is sufficiently complex that its zero level set cannot be obtained in closed form. However, it can be obtained asymptotically for short and long wavelengths and solved numerically in the intermediate regime.

(a) Compression
In the case of pure compression, we set γ = 1 and consider λ as our bifurcation parameter. The determinant of M can be written in the form where each p i is a polynomial in its arguments and For large values of k, exp(kζ 1 ) is the dominant term and thus, p 1 must vanish in order for the determinant to vanish in that limit. This polynomial-which has total degree 34-vanishes whenever λ is equal to either a particular root of the equation given by or a particular root λ * (β) of a polynomial given by the equation  is always strictly less than λ biot . Thus, λ biot provides a lower bound for the critical compression ratio required to cause the emergence of non-trivial solutions.
To better understand the solution set, we can solve the determinant relation numerically. We can fix a stiffness ratio β and find the compression ratio λ as a function of the wavenumber k. An example of such a dispersion curve is shown in figure 2 for the particular value β = 10. From this, we can deduce that if we were to gradually decrease the compression ratio λ from 1, we would expect to see non-trivial periodic solutions emerging at λ cr ≈ 0.89 with wavenumber k cr ≈ 0.61. We can repeat this process and track the position of this critical point as we vary the value of β, as shown in figure 3. As β decreases towards 1, λ cr approaches λ biot . For values of β infinitesimally above 1, a finite wavenumber k ≈ 0.941 is selected, but at β = 1, all wavenumbers are possible. For β < 1, we see the reappearance of a critical point, but it is in fact a local minimum rather than a local maximum. Hence, surface instability appears first for all values of β < 1.
For large values of β corresponding to a stiff film on a soft substrate, the selected wavenumber becomes vanishingly small and the critical compression ratio on an infinite domain approaches 1, corresponding to the Euler buckling instability. A standard asymptotic analysis reveals the following approximations (illustrated in figure 3): We recover the well-known dependence for the wavelength with a β 1/3 scaling that was already established by Biot [41] and has been recovered numerous times since then (see [42], for example).

(b) Growth
The growth case displays many similarities to the compression case. Considering large values of k once more reveals the existence of a Biot-type wrinkling instability for the system as described in figure 1. As before, the determinant can be written in the form  where each p i is a polynomial in its arguments and (ζ i ) 4 i=0 = (0, 1 + γ 2 , −1 + γ 2 , −1 − γ 2 , 1 − γ 2 ). For large enough k, exp(kζ 1 ) is the dominant term and hence in order for the determinant to vanish in that limit, p 1 must vanish. We find that polynomial p 1 vanishes whenever γ is equal to either a particular root of the equation given by or a particular root γ * (β) of the equation Further examination reveals that we have γ * (β) > γ biot for all values of β > 0. Thus, γ biot provides an upper bound on the critical growth factor required in order to achieve non-trivial periodic solutions. Solving the determinant relation numerically once more, we can fix a stiffness ratio β and find the growth factor γ as a function of the wavenumber k. Examples of such dispersion relations are shown in figure 4.
In a thought experiment where we gradually increase the growth factor γ from 1, we expect to see nontrivial periodic solutions emerging at γ cr . As before, by repeating this process, we can track the position of this critical point as we vary the value of β, as shown in figure 5. As β decreases, γ cr approaches γ biot and we see that the value of k cr increases without bound, demonstrating the aforementioned instability. The value of β = β min at which the wavenumber first diverges can be found exactly (but is not given analytically here) and is β min ≈ 1.90379.   Euler-type (large wavelength) and Biot-type (small wavelength) instabilities in the system. To this end, we repeat the linear analysis found in figure 3, adding additional insights where necessary.
Since the method has already been described at length, we briefly explain the new aspects of the problem without details.

(a) Bilayer with surface tension
The first modification we consider is the addition of surface energy. In elastic solids, there is an energetic cost to maintaining a surface that we must incorporate into our variational formulation when the material is sufficiently soft or to model the effect of a small layer on top of the material surface. To do this, we add another term to the energy functional in (2.9) to represent the surface energy at the interface between the layers and/or at the top of the upper layer. Following [11], this contribution takes the form where d is a surface energy density and Γ is a subset of ∂B f ∪ ∂B s . The addition of this term has no effect on the bulk Euler-Lagrange equations (2.17), but instead modifies the boundary conditions. In particular, (A 9a) becomes This extra term adds dependence on the surface energy parameter d to the system of linear equations (2.29) so that it is now of the form As before, this homogeneous system of linear equations has non-trivial solutions precisely when the determinant of the matrixM vanishes.

(i) Compression
In the compression case (γ = 1), we can write the determinant in the form where eachp i is some polynomial in its arguments and (ζ i ) 4 i=0 = (0, λ −1 + λ, −λ −1 + λ, λ −1 − λ, −λ −1 − λ). First, we remark that k = 0 is always a solution for λ = λ biot . Second, for large values of k, we have again that exp(kζ 1 ) is the dominant term and hence for large k,p 1 must vanish. There is no longer a zero of this polynomial at λ biot for all k, but there is still one at λ * (β) (for β = 1). For this root, we have λ * (β) < λ biot for all β > 0 as shown in figure 6. Hence, we conclude that as λ decreases, it eventually reaches λ biot at k = 0 which becomes the first instability.
When we compute the position of the critical growth and wavenumber as a function of β −1 (plotted in figure 7), we see a dramatic change in the qualitative behaviour of both quantities. Firstly, we see the disappearance of the critical point for values of β 2.1. However, the critical point ceases to be a global maximum before this occurs: for values of β 2.6 the global maximum of the dispersion curve occurs at k = 0 with a selected compression ratio of λ biot . Hence, the addition of surface tension prevents the Biot instability from occurring, replacing it with a Eulertype instability: if the film is sufficiently soft then the whole system buckles in a similar manner to a beam instead of displaying periodic fine wrinkling.
With the addition of another parameter, we can also fix the value of β and track the change in the critical growth and wavenumber as d varies. As one might expect, figure 8 demonstrates that the higher the surface energy density, the lower the compression ratio required to induce wrinkling and the lower the wavenumber of the wrinkling.  As before, for large values of β corresponding to a stiff film on a soft substrate, the selected wavenumber becomes vanishingly small and the critical compression factor approaches 1. A standard asymptotic analysis reveals the following approximations (illustrated in figure 7) that demonstrate the influence of the surface energy parameter d on the critical growth factor and In the growth case, we can write the determinant in the form where eachp i is some polynomial in its arguments and (ζ i ) 4 i=0 = (0, 1 + γ 2 , −1 + γ 2 , −1 − γ 2 , 1 − γ 2 ). Consideration of the dominant term in the large k limit yields an asymptote at γ = γ * (β), which approaches γ biot from above in the large β limit. As in the compression case, this asymptote is independent from d. The dispersion curve is similar to figure 2 with γ * (β) replacing γ biot . Echoing the results from the compression case, the critical growth factor is significantly increased, but occurs at a smaller wavenumber.
When we compute the position of the critical growth and wavenumber as a function of β −1 (plotted in figure 9), we see a dramatic change in the qualitative behaviour of both quantities. In particular, we no longer see a blow-up in the wavenumber as β decreases and we see an apparent increase in γ cr without bound. However, the critical point that we are computing stops being a global minimum of γ for sufficiently small values of β. For β under this threshold, we would again expect a Biot-type instability.
Plotting the critical growth factor and wavenumber as a function of d (shown in figure 10) reveals that the higher the surface energy density, the higher the growth factor required to induce wrinkling and the lower the wavenumber of the wrinkling.
Another standard asymptotic analysis for large values of β gives the following approximations (illustrated in figure 9) for the correction that the surface energy parameter d induces on the critical growth factor: (1) Following the same method as before, we obtain a homogeneous linear system of eight equations in the eight unknownsc : . This leads to the solvability conditioñ M(k, γ , λ, β f , β t )c = 0, (4.12) which only has nontrivial solutions if detM = 0.

(i) Compression
In the compression case, on numerically plotting the solution set of the determinant relation, we find a similar dispersion curve compared to the unmodified problem (figure 2) but with λ biot replaced by an asymptote λ * that depends on both β f and β t . The addition of another elastic layer decreases the critical compression ratio and the compression threshold for large k while increasing the critical wavenumber. For a given, fixed stiffness ratio β f β −1 t , from figure 11 we can see that as β f and β t decrease, the critical compression ratio approaches the previously discussed threshold. However, we now find that for β < 1, we have a Euler-type buckling instability where the wavenumber k = 0 is selected.
By contrast, figure 12 shows us that if we fix β f and vary β t , we observe a gradual increase in λ cr and a gradual decrease in k cr as β t decreases, with the wavenumber remaining well determined. Thus, the addition of an upper layer decreases the critical compression ratio and increases the critical wavenumber selected in the system. Finally, figure 14 demonstrates that if we again fix β f and decrease β t , γ cr and k cr both decrease with no apparent blow-up behaviour. Hence, the addition of an upper layer increases the critical growth factor and critical wavenumber selected in the system. In particular, as the stiffness of the upper layer approaches that of the film from below, the critical growth and wavenumber increase without bound.

(c) Pressure
We can derive the effect of a normal pressure of magnitude p acting on top of the film layer by directly imposing this constraint on the surface. A pressure p on the surface can be expressed in terms of the Cauchy stress tensor T as T · n = pn for all points on the top surface. To express this condition, we compute, to first order, the normal vector field and the Cauchy stress. First, we recall that for a sufficiently regular deformation, the unit surface normal vector in the deformed configuration is given by  periodic decomposition, it reads: Similarly, the traction boundary conditions at the film-substrate interface become (4.26b)

(i) Compression
In the compression case, the addition of fibres initially seems to have a limited effect. Compared to the unmodified case with the same large stiffness ratio, adding fibres with 0 < m < 1 causes the critical wavenumber to increase, the critical compression ratio to decrease and has no effect on the large k asymptote. However, as we vary the stiffness ratio, we see some markedly different behaviour in the evolution of the critical point as a function of β (demonstrated in figure 15). For stiffness ratios β < β c ≈ 8, the critical wavenumber rapidly decreases. This local maximum close to k = 0 persists even when β > 1 (at β = 1 we have no length scale and the wavelength is again undetermined). The position of the critical compression ratio does not appear to degenerate to λ biot for small values of β.
Increasing the fibre stiffness parameter m (shown in figure 16) has a similar effect; for fibres stiffer than m = m c ≈ 0.67, the critical wavenumber becomes close to k = 0.
To summarize, for fibres significantly stiffer than the elastic substrate in which they are embedded, we see a lower wavenumber wrinkling pattern emerge.

(ii) Growth
In the growth case, the addition of fibres causes both the critical wavenumber and the critical growth factor to increase slightly but the large k behaviour of the system is unchanged. The critical growth factor and wavenumber have a similar qualitative behaviour compared to the unmodified  case with the notable characteristic that the stiffness ratio 1/β min at which the wavenumber blows up is significantly reduced. Plotting the dependence of the critical growth factor and wavenumber on the fibre stiffness parameter m as in figure 17 shows a gradual increase in both quantities as the fibres become stiffer. For a fixed β, there exists a finite (but extremely large) m such that γ cr = γ biot and k cr becomes infinite.

Conclusion
We have presented a complete linear analysis for the plane-strain wrinkling of a film on an elastic substrate in the case of lateral compression and film growth. The analysis does not make any approximation on the thinness of the film or the relative stiffness ratio between substrate and film. Hence, it can be used as a general benchmark for approximate theories and identify their domain of validity. We also considered the role of secondary effects such as surface tension, pressure and fibres. Our analysis further establishes that for films that are much stiffer than the substrate, a regular asymptotic expression in powers of 1/β leads to accurate predictions for the critical parameter and critical wavenumber selected at the wrinkling instability even when supplementary effects are considered. A rule of thumb is that for β 10, a 3-term expansion is sufficient in all cases to capture the correct behaviour. It also suggests that in this regime, approximate theories (beams and plates) may be sufficient as long as they correctly model the effect of the substrate. Our analysis can be used to gauge this calibration by matching the asymptotic behaviours of a plate or beam to the ones derived here.
As β decreases, a number of different effects appear that make general conclusions harder to reach. Depending on both the loading and the effect considered, qualitatively different behaviours are observed. For instance, the addition of any surface tension in compression changes the Biot surface instability (k cr → ∞) to a Euler-type instability (k cr → 0). Yet, the Biot instability is still the first selected for a growing film. Similarly, the minimal value β min at which a linear instability is found depends greatly on both the loading and extra surface effects. It is, therefore, harder to obtain a general picture for the bifurcation of soft films on substrate. Yet, the linear analysis may not even be relevant in that regime for two reasons.
First, the film may undergo a creasing instability for values of the axial stretch around λ = 1/γ ≈ 0.64 [35]. Hence, the linear unstable wrinkling mode may not be observed past that critical value. Whether this instability is universally observed in bilayers and always selected is still an open problem.
Second, the analysis performed here is only linear and does not allow us to conclude about the existence of periodic solutions past the bifurcation point. The main problem is that the wrinkling instability may be supercritical or subcritical depending on the stiffness ratio [9,[45][46][47][48]. Previous studies suggest that for sufficiently stiff films, the wrinkling instability is supercritical. The question is then to determine the value of β at which this supercritical bifurcation becomes subcritical and whether this value occurs before or after the Biot instability or the wrinkling instability. However, to answer this question and have a full picture of the wrinkling instability requires a full weakly nonlinear analysis of the bifurcation we have identified here. We leave this exercise for the second instalment of this work.
Data accessibility. This article has no additional data. Competing interests. The authors declare that they have no competing interests.

Funding. The support for Alain Goriely by the Engineering and Physical Sciences Research Council of Great
Britain under research grant no. EP/R020205/1 is gratefully acknowledged. The research of Hamza Alawiye was supported by UK Engineering and Physical Sciences Research Council grant no. EP/L015811/1. This work was supported by the NSF grant no. CMMI 1727268 to E.K.