Wrinkling composite sheets.

We examine the buckling shape and critical compression of conﬁned inhomogeneous composite sheets lying on a liquid foundation. The buckling modes are controlled by the bending stiﬀness of the sheet, the density of the substrate, and the size and the spatially dependent elastic coeﬃcients of the sheet. We solve the (linearized) F¨oppl–von K´arm´an equations describing the mechanical equilibrium of a sheet when its bending stiﬀness varies parallel to the direction of conﬁnement. The case of a homogeneous bending stiﬀness exhibits a degeneracy of wrinkled states for certain sizes of the conﬁned sheet. This degeneracy disappears for spatially dependent elastic coeﬃcients. Medium length sheets buckle similarly to their homogeneous counterparts, whereas the wrinkled states in large length sheets localize the bending energy towards the soft regions of the sheet.


I. INTRODUCTION
From surfactant monolayers to tectonic plates, the deformation of elastic sheets underlies a vast landscape of problems in science and engineering. Not only do their patterns evoke a great aesthetic appeal, but they underlie both function and form in the world around us [1][2][3][4].
The buckling of a homogeneous sheet on top of a supporting foundation has been the object of soft matter studies for decades [see e.g., 5-7, and Refs. therein]. In a free-standing confined sheet, the buckling compressive load (which we refer to as "compression" throughout the paper) and wavelength are set by the size 2a * of the confined sheet. However, in the presence of a supporting foundation, an intrinsic length scale ℓ * is set by the mismatch between the elasticity of the sheet and that of the foundation. When a resting sheet with clamped ends is confined longitudinally, a uniaxial compression builds up and wrinkles of the surface emerge [8]. Such wrinkle patterns have been observed in, among many other settings, thin polymer sheets resting on the surface of water [9], bacterial biofilms [10], granular rafts [11], the mineral veins of rocks [3], strained epitaxial films [12], and in glaciology [13,14]. From a theoretical standpoint, wrinkles of confined elastic sheets are predicted by a simple set of geometric rules [15].
In many settings the deforming material is treated as a homogeneous solid, despite a potentially important intrinsic composite structure. Indeed, whether or not the composite structure can be treated as effectively homogeneous requires appropriate quantitative analysis. When the host material is stiff, such as ice, and the inclusions are soft, the theory of Eshelby [16] shows that the effective Young's modulus, and hence the bending stiffness is reduced relative to the case of the stiff host material alone. When the host solid is soft and the inclusions are uniformly distributed [17], two possibilities exist -composite softening or stiffening-depending on the size of inclusions relative to the elastocapillary length, l * ≡ γ * /E * 0 , where γ * is the surface tension of the inclusion host interface and E * 0 is the Young's modulus of the host material. However, by controlling the distribution of inclusions we can prepare inhomogeneous composites, having for example a gradient of bending stiffness, and thus manipulate the failure properties of such elastic composites.
Here we examine the deformation of composite elastic sheets floating on a liquid foundation. We treat the case when there is a gradient of bending stiffness parallel to the direction of confinement. This gradient is envisaged to be created by controlling the distribution of liquid inclusions.
First, we present key results for the critical compressive load and wrinkle patterns in homogeneous floating sheets. The interplay between the intrinsic length scale, ℓ * , and the size of the confined sheet, 2a * , determines the wrinkled state that has the lowest compression. The observed wrinkles are then either symmetric or antisymmetric about the center of the sheet. Importantly, at certain sheet sizes the wrinkled states are degenerate: instead of having a well-defined symmetry, they can be the superposition of any symmetric and antisymmetric states. These degeneracies are a striking feature that enables the manipulation of the wrinkles, and hence a promising range of practical applications, such as the propagation of mechanical signals in soft media, which is an area of current interest [e.g., 18,19].
While the homogeneous case provides important insights into the wrinkle pattern, our main interest is in the study of inhomogeneous sheets. Taking advantage of the effective medium behavior of composites [16,17,20], we examine the wrinkles in composite sheets.
We impose a spatial distribution of liquid inclusions that translates into spatially varying elastic moduli (Young's modulus and Poisson ratio) and hence a spatially varying bending stiffness. The inhomogeneous stiff- ness is responsible for breaking the symmetry of the wrinkle patterns and eliminating the degeneracies. For sheets large compared to the intrinsic length scale, this gradient enhances the localization of wrinkles. This localization happens at the onset of wrinkling, and hence it does not require any non-linear behaviour. Inhomogeneous sheets are technologically appealing for they offer a new means to control fracturing via the localization of the wrinkles.

II. MECHANICS OF A FLOATING COMPOSITE SHEET
We start with a presentation of the governing equations for a composite thin sheet floating on a liquid, including the effects of in-plane forces. For more detail the reader is referred to the book by Mansfield [4].

A. Problem formulation
Let x * , y * , z * (all dimensional quantities are denoted with a star, * , superscript) be the Cartesian coordinates in the horizontal direction, into the page, and in the vertical direction, respectively. We consider a sheet with thickness h * and longitudinal undeformed length L * , such that h * ≪ L * . The sheet (of density ρ * s ) rests on a liquid of density ρ * , and in the absence of all forces apart from gravity, it floats with its mid-line at a height z * = h * (1/2 − ρ * s /ρ * ) above the free surface. We measure all vertical displacements relative to this equilibrium level [21], as in Fig. 1 where the origin of the vertical axis is set at the liquid level. The sheet is confined lengthwise to a distance d * , and has clamped edges located at x * = ±a * in the longitudinal direction, with a * ≡ (L * − d * )/2 (see Fig. 1).
The geometric incompatibility between the sheet's undeformed length and the confinement length yields a mechanical instability producing a vertical displacement w * (x * , y * ). This displacement also takes into account the restoring pressure p * = −ρ * g * w * that the liquid foundation exerts on the sheet. An in-plane compression τ * (with dimensions of force per unit length) is calculated as an emergent property, rather than being imposed.
The static equilibrium of a confined floating sheet is determined by the balance between elastic forces and the hydrostatic pressure. The former arise directly from the in-plane forces at the edges due to confinement, and the latter is purely normal (i.e. the in-plane components are zero), and hence the two in-plane equilibrium equations are reduced to one by introducing a force function χ * . Thus, the internal in-plane forces per unit length are obtained by double differentiation of χ * [4]. Following the standard derivation of the plate equation, see [4], but accounting for the possibility that ν = ν(x * ), we find that normal displacements of the sheet satisfy: where B * (x * ) is the bending stiffness (or flexural rigidity) of the sheet and the von Kármán operator is Importantly, the bending stiffness, which has the dimension of energy, is a combination of the thickness h * and the elastic properties of the material, and hence it is a function of x * for composite materials whose elastic properties are inhomogeneous; viz. the Young's modulus, E * (x * ), and the Poisson ratio, ν(x * ). Such inhomogeneous sheets can be made by embedding one material in the other to control the resulting composite structure. Despite a varying composition, we neglect any variation of the density ρ * s of the sheet. This is justified for stiff composites in the dilute regime, and for soft composites, because the matrix and the inclusions have approximately the same density. A varying bending stiffness would also correspond to a single component material sheet with varying thickness, but we do not treat this case here (because it would complicate the free-floating equilibrium). For a midplane displacement with components (u * , v * , w * ), the sheet's in-plane strains, ǫ ij , are given by where, for example, (w * ) ,x * denotes the partial derivative of w * with respect to x * . Note that the displacements u * and v * may be eliminated from these relationships by using the condition of compatibility [see pg. 13 of Ref. 4].
Relating the strains to the in-plane forces, and hence to the derivatives of the force function χ * , one finds which gives the equilibrium in the plane of the sheet. The in-plane forces per unit length are then derived from χ * by double differentiation [4].
In the absence of inclusions, namely for an homogeneous elastic sheet, the elastic coefficients are constant (Young's modulus, E * 0 , and Poisson ratio, ν 0 ; we denote the elastic constants of the homogeneous sheet with subscript "0") and hence so too is the bending stiffness; in particular, the baseline bending stiffness using Eq. (3). The intrinsic length scale of the displacements, or wrinkles, ℓ * ≡ (B * 0 /ρ * g * ) 1/4 , plays a central role in the problem. Thus, rescaling lengths with ℓ * , and the force function χ * with the bending stiffness, we make equations (1) and (5) dimensionless. Dimensionless quantities are unstarred so that x = x * /ℓ * , χ = χ * /B * 0 , etc., and hence equations (1) and (5) become respectively, where the dimensionless elastic modulus is E(x) ≡ E * (x * )/E * 0 , and the composite bending stiffness is . These equations are commonly referred to as the Föppl-von Kármán equations [4]. Their unknowns are the deflection w(x, y) and the force function χ(x, y). Importantly, the terms on the right hand side of Eqs. (7) are non-linear: [χ, w] couples the curvature with the in-plane compression in the balance of normal forces, and [w, w]/2 is the Gauss curvature. Note that in Eq. (7b) the dimensionless ratio is a measure of the relative ease of bending to stretching of the sheet, or the stretching-stiffness. Because B * 0 ∼ E * 0 h * 3 so that S ∼ [E * 0 /(ρgh * )] 1/2 , for thin sheets we expect that S ≫ 1, and hence deformation can be accommodated more easily by bending than stretching. For example, a 10 cm thick sheet of ice, for which E * 0 ∼ GPa and ν 0 ≈ 0.3, has S ∼ 1000 [2]. A soft material, such as PDMS, for which E * 0 ∼ 100 kPa and ν 0 ≈ 0.5, has S ∼ 100 for h * on the order of cm. However, thin sheets do not always bend more easily than they stretch. Provided that S is finite, sheets may accommodate some of the imposed deformation by compressing (negative stretching), particularly when the confinement is small and the supported sheet resists buckling. We therefore refer to S as the 'inextensibility' of the sheet. Before returning to this idea, we introduce buckling in two dimensions where the equations of the mechanical equilibrium (7) simplify.

B. Two-dimensional buckling of a sheet
In the two-dimensional buckling problem shown schematically in figure 1, there are no variations in the y direction (i.e. into the page). Therefore, the vertical displacement of the sheet, w, is independent of y and equations (7a)-(7b) simplify to the following system; and . However, χ ,xy is just the traction exerted on the sheet in the y direction. This traction is zero for compression purely in the x direction, so we have f ′ (x) = 0, and hence f (x) is a constant. Since the compressive load at the boundary x = a (a ≡ a * /ℓ * ) is in the x direction, and has magnitude τ ≡ τ * (ℓ * ) 2 /B * 0 , we must have and so that f (x) = −τ . Thus Eq. (9a) becomes which is Euler's linearized elastica equation [see e.g. §20 of Ref. 22] with a lateral load due to the hydrostatic pressure in the liquid foundation and varying elastic properties along the axis of confinement. This is analogous to the small-deflection equilibrium of a compressed inhomogeneous column [23]. The boundaries of the thin sheet at x = ±a are clamped so that w ,x (±a) = w(±a) = 0. Recall that the bending stiffness is a function of the longitudinal coordinate x through the elastic coefficients E(x) and ν(x). These functions account for the variation of the liquid volume fraction φ; we treat composites of isotropic elastic host solids into which liquid inclusions are embedded, as shown in Fig. 2, using a linear distribution as where φ(x = −a) ≡ φ −a . The elastic response of a composite is characterized by the elastocapillary length l * introduced in §I. When inclusions of radius R * are smaller (larger) than l * , the composite is difficult (easy) to deform because surface tension can maintain (cannot maintain) their sphericity. For example, for a typical liquid inclusion with γ * = O(0.01N m −1 ) [17], and soft materials with a range of E * = O(1−100kPa), l * ranges from 1 to 100µm. Thus, surface tension effects are important for micron-sized inclusions. A summary of different theoretical models that predict the bulk mechanical properties of a composite is included in Appendix A.

III. BUCKLING OF A HOMOGENEOUS SHEET
The force balance of an elastic sheet with homogeneous properties, namely B(x) = 1 in Eq. (11), is given by The boundary conditions corresponding to a clamped sheet are: Again, the subscript "0" denotes those physical quantities of the homogeneous case.

A. Buckling profiles
Eqs. (13) and (14) define an eigenvalue problem whose detailed solution is given in Appendix B. As expected from the reflection symmetry about x = 0, there are two distinct types of solutions with opposite parity: even functions of the form w profiles. For both cases the wavenumber k satisfies Therefore, there are two possible wavenumbers, k ± , given by that are plotted in Fig. 3. Importantly, the domain of k ± (τ 0 ) (such that k ± (τ 0 ) ∈ R) is bounded from below, where τ 0 = 2 gives the value where the two branches of the solutions meet, k + (τ 0 ) = k − (τ 0 ). We show in Appendix C that there are no solutions of the eigenvalue problem (Eqs. 13,14) for complex k ± . In general w 0 (x) will contain both of the wavenumbers given in (16). The condition of zero vertical displacement at x = ±a is satisfied by the symmetric and antisymmetric combinations viz., and w (a) where A s , A a are as yet to be determined constants. The remaining boundary condition that (w 0 ) ,x (±a) = 0 leads to the relations between k + and k − ; with which are associated the symmetric (17a) and antisymmetric (17b) modes respectively. Since k ± = k ± (τ 0 ), the solutions of Eq. (18a) determine the compression τ 0 (x) from Eq. (17b). For a given value of a, each relation (18a) and (18b) has an infinite number of solutions, the smallest of which is τ 0 = 2 (see Fig. 3). When τ 0 = 2, k + = k − , Eqs. (17) give w 0 (x) ≡ 0. Therefore, each value of τ 0 > 2 that solves either of Eqs. (18) corresponds to a different mode of buckling in the sheet. We shall consider only the lowest mode of buckling, which corresponds to the smallest value of τ 0 > 2 that solves either Eq. (18a) or Eq. (18b).
The non-zero lower bound for the compression reveals compression prior to buckling as discussed in §II A. We consider a sheet that is confined to a domain of size 2a slightly shorter than its undeformed length (L) such that its associated compression τ a is smaller than the lowest value of τ 0 (a). The subscript "a" distinguished the applied compression needed to confine the sheet to a size 2a from the eigenvalues of the mechanical equilibrium from Eq. (13). Since τ a < τ 0 (a) the sheet is stable in its flat configuration, i.e. w(x) = 0, and the strain along the direction of confinement, Eq. (4), is where we have used the fact that the compression is purely in the x direction, so that χ ,2y = −τ a and χ ,2x = 0. The integral of Eq. (19) between x = −a and a then gives L. Below that size, the sheet buckles and it tends to recover its undeformed arclength, bending rather than compressing.

B. Amplitudes
Having determined τ 0 by solving Eqs. (18), the solution for the shape of the buckled sheet is given by Eqs. (17) up to the multiplicative constants A s and A a . Unlike the infinitesimal amplitude of the buckling shape of an incompressible column, these amplitudes are finite at the onset of buckling because the sheet resists buckling by compressing when the confinement is small (for S ≫ 1). The value of the amplitude is determined by the undeformed length of the sheet, L, in its undeformed state. Because S ≫ 1, we neglect stretching of the sheet in the buckled state, and if we assume small deformations, the contour length of the deformed sheet must be equal to its undeformed length, viz., (20) Integrating Eq. (4) for the strain ǫ x * x * from x = −a and a, and using symmetry about x = 0, results in where we wrote the strain in terms of the force function and substituted a compression τ 0 in the x-direction, i.e. χ ,2y = −τ 0 and χ ,2x = 0. Using w 0 (x) from either Eq. (17a) or Eq. (17b) we find which can be substituted into Eq. (21) to give A s and A a . Since k ± and τ 0 are determined by the numerical solution of Eq. (18a) (or Eq. 18b), they depend only on the value of a. Therefore, for given values of a and L, the shape of the sheet can be completely determined numerically, including the amplitude, up to a sign. For convenience, we letã which denotes the horizontal projection of the confined sheet in terms of the number of cycles (or wavelengths) for a sinusoidal function with unit wavenumber.
We classify the sheets according to their size: Short length sheets are shorter than the instrinsic length scale of the wrinkles, ℓ * , (ã < 1), medium length sheets are a few length scales in size [ã = O(1, 10)], and large length sheets are many length scales in size (ã ≫ 1). This classification is relevant to the choice of the displacement parallel to the direction of confinement. Wrinkled surfaces in large length sheets are stable against further confinement up to a displacement d c = (2π) 2 /L [24]. This prediction is supported by experimental results of large length sheets, for which d c ≈ 0.3 [9]. Given the L-dependent critical displacement, we adopt d = d * /ℓ * = 0.2 (i.e. d * = ℓ * /5) for large length sheets, and d = L/50 for medium and short length sheets. Both cases meet the condition of stable wrinkled surfaces, d < d c .

C. Effects of confinement
Eq. (16) and the relations (18) show that the compression τ 0 is a function of the control parameterã. Thus, for different values ofã, we compute the two smallest numerical roots of each relation (Eqs. 18 have a countable infinity of solutions) and we plot them in Fig. 4 (a) for short and medium length sheets, whereas for large length sheets in (b), we only plot the two lowest compressions. We recall that a compression at a given value ofã is associated with a wrinkled state which is either symmetric or antisymmetric, denoted by the index 's' or 'a' respectively. The four curves τ 0 = τ 0 (ã) in Fig. 4(a) are arranged in pairs τ (s,j) 0 (ã), τ (a,j) 0 (ã) , labeled by index j; j = 1 is the lower and j = 2 the upper.
For short length sheets (ã < 1), in every pair, the symmetric state has a lower compression than the antisymmetric one. Furthermore, the compressions are inversely proportional toã 2 . This power-law behavior can be understood by dimensional analysis; as the size of the confined sheet decreases it becomes the smallest length scale in the problem and hence τ * We confirm the inverse square law in Appendix E. However, a fundamental assumption for the Föppl-von Kármán equations (7) is that the size of the sheet is larger than its thickness, and hence we will not explore the limit in which this condition is violated. In the large-length limit, plotted in Fig. 4 (b), the functions τ 0 (ã) have an almost periodic component (see Appendix E) and they asymptotically approach the minimum τ 0 = 2.
For every pair j, τ (s,j) 0 and τ (a,j) 0 are out of phase and there are specific values ofã at which they cross. Hence, the wrinkled state with the lowest compression is alternatively symmetric or antisymmetric. When the modes cross, the two relations (18) are satisfied by the same values of τ 0 andã. Using this condition, we find two families of crossings (see Appendix D). The crossing points of what we call type I are located at with l, j ∈ N ⋆ ; index l labels the crossings within each pair j -starting at l = 1 for the crossing at the smallestã, which is of type I, and increasing with the size of the confined sheet. Similarly, for what we call type II crossings we haveã = l(l + j) ≡ã These indices also explicitly determine the corresponding compressions [25]. Note that the crossings of type I are determined by the Indeed, after Eq. (17a), where we used the property (18a) of symmetric modes in the second equality. Imposing the condition (26), the vanishing of the bracket yields the property (18b) of antisymmetric modes. Analogously, we show that the crossings of type II are determined by the condition At crossings of type I (II), the wrinkled state with the lowest compression changes from symmetric to antisymmetric (from antisymmetric to symmetric) asã increases. A symmetric (antisymmetric) mode has an even (odd) number of nodes. Hence, a new node appears at every crossing when the size of the confined sheet increases. In that sense, the index l accounts for the number of nodes. These symmetric/antisymmetric transitions, and the subsequent appearance of new nodes, are intuitive if one considers the case of a very large sheet (a * ≫ ℓ * ). In that case the wrinkles are perfectly sinusoidal, and hence they define a unique wavelength λ * = 2πℓ * independent of the size of the system. Depending on the number of half wavelengths, λ * /2, that fit in the size of the confinement, 2a * , the number of nodes (including the two nodes at the clamped ends) is given by ⌊2a * /(λ * /2)⌋ + 1, where ⌊.⌋ denotes the floor function. When this number is even/odd, the wrinkles are symmetric/antisymmetric about the center of the sheet. Thus, the number of nodes changes by ±1, and the symmetry does change too, at values of a * such that 2a * /(λ * /2) = 2a * /(πℓ * ) is an integer.
For the finite sheets whose length is not very large compared to ℓ * , the symmetry changes occur at values of a * (cf. Eqs. 24 and 25) which are nearly, but not exactly, multiples of πℓ * /2. This behavior is concomitant with the almost periodic periodic component of the functions τ 0 (ã). The wrinkle pattern has then two wavelengths determined by the ratio a * /ℓ * . Fig.3 shows how the two wavelengths corresponding to the wavenumbers k ± (λ ± = 2π/k ± ) depart from 2π for shorter lengths (i.e. larger values of the compression τ 0 ).
The derivative τ ′ 0 ≡ dτ 0 /dã ≤ 0; see e.g., τ (s,1) 0 [ã], plotted near the crossingã (I) 11,1 in the inset of Fig. 4 (b), and its numerical derivative in the other inset. This indicates that the wrinkled state of a sheet (confined initially to a size 2a) is stable to small changes in the externally applied compression τ a = τ 0 (a) + δτ a where δτ a can be positive or negative, and is balanced by a small change in the confining size δa = δτ a /τ ′ 0 (a), such that τ a = τ 0 (a + δa). Intuitively, a negative τ ′ 0 corresponds to the increasingly strong response (τ 0 larger) of the sheet to shorter confining sizes 2a. This type of behavior is common in elastic solids, which usually display monotonic force-displacement curves. However, force-displacement curves whose first derivative changes sign are observed in buckling experiments [e.g., 26]. Those cases involve more subtle physics including hysteresis.
In both crossings of type I and II, two solutions coexist and any linear combination of them that we can generate by introducing a phase θ, will also be a solution. We plot in Fig. 5 the coexisting buckling profiles of a confined sheet of sizeã =ã (I) 2,1 , given by the linear combination in Eq. (29) with θ ∈ [0, π]. The nodes (depicted by the solid symbols in Fig. 5) are continuously displaced along the sheet, and hence we can gradually shift the wrinkles and create/anihilate one node. Ideally, this transformation is done without any external energy supply to the system and the compression of the sheet remains constant. This features the possibility of using buckled floating sheets as supports to transmit signals, exploiting multistability similarly to the elastomeric bistable beam elements in the work by Raney et al. [18] and the bistable shells by Vasios et al. [19].

IV. BUCKLING OF INHOMOGENEOUS SHEETS
We now imagine adding liquid inclusions to our soft host and solve the problem numerically. We use the numerical methods from Chebfun [27] to simulate the mechanical equilibrium Eq. (11). We consider a sheet with dimensionless stretching stiffness S ≈ 300, that corresponds to a PDMS matrix (E * 0 = 10 5 Pa, ν = 0.5, ρ * PDMS = 0.95 · 10 3 kg · m −3 ) of thickness h * = 1 mm. The sheet is confined and we use a linear distribution (Eq. 12), parallel to the direction of confinement, of stiffening (γ′ = 10, with γ′ ≡ l * /R * as per the Appendix A; R * denotes the size of the inclusions), incompressible liquid inclusions. The scaled bending stiffness B(x) equals the dimensionless effective Young's modulus for stiffened soft composites given by Style et al. [17,28]: where φ is a function of x given by Eq. (12). The control parameters here are the liquid volume fraction at x = −a, φ −a , and the size of the confined sheet,ã.

A. Critical compression of a soft composite sheet
For homogeneous sheets, we saw in §III how the symmetry about x = 0 enforces solutions w(x) with a well defined symmetry (see e.g. Fig. 1) and the degeneracy of the associated compressions (see Fig. 4). The gradient of bending stiffness created by varying the volume fraction of liquid inclusions (cf. Eq. 12) breaks the symmetry of the solutions of Eq. (11), and the corresponding eigenvalues do not cross [29,30]. Thus, the compressions corresponding to the different solutions Eq. (11) are ordered from smaller to larger for any finiteã, and the intersecting compressions in confined homogeneous sheets (whose size is given by Eqs. 24 and 25) split.
We show in Fig. 6 the splitting between the two lowest compressions of the inhomogeneous sheet, denoted τ (1) (ã) and τ (2) (ã), with τ (1) (ã) < τ (2) (ã). Let  (1) be the separation between neighbouring compressive loads. For a homogeneous sheet, it vanishes at certain values ofã given by the crossing points from Eqs. (24) and (25). We see in Fig. 6 that ∆τ (2−1) never vanishes in the inhomogeneous case, but still has a local minima around every crossing point. We next examine the growth of ∆τ (2−1) with φ −a . We fix φ −a and then we find numerically the local minima of ∆τ (2−1) with respect toã, which we denote ∆τ (2−1) min . The values ofã that minimize ∆τ (2−1) start at the crossing points when φ −a = 0, and increase with φ −a . The first eight results ∆τ (2−1) min (corresponding to the first eight crossings in the homogeneous case) are plotted in Fig. 7 for different values of φ −a /a, where a is the corresponding sheet size that minimizes ∆τ (2−1) . All curves ∆τ (2−1) min collapse into one, therefore we conclude that ∆τ (2−1) min grows linearly with φ −a /a, which is twice the amplitude of the gradient from Eq. (12).

B. Buckling profiles of medium length sheets
For medium length sheets [ã = O(1, 10)], we examine the buckling profiles corresponding to the lowest compression, τ (1) . The contour plot of w(x) in Fig. 8 shows how the buckling profiles are no longer symmetric and the largest deflections (i.e. those brighter/darker areas) tend to be located towards the softer side of the sheet (x > 0). This intuitive feature is also observed in the buckling profiles of an inhomogeneous column [23]. Despite the lack of symmetry, the buckling profiles of a medium length inhomogeneous sheet are similar to the fundamental symmetric/antisymmetric modes characteristic of the homogeneous sheets. We project the numerical solutions for the inhomogeneous sheet, w(x), onto the modes of the homogeneous sheet (Eqs. 17), i.e. the basis of eigenfunctions of the eigenvalue problem given by Eq. (13) [31]. More concretely, we expand the solution of Eq. (11) as the infinite linear combination where we use the same index j to denote each pair of solutions -symmetric (s) and antisymmetric (a)-as in the homogeneous sheet. We redefine the eigenfunctions (17a) and (17b) as and w The coefficients a 1 are plotted in Fig. 9. The first pair of modes -either the symmetric or the antisymmetric -are dominant in the linear combination that describes w(x). This dominant term changes in the vicinity of the crossing points (depicted by vertical dashed lines in Fig. 9). However, for larger and more inhomogeneous sheets the sum a This indicates a reduction in the projection of w(x) onto the first pair of modes. Higher order modes are necessary to account for the buckling profile of longer and more inhomogeneous sheets. Thus the difference between the wrinkles in homogeneous and inhomogeneous sheets increases. A rationale behind this result is found in the compression graph [ Fig. 4 (a)] of a homogeneous sheet: all buckling modes converge towards τ = 2 in the limit of large sheets (ã ≫ 1), and so they contribute similarly to the buckling profile. As τ → 2 the wavenumbers defining the different modes in the homogeneous sheet tend to accumulate at unity, see Fig. 3.  24) and (25).

C. Large length sheets: Localization and fracture in hard composites
The results of Fig. 10 highlight a fundamental difference between inhomogeneous sheets of medium and large lengths. While the buckling profile of medium length sheets is described in terms of the first pair of modes, many modes combine to produce the wrinkles in long sheets. Localization of the wrinkles originates in this superposition of modes which accumulate at the softer end of the sheet, as seen in Fig. 11 for the long inhomogeneous sheet (S ≈ 300, d * = L * /5) with stiffening inclusions (γ′ = 10 in Eqs. 12 and 30). Localization is also observed in the non-linear regime when a floating sheet is confined beyond the critical compression and wrinkles attenuate, except near the center of the sheet, until a fold emerges [9,24,33,34].
For brittle materials such as ice, the localization of bending is crucial in characterizing their failure. Because ice is a hard material, we model the composite structure of water inclusions with the effective Young's modulus and the Poisson's ratio of Eshelby [16] as (see Appendix A). A model for fracture in ice floes is given in Vella and Wettlaufer [2]; when a crack is initiated, the stresses exceed a critical value σ * m [35]. For elastic sheets, stresses vary linearly through the thickness of the sheet, and so the maximum stress is achieved at the surface of the sheet. This stress is related to the maximum bending moment per unit length acting on an element of the sheet [see pg. 5 of Ref. 4]: where the maximum bending moment per unit length is w * ,2x * (x * ) is the curvature of the sheet in the small slope approximation (which is implicit in the derivation of the Föppl-von Kármán equations (7a) and (7b), [4]). This maximum represents a trade-off between the bending stiffness, which peaks at the stiffer end of the sheet, and the curvature of the buckling profile, which is larger towards the softer side of the sheet.
We rescale the bending moment per unit length by B * 0 (ℓ * ) −1 , and the stress by B * 0 (h * ) −2 (ℓ * ) −1 . Now, considering a sheet with constant thickness h * and varying Young's modulus, the dimensionless failure criterion is We plot in Fig. 12 (a) the buckling profiles of a thin floating sheet with S ∼ 10 4 and a decreasing volume fraction of softening inclusions (Eq. 12 and elastic moduli given by Eqs. 35,36) along the direction of confinement (ã = 10, d = L/50). The vertical axis in the contour plot denotes an increasing volume fraction at x = −a, φ −a . The choice for the dimensionless stretching-stiffness models the behaviour of a fresh ice sheet with ρ * ice = 0.9 kg · m −3 , E * 0 = 1 GPa, h * = 1 mm and ν 0 = 0.33 [36,37], floating on water. In addition to the wrinkle localization on the softer end of the sheet, which is more prominent at larger values of the volume fraction φ −a , there is a deviation between the position of the largest deflection (green dots) and the coordinate of |M max | (red dots). To wit, the sheet's maximum bending moment is not at the largest wrinkle, but at wrinkles on the stiffer side.
The maximum bending moment and the largest deflection in Fig. 12 (a) move stepwise towards the soft end as φ −a increases, and the separation between their positions changes with φ −a . Now, because the largest curvature occurs at the largest deflection, and |M max | is proportional to the curvature (see Eq. 38) then |M max | is a discontinuous function of φ −a . In Fig. 12 (b), the discontinuities correspond to when the maximum bending moment jumps to the left and hence approaches the location of the largest deflection. Finally, we note that the maximum bending moment, plotted in Fig. 12 (b), increases in inhomogeneous sheets.

V. CONCLUSIONS
We have examined the two-dimensional buckling and wrinkle patterns in floating homogeneous and inhomogeneous thin elastic sheets. With two control parameters, the size of the confined sheet and the gradient of bending stiffness, we quantified the wrinkled states using the Föppl-von Kármán theory of thin sheets. A central test of the results is to vary the bending stiffness by varying the volume fraction of inclusions in the host solid.
In homogeneous sheets, the only control parameter is the size; the buckling profile is determined by the lowest compressive load. Thus, the mode that will be observed in a buckling experiment corresponds to the lowest compression. However, for some particular sizes of the confined sheet, the same compression is associated with two different wrinkling modes that can be continuously displaced along the sheet (see e.g. Fig. 5).
In contrast, this degeneracy is not observed in inhomogeneous sheets. Indeed, the otherwise crossing compressions of the homogeneous case split when a gradient of stiffness is applied parallel to the direction of confinement. The size of this splitting grows linearly with the magnitude of gradient of the volume fraction at all the crossing points. Importantly, the wrinkled states of confined inhomogeneous sheets sensitively depend on their size. While medium length sheets buckle very much like their homogeneous counterparts, the wrinkled states in large length sheets are a superposition of many modes. This feature of large length sheets allows for the localization of the bending energy, which is crucial in establishing a failure criterion, with particular relevance in glaciology.
Finally, the results presented here for floating sheets are also relevant for sheets on a linear soft elastic foundation [see e.g. Ref. 38]. A more complex behavior is expected in the more general case of a linear elastic foundation, whose response is expected to be geometrically nonlinear, in which localization of buckling occurs [see e.g., Ref. 39,40]. In addition to the range of applications of interest, from soft composites of biological relevance [41] to hard composites of engineering or geophysical importance, a thorough mathematical analysis, rather than the numerical study given here, may provide additional insights. The foundational Eshelby theory of solid composites [16] describes the elastic behavior of rigid composites with a dilute dispersion of noninteracting incompressible inclusions. Stiff-matrix materials such as ice, glass, ceramics and steel have E * = O(GPa) and ν ∼ 0.3, and thus have subnanometric elastocapillary length. Therefore, for typical inclusion sizes the effect of surface tension is negligible and we can use Eshelby theory [16] to compute the effective elastic moduli of compression, κ * , and rigidity, µ * , which are We denote the host matrix's elastic constants with subscript "0", those corresponding to the inclusions with subscript "1", and symbols with no subscript denote the solid composite. The volume fraction of inclusions is φ.
The incompressible liquid inclusions have zero shear modulus, µ * 1 = 0, and infinite bulk modulus, κ * 1 = ∞ (due to incompressibility), so that the Young's modulus and Poisson's ratio of the composite are where we have used Eqs. (40) and expanded to first order in φ. For a stiff-matrix composite with ν 0 = 0.3 with liquid inclusions, the composite Young's modulus (Poisson's ratio) is less than (greater than) the host matrix; E * < E * 0 and ν > ν 0 . For soft composites, micron sized inclusions create non-negligible interfacial stresses with an effective Young's modulus given by where γ′ ≡ l * /R * captures the size regime where surface tension operates [17,28]. Eq. (30) assumes that the inclusion concentration is dilute and hence we refer to it as the dilute theory (DT). However, we note that this approach is quantitatively accurate up to φ ≈ 0.2 when γ′ > 2/3, which is in the stiffening regime where E * (φ, γ′ > 2/3)/E * 0 > 1 [42]. The constituents of the composite are incompressible and hence ν = 1/2 throughout.

