Liquid crystal elastomers wrinkling

When a liquid crystal elastomer layer is bonded to an elastic layer, it creates a bilayer with interesting properties that can be activated by applying traction at the boundaries or by optothermal stimulation. Here, we examine wrinkling responses in three-dimensional nonlinear systems containing a monodomain liquid crystal elastomer layer and a homogeneous isotropic incompressible hyperelastic layer, such that one layer is thin compared to the other. The wrinkling is caused by a combination of mechanical forces and external stimuli. To illustrate the general theory, which is valid for a range of bilayer systems and deformations, we assume that the nematic director is uniformly aligned parallel to the interface between the two layers, and that biaxial forces act either parallel or perpendicular to the director. We then perform a linear stability analysis and determine the critical wave number and stretch ratio for the onset of wrinkling. In addition, we demonstrate that a plate model for the thin layer is also applicable when this is much stiffer than the substrate.


Introduction
Liquid crystal elastomers (LCEs) are complex materials that combine the elasticity of polymeric solids with the self-organisation of liquid crystalline structures [31,38].Due to their molecular architecture, consisting of cross-linked networks of polymeric chains containing liquid crystal mesogens, these materials are capable of interesting mechanical responses, including relatively large nonlinear deformations which may arise spontaneously and reversibly under certain external stimuli (e.g.heat, solvents, electric or magnetic field).These qualities suggest many avenues for technological applications, such as soft actuators and soft tissue engineering, but more research efforts are needed before they can be exploited on an industrial scale [32,48,53,59,65,79,82,84,[90][91][92].
In particular, the occurrence of material microstructures and the formation of wrinkles in response to external forces or optothermal stimulation are important phenomena that require further investigation.A survey on surface wrinkling in various materials is provided in [57].The utility of surface instabilities in measuring material properties of thin films is reviewed in [23].Within the framework of elasticity, relevant experimental and theoretical studies of wrinkling responses in LCEs are as follows: (I) Wrinkling and shear stripes formation in LCEs: (I.1) Experimental results.In [39,54,74,99], stretched LCE samples were reported to display director re-orientation and microstructural shear stripes without wrinkling near the clamped ends, unlike in purely elastic samples.In [1,2], spontaneous formation and reorientation of surface wrinkles under temperature changes were observed for stiff isotropic elastic plates bonded to a monodomain nematic elastomer foundation.In this case, wrinkles formed either parallel or perpendicular to the nematic director, depending on the temperature at which the bilayer was prepared.It was further suggested that controlled surface wrinkling patterns could be employed to characterise the mechanical properties of thin polymer films deposited on the LCE.Surface wrinkles in thin liquid crystal polymer films on a glass substrate were reported in [49].In [83], snake-mimicking soft actuators composed of a bilayered LCE ribbon were described.(I.2) Models based on the linear elasticity theory.In [22,67], the formation of shear stripes microstructure and fine-scale wrinkles in stretched LCEs were studied using a Koiter-type theory.Namely, for physically relevant parameters, it was shown that the microstructure was finer than wrinkles, and that these instabilities occurred for distinct mesoscale stretches.For LCE plates with or without a compliant substrate, in [52], such instabilities were considered within the theoretical framework of Föppl-von Kármán.In [70], small-amplitude wrinkles of a compressed isotropic elastic film on a soft nematic elastomer substrate were treated.The formation and evolution of wrinkling in a solid liquid crystalline film attached to a compliant elastic substrate was explored numerically in [43,[93][94][95]98].Finite element simulations of shear striping were presented in [97].
(II) Wrinkling of isotropic elastic materials: (II.1) Models based on the linear elasticity theory.Wrinkling of a stretched isotropic elastic membrane, with zero bending stiffness, was examined in [66,71,72].A stretched Koiter-type plate was considered in [75].In [6-8, 16, 24, 46], a Föppl-von Kármán plate bonded to a compliant infinitely deep linearly elastic foundation under equibiaxial compression was analysed.Wrinkling of twisted ribbons was addressed in [50], inspired by [9].Large-amplitude wrinkling of a thin linearly elastic film attached to a pre-stretched substrate modelled by an Ogden strain-energy function was treated in [73].Multi-layer systems were studied in [55].(II.1)Models based on finite strain elasticity, with either compression or pre-stretch causing wrinkling.A numerical study of a stretched neo-Hookean sheet was presented in [65].
The first linear stability analysis of a nonlinear elastic half-space under compression in a direction parallel to its surface was developed in [13].The stability of the wrinkling mode experienced by a compressed half-space of neo-Hookean material was investigated in [20].In [17], the first nonlinear analysis of a bilayer system formed from a stiff nonlinear hyperelastic film bonded to a nonlinear hypelastic substrate and subjected to uniaxial compression parallel to the interface was provided.A two-layer system comprising a neo-Hookean film bonded to an infinitely deep neo-Hookean substrate, with the film/substrate stiffness ratio ranging from the limit in which the film and substrate have the same elastic modulus to a very stiff film on a compliant substrate, subjected to a combination of compression and pre-stretched, was examined in [21].The influence of the relative stiffness in a two-layer system of neo-Hookean materials was explored in [47].In [42], a uniaxially compressed neo-Hookean thin film bonded to a neo-Hookean substrate of comparable stiffness was studied.In [18,41], the exact nonlinear elasticity theory combined with an asymptotic perturbation procedure was employed to derive the critical stretch value at which period-doubling secondary bifurcation under uniaxial compression occurred.
(III) Wrinkling in growing biological tissues: (III.1) Models based on plane strain finite elasticity, where the physical system is reduced to a two-dimensional problem, with either compression or growth causing wrinkling.
A semi-analytical approach for the study of the nonlinear buckling behaviour of a growing soft layer was developed in [27].In [15], a model for a growing thin neo-Hookean layer bonded to an infinitely deep neo-Hookean substrate was proposed, where the stiffness ratio between the layer and substrate was relatively low while the thin layer was stiffer.A system formed from a neo-Hookean layer attached to a stiffer or softer neo-Hookean substrate undergoing differential growth was treated in [10].Bilayer film-substrate systems of neo-Hookean materials with different stiffness ratios, subjected to combined compression and growth conditions were studied systematically in [45].In [56], multi-layer systems were considered.A linear stability and a weekly nonlinear post-bifurcation analysis of a growing neo-Hookean film on a stiffer or softer neo-Hookean substrate, subjected to lateral compressive forces acting on the whole system were presented in [3,4], respectively.Notably, [3,4,10,27] employed a stream function [11,19,25] operating on a mixture of the reference and deformed coordinate systems to replace the use of Lagrange multipliers associated with the incompressibility constraint.(III.2) Models based on the three-dimensional finite elasticity.For a growing soft layer subjected to equi-biaxial compression, in [26], a weakly nonlinear analysis of the wrinkling instability was performed using the multiple-scale perturbation method applied to the incremental theory in finite elasticity.
At the constitutive level, for ideal monodomain LCEs, where the mesogens are uniformly aligned throughout the material, a general formulation is provided by the phenomenological neoclassical strain-energy function proposed in [14,86,89].This model is based on the molecular network theory of rubber elasticity [77], and its parameters are directly measurable experimentally or derived from macroscopic shape changes [87,88].Here, we adopt the stressfree state of the isotropic phase at high temperature as the reference configuration [28,[33][34][35][36][60][61][62], rather than the nematic phase in which the cross-linking was produced [5,14,81,85,86,89,96].Within this theoretical framework, the material deformation can be expressed as a composite deformation from a reference configuration to the current configuration, via an elastic distortion followed by a natural (stress free) shape change.The multiplicative decomposition of the associated gradient tensor is similar to those found in the constitutive theories of thermoelasticity, elastoplasticity, and growth [44,58], but it is fundamentally different since, for nematic materials, the stress free geometrical change is superposed on the elastic deformation, which is applied directly to the reference state.This difference is important since, although the elastic configuration obtained by this deformation may not be observed in practice, it may still be possible for the nematic body to assume such a configuration under suitable external stimuli.The resulting elastic stresses can then be used to analyse the final deformation where the particular geometry also plays a role [61,62].
In this study, we examine nonlinear elastic wrinkling of three-dimensional (3D) bilayer systems containing a monodomain LCE layer and a homogeneous isotropic incompressible hyperelastic layer.First, we present a general linear instability analysis, which is valid for a range of bilayer systems and deformations.We then apply this analysis to a biaxial stretch combining elastic and natural deformations.For rubber-like material, experimental testing under biaxial (also known as 'pure shear') stretch has been carried out since the original work by Rivlin and Saunders (1951) [68].When a nematic LCE sample is subjected to biaxial stretch, depending on the stretch ratios, the deformation induces phases of shear striping or wrinkling [22,30,34,36,60].We consider the following two cases: a film of hyperelastic material attached to a LCE substrate, and a LCE film bonded to a hyperelastic substrate.In each case, assuming that the nematic director is uniformly aligned parallel to the film plane, and that biaxial forces act either parallel or perpendicular to the director field, we determine the critical wave number and stretch ratio for the onset of wrinkling.We find that: (a) When the LCE is at high temperature then cooled, it extends along the director and contracts in the perpendicular direction causing wrinkling parallel to the director, as sketched in figure 1; (b) When the LCE is at low temperature then heated, it contracts along the director and extends in the perpendicular direction causing wrinkling perpendicular to the director, as sketched in figure 2; (c) At constant temperature, if the nematic layer is compressed along the director, it can lead to reorientation of the director until it aligns perpendicular to the compressive force, after which, further compression produces wrinkling parallel to the rotated director.
In section 2, we present the constitutive description for LCEs.In section 3, we introduce the bilayer system and the minimisation problem for its equilibrium state.Section 4 is devoted to the infinitesimal incremental perturbation of the displacement field with respect to a homogeneously deformed state.Critical values for the formation of wrinkles in specific bilayer systems under biaxial stretch are obtained in section 5.In section 6, we briefly consider the case when the film deformation is approximated using a plate theory, and demonstrate the validity of this approximation when the film is much stiffer than the substrate.In the final section, we draw concluding remarks and outline some potential future developments.Wrinkling parallel to the nematic director n when a > 1 and the system is at high temperature then cooled, so it extends along the director and contracts in the perpendicular direction [1,2].Wrinkling perpendicular to the nematic director n when a < 1 and the system is at low temperature then heated, so it contracts along the director and extends in the perpendicular direction [1,2].

The neoclassical model
The neoclassical strain-energy density function describes an ideal LCE and takes the form where F denotes the deformation gradient from the isotropic state, n is a unit vector field, known as the director, and W(A) represents the strain-energy density of the isotropic polymer network, depending only on the (local) elastic deformation tensor A. The tensors F and A satisfy the relation (see figure 3) where is the spontaneous deformation tensor defining a change of frame of reference from the isotropic phase to a nematic phase.Note that, the spontaneous tensor G given by ( 3) is symmetric, i.e.G = G T (where the superscript 'T' denotes the transpose operation), but the elastic tensor A may not be symmetric.In (3), a > 0 is the temperature-dependent stretch parameter, ⊗ is the usual tensor product of two vectors, and I = diag(1, 1, 1) is the identity tensor.We assume that a is spatially-independent (i.e.no differential swelling).This parameter can be estimated by examining the thermal or light-induced response of the nematic elastomer with uniform planar alignment [51].The ratio r = a 1/3 /a −1/6 = a 1/2 represents the anisotropy parameter, which, in an ideal nematic solid, is the same in all directions.In the nematic phase, both the cases with r > 1 (prolate molecules) and r < 1 (oblate molecules) are possible, while when r = 1, the energy function reduces to that of an isotropic hyperelastic material [32].
Monodomains can be synthesised with the parameter a taking values from 1.05 (glasses) to 60 (nematic rubber), which would correspond to changes in natural length between 7% and 400% (spontaneous extension ratio a 1/3 between 1.02 and 4) [29].
In (1), the elastic strain-energy function W is minimised by any deformation satisfying AA T = I [64,78], whereas the nematic strain-energy W (nc) is minimised by any deformation satisfying FF T = G 2 .Hence, every pair (GR, n), with R an arbitrary rigid-body rotation (i.e.R −1 = R T and det R = 1), is a natural (i.e.stress free) state for this material model.By (3), the following identity holds The director n is an observable (spatial) quantity.Denoting by n 0 the reference orientation of the local director corresponding to the cross-linking state, n may differ from n 0 both by a rotation and a change in r.We confine our attention to the case when the hyperelastic strain-energy function in (1) describes an incompressible neo-Hookean material [76], i.e.
where 'tr' denotes the trace operator, and μ > 0 represents the constant shear modulus at small strain.
For a finite elastic deformation with gradient tensor A, the left and right Cauchy-Green tensors are defined, respectively, by Using these deformation tensors, the elastic Almansi strain tensor is equal to [64, pp 90-91] e = 1 2 and the elastic Green-Lagrange strain tensor is [64, pp 89-90] Based on the hyperelastic model defined by ( 5), we construct the following strain-energy function of the form (1), where we assume that the nematic director is 'free' to rotate relative to the elastic matrix, i.e.F and n are independent variables [28,34,36].For this case, the relations between the stress tensors for the LCE described by ( 9) and the base polymeric network modelled by ( 5) are presented in appendix A (see also [61]).In particular, the following additional condition related to the general lack of symmetry of the Cauchy stress tensor in LCEs holds [5,96], We note that many macroscopic deformations of LCEs induce a director re-orientation, such that the director becomes parallel to the direction of the largest principal stretch.

The nonlinear bilayer system
We now consider a system composed of two solid layers, namely a thin film and a thick substrate, attached continuously to each other across a contact interface.In the undeformed reference configuration, each body occupies a compact domain Ωk ⊂ R 3 , k ∈ {f, s}, with the index 'f' for the film and 's' for the substrate, such that its interior is an open, bounded, connected set Ω k ⊂ R 3 and its boundary ∂Ω k = Ωk \Ω k is Lipschitz continuous (in particular, we assume that a unit normal vector N exists almost everywhere on ∂Ω k ).We denote by Ω = Ω f ∪ Ω s the domain occupied by this system, such that We denote by X = (X 1 , X 2 , X 3 ) and x = (x 1 , x 2 , x 3 ) the Cartesian coordinates for the reference and the deformed state, respectively.We then designate the plane formed by the first two directions as the film's plane and the third direction as the film's thickness direction, and take such that the thickness of the film is much smaller than that of the substrate, i.e. h H.We assume that each layer is made of either a purely elastic material characterised by a neo-Hookean strain-energy function of the form (5), or a LCE described by a neoclassical strain-energy function of the form (9). Therefore, in the reference configuration, the material in each layer is stress-free and isotropic.We denote by μ f and μ s the respective constitutive coefficients, and their ratio by β = μ f /μ s .
For the two-layer system, the kinematically admissible displacement fields u = u(X) = x − X have gradient tensor ∇u = F(u) − I, and satisfy the incompressibility condition det F(u) = 1, the essential boundary conditions on Γ D , and the continuity condition [u] = u + − u − = 0 across the interface Γ C .In addition, for the LCE layer, the kinematic constraint (2) holds.
For the hyperelastic material, we define the function where W k (F(u)) takes the form ( 5) with μ = μ k and A = F(u) = I + ∇u, and p k is the Lagrange multiplier for the incompressibility constraint det Similarly, for the LCE, we define where given by ( 9) with μ = μ k and F = F(u) = I + ∇u, and p k is the Lagrange multiplier for the incompressibility constraint det The elastic energy stored by each individual body is equal to (see appendix B for details) with Ψ k described by either (11) or (12).Thus, the energy functional is E k = E k (u, p k ) for the hyperelastic material, and Assuming that the two layers are subjected to dead-load traction f N on Γ N ⊂ ∂Ω\Γ C , we are interested in the minimisation of total potential energy over all kinematically admissible displacement fields u = [u 1 , u 2 , u 3 ] T .
Recalling that (10) holds, the eight Euler-Lagrange equations for the minimisation problem in Ω k , k ∈ {f, s} are: where u i, j = ∂u i /∂X j , i, j = 1, 2, 3.These equations are complemented by the following conditions across the interface: , the jump in displacement is equal to zero, i.e.
[u] = 0, (19) and the surface tractions are equal in magnitude and opposite in orientation, i.e.
Equations ( 15)-( 18) and contact conditions (19) and ( 20) are completed by the following external boundary conditions: where λ 1 > 0 is fixed.• On the front and back surfaces where λ 2 > 0 is fixed.• The top and bottom surfaces In addition to the specified contact and boundary conditions, as a general rule, at a corner where one of the adjacent edges is subject to essential conditions and the other to natural conditions, the essential conditions take priority, and when both edges meeting at a corner are subject to natural conditions, these conditions are imposed simultaneously at the corner.Also, in order to avoid rigid-body rotations of the entire system, the lower left-hand corner is clamped, i.e. u = 0 at that corner.

Incremental perturbation
Next, we assume an infinitesimal incremental perturbation of the displacement field with respect to the homogeneously deformed state described by ( 24)-( 26).The perturbed displacement field takes the form where u (0) is the displacement corresponding to the homogeneously deformed state, u (k) is the infinitesimal incremental displacement, and 0 < 1 is a small parameter.The associated Lagrange multipliers are equal to For the incremental perturbation, the continuous contact and external boundary conditions are as follows: • The top surface • The bottom surface We search for solutions of the form where the prime denotes differentiation and q > 0 is the wave number.We first take a system consisting of a hyperelastic film and a LCE substrate [2].Then the opposite configuration of a LCE film atop a hyperelastic substrate can be directly obtained by swapping variables, as will be explained later.
From the Euler-Lagrange equations, we obtain the following fourth-order ordinary differential equations, where ζ k , ζ k , ζ k , and ζ iv k are the derivatives of order 1, 2, 3 and 4 of ζ k with respect to x 3 , respectively, and k ∈ {f, s}.The corresponding Lagrange multipliers take the form When λ 2 1 λ 2 = 1, the following general solutions satisfy equations ( 36) and ( 37), respectively, To determine the eight unknown coefficients, {c (f) i , c (s) i } i=1,...,4 , in the above equations, we substitute expressions (40) and (41) in the contact and external boundary conditions ( 29)- (34).We obtain the following homogeneous system of eight linear equations in the eight unknown coefficients: When the substrate is infinitely thick, the external conditions on the bottom surface Γ 3s are replaced by the assumption that u (s)  i → 0, i = 1, 2, 3, as x 3 → −∞.In this case, c (s) 3 = c (s) 4 = 0, and this homogeneous algebraic system reduces to where T and the entries of the 6 × 6 coefficient matrix M = (M i j ) i, j=1,...,6 are explicitly given in table 1.

Wrinkling under biaxial stretch
The general equations for the linear stability analysis in section 4 are valid for all bilayer systems and deformations described in section 3.As a particular application, we now examine the wrinkling of our bilayer system with an infinitely thick substrate, subject to biaxial stretch, such that [60] and The biaxial deformation described by (51) was previously addressed in [60,62] (see also [30,35,36]) where shear-striping pattern formation and the accompanying 'soft elasticity' phenomenon [80] were analysed in detail.We will rely those results for our subsequent analysis.
For a systematic treatment of the effect of biaxial forces acting either parallel or perpendicular to the nematic director, while the director remains parallel to the contact surface between the film and the substrate throughout deformation, we further require that which, together with (51), implies This condition prevents the director from rotating out-of-plane, and also assumes that the largest principal stretch is in the second direction.The former is a reasonable physical assumption, while the latter is imposed without restricting the generality of the problem if we consider the reference director to be uniformly aligned either parallel or perpendicular to the second direction.The orthogonal situation, with λ 1 > λ 2 , then leads to analogous results when the reference director aligns either parallel or perpendicular to the first direction.We are interested in the solutions λ = λ cr of equation that satisfy (54).
The associated critical second Piola-Kirchhoff stress tensors (see appendix B) then take the form cr a 2/3 , 0 for the elastic layer, ( 56) The general formulae for the above stress tensors are given by (A.3) and (A.8), respectively.

Hyperelastic film on LCE substrate
When the bilayer system is formed by a hyperelastic film and a LCE substrate, the entries of the 6 × 6 matrix M in (50) are given in table 2 where we take h = 1 without loss of generality.We are interested in solving equation (55).Fixing a, for each finite value of β, the solutions of det M = 0, viewed as an equation in λ, provide solution curves λ = λ(q).In compression, the critical stretch ratio is the maximal value of λ of this function.We refer to this extremal value and the value q where it is achieved as the pair (λ cr , q cr ) = (λ cr (β), q cr (β)).We have made the dependence on β explicit to indicate that we are primarily interested in the change of these critical values with the stiffness ratio β.Equivalently, the pair (λ cr (β), q cr (β)) is obtained as a solution of the system To find these values asymptotically, we first consider the case of very large β.In the limit ζ = 1/β → 0, the critical wave number is q cr = 0 and λ cr = λ 0 is obtained as a solution of det M = 0.Then, we substitute λ cr = λ 0 + λ 1 ζ ν 1 , q cr = q 1 ζ ν 2 and expand E in powers of ζ (assuming ν 1 and ν 2 positive) to find the correct balance of exponents (for our problem, we found the scaling ν 1 = −2/3 and ν 2 = −1/3).Once these are obtained, a regular asymptotic expansion in powers of ζ is substituted in E, and solved to arbitrary order.However, this solution is only valid locally near ζ = 0. To obtain the entire curve numerically for a given value of a, we fix β and solve E for λ and q simultaneously, by Newton's method, to derive their critical values.Explicitly, we start with a large value of β ≈ 100 for which the asymptotic solution provides an initial guess that is very close to the numerical root.Once the numerical solution with this value of β is known, we use it to solve for nearby values and, by regular continuation method, obtain the entire curve.Case 1. First, we assume that the director is initially aligned in the second direction, such that g 1 = a −1/6 and g 2 = a 1/3 .For all a > 0, in the asymptotic limit of large stiffness ratio β = μ f /μ s , we obtain the following critical stretch ratio λ cr and critical wave number q cr , respectively: For finite values of 0 < q < ∞, we can solve equation ( 55) numerically (for q → ∞, see appendix C). Figure 4 illustrates the critical stretch ratios λ cr and critical wave numbers q cr as functions of the relative stiffness ratio β: Table 2.The coefficients {M i j } i, j=1,...,6 of the homogeneous system of linear equation ( 50), for a hyperelastic film on a LCE substrate, with {λ i } i=1,2,3 given by ( 51) and {α i } i=1,2,3 given by (52).
• When a < 1, wrinkling parallel to the director is generated at the critical stretch ratio 0 < λ cr < a 1/3 ; • When a > 1, wrinkles parallel to the director is produced at the critical stretch ratio 0 < λ cr < a 1/12 .
For example, when a > 1 and the LCE is at high temperature then cooled, it extends along the director and contracts in the perpendicular direction, causing wrinkling parallel to the director (see figure 1).A similar behaviour produced experimentally is discussed in [1,2].
At constant temperature, microstructural shear striping patterns can also occur in the prolate LCE substrate if the director rotates [30,34,36,37,40,60,62,69,81,96].To verify whether shear striping is likely in our system, we assume that the LCE substrate is in its natural state before the film is attached, i.e. λ = a −1/6 in the reference configuration.Starting from the reference state with g 1 = a −1/6 and g 2 = a 1/3 , analysis similar to that carried out in [60,62] shows that this state is unstable and the director rotates, while shear stripes develop, as λ increases within the interval a −1/6 , a 1/12 .Then, decreasing the value of λ < a −1/6 will not produce rotation of the director in the substrate, i.e. no striping pattern will form.
Case 2. Second, we assume that the director is initially aligned in the first direction, such that g 1 = a 1/3 and g 2 = a −1/6 .For large stiffness ratio β = μ f /μ s , the asymptotic results read: q cr = a −2/9 a 1/2 + 1 2 1/3 3 Figure 5. Case 2 (elastic film on LCE substrate, with director aligned in the first direction): the critical stretch ratio 0 < λ cr < min a 1/12 , a 1/3 and critical wave number q cr as functions of relative stiffness ratio β = μ f /μ s when g 1 = a 1/3 and g 2 = a −1/6 .For different magnitudes of the nematic parameter a, the solid lines represent critical values, while the dashed lines show the value min a 1/12 , a 1/3 , respectively.
For finite values of 0 < q < ∞ (see appendix C for q → ∞), figure 5 presents the critical stretch ratios λ cr and critical wave numbers q cr as functions of the stiffness ratio β: • When a < 1, wrinkling is produced at the critical stretch ratio 0 < λ cr < a 1/3 .For example, when a < 1 and the LCE is at low temperature then heated, it contracts along the director and extends in the perpendicular direction, causing wrinkling perpendicular to the director (see figure 2).A similar behaviour obtained experimentally is discussed in [1,2].
• When a > 1, wrinkling perpendicular to the director is induced at the critical stretch ratio 0 < λ cr < a 1/12 .
At constant temperature, we assume that the LCE substrate is in its natural state before the film is attached, i.e. λ = a 1/3 in the reference configuration.This state is unstable and the director rotates, while shear stripes form, as λ decreases within the interval a 1/12 , a 1/3 [60].For λ < a 1/12 , ( 54) is satisfied and the problem reduces to that of case 1. Visual inspection of the critical stretch ratios in figures 4 (top-left) and 5 (top-left), and comparison of the asymptotic values given by ( 59) and ( 61) further suggest that the critical stretch ratios λ cr in case 1 are smaller than those in case 2. We interpret this as 'delay' in wrinkling due to director rotation.

LCE film on hyperelastic substrate
For a bilayer system consisting of a LCE film on a hyperelastic substrate, the entries of the 6 × 6 matrix M in (50) are given in table 3, again with h = 1.
−qλa 1/6 / g 2 1 g 2 0 0 Figure 6.Case 3 (LCE film, with director aligned in the second direction, on elastic substrate): the critical stretch ratio 0 < λ cr < min a 1/12 , a 1/3 and critical wave number q cr as functions of relative stiffness ratio β = μ f /μ s when g 1 = a −1/6 and g 2 = a 1/3 .For different magnitudes of the nematic parameter a, the solid lines represent critical values, while the dashed lines show the value min a 1/12 , a 1/3 , respectively.
Case 3. First, the director is initially aligned in the second direction, such that g 1 = a −1/6 and g 2 = a 1/3 .The asymptotic limit of large stiffness ratio β = μ f /μ s gives For finite values of 0 < q < ∞ (see appendix C for q → ∞), figure 6 displays the critical stretch ratios λ cr and critical wave numbers q cr as functions of the stiffness ratio β: • When a < 1, wrinkles parallel to the director form at the critical stretch ratio 0 < λ cr < a 1/3 ; • When a > 1, wrinkles parallel to the director are obtained at the critical stretch ratio 0 < λ cr < a 1/12 .
At constant temperature, the discussion is similar to that of case 1. Namely, if the prolate LCE film is in its natural state prior to being attached to the elastic substrate, then decreasing the value of λ < a −1/6 will not cause shear striping in the film.
Case 4. Second, the director is initially aligned in the first direction, such that g 1 = a 1/3 and g 2 = a −1/6 .In the asymptotic limit of large stiffness ratio β = μ f /μ s , we find Figure 7. Case 4 (LCE film, with director aligned in the first direction, on elastic substrate): the critical stretch ratio 0 < λ cr < min a 1/12 , a 1/3 and critical wave number q cr as functions of relative stiffness ratio β = μ f /μ s when g 1 = a 1/3 and g 2 = a −1/6 .For different magnitudes of the nematic parameter a, the solid lines represent critical values, while the dashed lines show the value min a 1/12 , a 1/3 , respectively.
q cr = a 2/3 + a 1/6 2 1/3 3 For finite values of 0 < q < ∞ (see appendix C for q → ∞), figure 7 shows the critical stretch ratios λ cr and critical wave numbers q cr as functions of the stiffness ratio β: • When a < 1, wrinkles perpendicular to the director appear at the critical stretch ratio 0 < λ cr < a 1/3 ; • When a > 1, wrinkles perpendicular to the director form at the critical stretch ratio 0 < λ cr < a 1/12 .
At constant temperature, the discussion is similar to that of case 2.Then, visual inspection of the critical stretch ratios in figures 6 (top-left) and 7 (top-left), and comparison of the asymptotic values given by ( 63) and (65) suggest that the critical stretch ratios λ cr in case 3 are smaller than in case 4, which can be interpreted as wrinkling being 'delayed' while the director rotates.

Approximate results based on plate models
In [61], we derived a Föppl-von Kármán-type constitutive model for a liquid crystalline solid which is sufficiently thin so that it can be approximated by a plate equation.In our derivation, we relied on the following assumptions: (i) surface normals to the plane of the plate remain perpendicular to the plate after deformation; (ii) changes in the thickness of the plate during deformation are negligible; (iii) the stress field in the deformed plate is parallel to the midsurface.For a wrinkling film, the plate equation takes the general form [24] where ξ and ξ iv are the derivatives of order 2 and 4, respectively, of the normal deflection ξ with respect to x 1 , S is the plane stress component in the first direction of the uniformly deformed film, and f s is the normal pressure exerted in the third direction by the (infinitely deep) substrate.
Taking the wrinkling solution [12,24] ξ = ξ 0 cos(qx 1 ), yields The critical wave number, which minimises the above stress, and the corresponding critical stress then take the form (see also [1]) where β = μ f /μ s , as before.We recover here the scaling with respect to β given in [52] based on an LCE plate theory.Equation (67) and the critical values given by ( 70) are valid for both elastic and LCE films.The difference between these two types of films lies in the material constitutive law, and hence, in the formula for the stress S cr in terms of the deformation (see [61] for details).Specifically, in each of the four cases described in the previous section, the critical plane stress S cr and the critical stretch ratio λ cr are related as follows: Cases 1 and 2. The film is made of an isotropic elastic material, hence, and, by (70), Case 3 and 4. The film is made of a monodomain LCE, hence, and, by (70), In each case, the critical stretch value for the plate model approximates well the corresponding value for the 3D models in the asymptotic limit of large stiffness ratio, β → ∞.We conclude that the plate model is applicable when the film is much stiffer than the substrate.

Conclusion
We conducted a linear stability analysis for the onset of wrinkling in nonlinear two-layer systems formed from a hyperelastic film on a LCE substrate, or a LCE film on a hyperelastic substrate.We assumed that the hyperelastic material is described by a neo-Hookean strainenergy function, while the LCE is characterised by the neoclassical strain density.We also assumed that the substrate is infinitely deep and that the nematic director is uniformly aligned either parallel with or orthogonal to the biaxial forces causing wrinkling.We note that bilayer systems where the substrate is of arbitrary thickness can also be analysed by the same method if the original 8 × 8 algebraic system, which we derived, is used instead of its reduced 6 × 6 form corresponding to the case with infinitely deep substrate.Typically, in wrinkling problems where the substrate is about 10 times thicker than the film, depth effects are negligible.Hence, our analysis applies to this situation.
To illustrate the general theory, we specialised to a biaxial stretch resulting from a combination of elastic and natural deformations.In this case, we showed that, if the LCE is at high temperature then cooled, it extends along the director and contracts in the perpendicular direction causing wrinkling parallel to the director, and if it is at low temperature then heated, it contracts along the director and extends in the perpendicular direction causing wrinkling perpendicular to the director.At constant temperature, we found that, if the compressive force acts along the director, it can first cause the director to rotate and align perpendicular to the compressive force, then wrinkling parallel to the rotated director can form.We further considered the case when the film deformation is approximated by a plate theory, and demonstrated the validity of this approximation when the film is much stiffer than the substrate.
From classical wrinkling problems, the linear stability analysis provides an accurate picture of the material behaviour when the film is much stiffer than the substrate, i.e. when the parameter β is large, as the instability is supercritical [3].For smaller values of the stiffness ratio β, we expect that this instability will be replaced by a subcritical bifurcation and, possibly, a creasing instability [20].However, a full weakly nonlinear analysis needs to be carried out to study this regime.Such an analysis will be more meaningful if there is a suitable experimental situation of interest.eigenvectors) with the left Cauchy-Green tensor given by ( 6), and with the Almansi strain tensor given by (7).
The corresponding first Piola-Kirchhoff stress tensor (representing the internal force per unit of undeformed area acting within the deformed solid) is equal to where Cof(A) = (det A) A −T is the cofactor of A. This stress tensor is not symmetric in general.
The associated second Piola-Kirchhoff stress tensor is and this is symmetric and coaxial with the right Cauchy-Green tensor defined by ( 6) and with the Green-Lagrange strain tensor defined by (8).
Next, we derive the stress tensors of a deformed LCE, with the strain-energy function described by (9), in terms of the stresses in the base polymeric network when the nematic director is 'free' to rotate relative to the elastic matrix (see also [61]).In this case, F and n are independent variables, and the Cauchy stress tensor for the LCE is calculated as follows, where T is the elastic Cauchy stress defined by (A.1), J = det F, and the scalar p (nc) (the hydrostatic pressure) represents the Lagrange multiplier for the internal constraint J = 1.
Because the Cauchy stress tensor T (nc) given by (A.4) is not symmetric in general [5,70,96], the following additional condition must hold [5,96], By the principle of material objectivity stating that constitutive equations must be invariant under changes of frame of reference, this is equivalent to (see [5] for details), where T (nc) − T (nc) T /2 represents the skew-symmetric part of the Cauchy stress tensor.The corresponding first Piola-Kirchhoff stress tensor for the LCE is equal to where P is the elastic first Piola-Kirchhoff stress given by (A.2).The associated second Piola-Kirchhoff stress tensor is where S is the elastic second Piola-Kirchhoff stress tensor given by (A.3).
The strain-energy density (5) takes the equivalent form where {α 2 i } i=1,2,3 are the eigenvalues of the Cauchy-Green tensor given by (6).Then, for LCEs, where the nematic director can rotate freely relative to the elastic matrix, the strainenergy function given by ( 9) takes the equivalent form where {λ 2 i } i=1,2,3 and {e i } i=1,2,3 denote the eigenvalues and eigenvectors, respectively, of the tensor

Appendix B. The 3D equilibrium equations
In this appendix, we show that the elastic energy of an ideal LCE is equal to the elastic energy of its polymeric network (see also [61]).We consider a solid LCE body characterised by the strain-energy function defined by ( 9), and occupying a compact domain Ω ⊂ R 3 , such that the interior of the body is an open, bounded, connected set Ω ⊂ R 3 , and its boundary ∂Ω = Ω\Ω is Lipschitz continuous (in particular, we assume that a unit normal vector N exists almost everywhere on ∂Ω).The elastic energy stored by the body is equal to [64, p 205] where p (nc) (det F − det G) enforces the condition that det F = det G, with det G representing the local change of volume due to the 'spontaneous' deformation.
Assuming that, on part of its boundary, Γ N ⊂ ∂Ω, the body is subject to a surface force (traction) f N that is also derivable from a potential function, we are interested in the minimisation of total potential energy subject to det F = det G, over all displacement fields u = u(X) = x − X, with gradient tensor ∇u = F − I, such that the kinematic constraint (2) holds.Following a standard procedure of computing variations δE in E for infinitesimal variations δu in u, the following first variation in E t is obtained (see [64, p 310-312]), where P (nc) is the first Piola-Kirchhoff stress tensor (given explicitly by (A.7) in appendix A).
The principle of stationary potential energy then states that, of all admissible displacements fields that the body can assume, those for which the total potential energy assumes a stationary value, such that δE t = 0, correspond to static equilibrium state described by together with the traction boundary condition causing the deformation, P (nc) N = f N , where N is the outward unit normal vector to the external surface Γ N .
In order to simplify the mechanical analysis, it is convenient to rewrite the above equations equivalently in terms of the elastic stresses for the base polymeric network.As the elastic deformation is applied directly to the reference state, using the relations between the stress tensors for the LCE, we obtain [61]

Appendix C. Analysis at q → ∞
The special limiting case where the wave number goes to infinity is interesting as it corresponds to the Biot instability and arises when the stiffness ratio β = μ f /μ s is sufficiently small.In practice, this instability is not observed because creasing, folding, or a subcritical wrinkling instability occur before that limit.Nevertheless, it is an important reference point to organise the bifurcation diagram, and we include this here for completeness.

Figure 1 .
Figure 1.Wrinkling parallel to the nematic director n when a > 1 and the system is at high temperature then cooled, so it extends along the director and contracts in the perpendicular direction[1,2].

Figure 2 .
Figure 2.Wrinkling perpendicular to the nematic director n when a < 1 and the system is at low temperature then heated, so it contracts along the director and extends in the perpendicular direction[1,2].

Figure 3 .
Figure 3. Schematic of the composite deformation of a nematic solid.
and by Γ C = ∂Ω f ∩ ∂Ω s ⊂ ∂Ω the open subset representing the interface between the two layers.The external boundary Γ = ∂Ω\Γ C is partitioned as Γ = Γ D ∪ Γ N , such that on Γ D , essential (displacement) boundary conditions are prescribed, while on Γ N , natural (traction) boundary conditions are imposed.