APPENDIX B: THE EIGENVALUE-VALUE PROBLEM FOR A HOMOGENEOUS SHEET
Here we solve the linear fourth order differential equation for the vertical displacement w 0 , where the eigenvalue τ 0 is determined by the boundary conditions w 0 (±a) = 0 and dw 0 dx ±a = 0.
Since equation (45) has constant coefficients, we seek solutions in the form w 0 (x) ∝ e ikx , yielding and the solutions for k 2 are Thus for the wavenumber k to be real, we need τ 0 ≥ 2. The eigenfunctions w 0 (x) are linear combinations of e ±ik+x and e ±ik−x , but we are interested in real solutions. Thus, we write the general solution of (45) as where the k ± are the positive roots of equation (48) and the C i (i = 1, 2, 3, 4) are real constants. The problem is linear and has homogeneous boundary conditions, so those constants can only be determined up to common factor. Enforcing (46), we obtain The vanishing determinant of the system (50) can be factorized as cos(k − a) k + sin(k + a) k − sin(k − a) sin(k + a) sin(k − a) k + cos(k + a) k − cos(k − a) = 0.
Hence, for a given a we have two possible relations between k ± and a: where 'or' means that either (52) or (53) is fulfilled, or both. Note that when only (52) is fulfilled, we must set C 3 = C 4 = 0 for the last two equations in (50) to be satisfied. Thus, in such a case the eigenfunction is even (see Eq. 49). Similarly, when only (53) is fulfilled, we must set C 1 = C 2 = 0 and the eigenfunction is odd. In Appendix D we ask for which values of a both (52) and (53) are satisfied.
We are interested in the dependence of the eigenvalue τ 0 on the size of the confined shseet a. On the one hand, we can rewrite the characteristic equation (47) as where k can be either k + or k − , noting the general property of the product of roots of a quadratic equation; k 2 + k 2 − = 1. Since we defined k ± to be positive, we have It is clear from equation (54) that τ 0 ≥ 2, as required, and the minimum is reached when k + = k − = 1. On the other hand, using Eq. (55) the relations (52) and (53) take the form and respectively. Thus for any a, k = 1 is a solution of (56) and (57). Moreover, if k is the solution of (56) or (57), so too is 1/k. Equations (56) and (57) actually have a countable infinity of pairs of solutions (k, 1/k), for each a. We call the solutions of Eq. 56 (Eq. 57) symmetric (antisymmetric) modes, because the associated eigenfunctions are even (odd).
Our main focus is the mode with the smallest eigenvalue. Therefore, for a fixed a we independently solve Eqs. (56) and (57), for each of which we seek the pair of solutions yielding the smallest value of τ 0 . In both cases this pair is that closest to (1, 1), but we do not know a priori whether the symmetric or antisymmetric mode has the lowest eigenvalue, a feature that depends on the value of a. We find that the different modes are arranged in pairs, and that each pair contains a symmetric mode and an antisymmetric mode. Moreover, these pairs are well separated from each other (see Fig. 4).

APPENDIX C: THE BUCKLING WAVENUMBER FOR HOMOGENEOUS SHEETS IS REAL
In §III we assumed that the wavenumber k observed in buckling is purely real. Here, we demonstrate this is the case by supposing that Eq. (47) has complex roots. Since the tension τ 0 is real and Eq. (47) contains only even powers, there must be two complex conjugate pairs of solutions. We may write one pair as k ± = k r ± ik i , with k r ≥ 0, and hence the other pair will be −k ± . We extend sine and cosine to the complex plane using analytic continuation, and thereby extend the boundary conditions detailed in Appendix B to obtain the counterparts of Eqs. (52) and (53) as (k r + ik i ) tan(k r + ik i )a = (k r − ik i ) tan(k r − ik i )a (58) and (k r + ik i ) cot(k r + ik i )a = (k r − ik i ) cot(k r − ik i )a, where for both (58) and (59) the right-hand side is the complex conjugate of the left-hand side, so the imaginary part of either side must be zero. This condition takes the form f (k r a) = −g(k i a) and f (k r a) = g(k i a), for Eqs. (58) and (59), respectively, where f (x) ≡ x sin(x) cos(x) and g(x) ≡ x sinh(x) cosh(x) .
By plotting the functions f , g and −g one finds that their ranges do not overlap (though f and g have the same limit as x tends to zero) and hence there is no solution of Eq. (60). Therefore, k is real.

APPENDIX D: CROSSINGS OF SYMMETRIC AND ANTISYMMETRIC MODES
Here, we determine the values of a for which both Eqs. (52) and (53) are satisfied and thus when a symmetric mode and an antisymmetric mode have the same eigenvalue, which we refer to as a crossing.
We refer to the eigenvalues at these crossing points as type II, and they are given by The symmetric and antisymmetric modes come in pairs, which we label with the index j ∈ N ⋆ such that j = 1 corresponds to the pair with the smallest eigenvalues. We find that the crossings within the pair j occur for (m j , n j ) = (l − 1, l − 1 + j) and (m j ,ñ j ) = (l, l + j), l ∈ N ⋆ .
The index l labels the crossings within a given pair, such that l = 1 corresponds to the smallest size for which a crossing of either type occurs. Hence, we obtain two matrices of crossing points given by a (I) l,j = π 2 (2l − 1)(2l + 2j − 1) and a (II) l,j = π l(l + j) , l, j ∈ N ⋆ .
The corresponding eigenvalues are For very small values of a, within a pair the symmetric mode always has the smallest eigenvalue. As a increases, the first pair crosses at a = √ 3π/2 (type I), corresponding to τ 0 = 10/3. The next crossing (type II) occurs at a = √ 2π, corresponding to τ 0 = 5/2. Between those crossing points, the antisymmetric mode has the smallest eigenvalue. As a increases this pattern repeats infinitely many times (see figure 4).