The smectic order of wrinkles (2024)

Energy of wrinkled patterns

Before introducing the physical model, we must start with choosing a convenient coordinate system to represent the configuration of a thin wrinkled sheet. Such representation must be able to account for long-wavelength deformations to capture the overall geometry and, as we shall see, short wavelength wrinkles that take up the excess area. For this purpose we choose surface coordinates in which the configuration, assumed to be lying near the plane z=0, takes the form

where a and φ are scalar fields that in general depend on u and v. While this choice of Cartesian-like coordinates {u,v} may seem unintuitive, it can be rationalized by considering a uniform wrinkle pattern. Let the amplitude a(u,v) and wavenumber ∇φ(u,v) be constant, that is, a(u,v)=a and φ(u,v)=qu. Then the arclength of a curve along u is just

When the amplitude is small , ds/du is approximately constant while still capturing the overall excess length taken up by the wrinkle. With the representation in equation (1), slow variations in a and ∇φ do not alter the local wrinkle features. Calculations in terms of u and v are greatly simplified because the baseline, wrinkled metric does not have a short wavelength dependence on u, v. This is the major advantage of using these coordinates, and the reason for introducing the O(a²) correction to equation (1). In addition, it is clear from equation (1) that for small amplitudes the coordinates {u, v} are close to the Cartesian {x, y} coordinates of the projection of f onto the plane.

We will subsequently assume that the amplitude a(u,v) and phase φ(u,v) represent a slowly varying envelope of the wrinkled pattern, namely and . We further assume small amplitude wrinkles, namely the excess length to be . These assumptions directly relate to properties of the solution rather than to controlled parameters in the problem, but we will later show that they are consistent with results obtained from the model in a wide parameter regime and with solutions observed in experiments and simulations. In the following text we use the ∼ notation to indicate equality up to higher order terms in one or more of these small parameters.

With the above assumptions, the 2D metric and curvature tensors of configuration (1) take the form

with i, j∈{u,v}, and n the unit vector normal to the surface. Keeping our assumed scale separation in mind, we will be interested in looking at coarse-grained versions of different fields, namely ones which are averaged over oscillations at the wrinkle length scale |∇φ|⁻¹, and have only features which persist across at least several wrinkles. For that purpose we define a coarse-graining operator , where Ω is a region of radius ≳|∇φ|⁻¹ around each point and X can be any tensor field. Scale separation renders this definition unambiguous; the sheet’s reference metric (and of course the planar metric) can be approximated as Euclidean at the scale of Ω (since rad(Ω)<<D) and therefore integration, even of tensor fields, is well-defined. Since on the scale of Ω the envelope of a, ∇φ can be considered constant, while φ itself changes significantly, this coarse-graining operation will be manifested in the vanishing of all high-frequency modes ∝e^inφ(u,v), where n≥1, in X. For tensor fields this can be done element-wise in coordinates. Applying this to the metric (3) and curvature (4) tensors yields

Let us now assume a scenario in which a thin elastic non-Euclidean shell is put atop a bulk of liquid (Fig. 2). The system energy H takes contributions from both the elastic energy of the thin sheet—stretching and bending^17,18—and from the gravitational potential energy associated with displacing liquid elements (perfect wetting of the sheet is assumed). It can therefore be written in the form

where

For an isotropic material the stretching modulus is (with E Young’s modulus, t the thickness and ν Poisson’s ratio), bending modulus is , and substrate modulus is K=ρg (fluid density and gravitational acceleration, respectively). The parameters and are the shell’s reference metric and curvature tensors, respectively, and where and tr denotes contraction with the contravariant reference metric tensor . Note that the gravitational potential is not proportional to , and that it is not fully determined by the surface’s fundamental forms, because of the non-covariance of the liquid (in the 2D geometry). The difference between its proper area form dx dy and our coordinate area form du dv is high-order in our small parameters and will be neglected in the following. Substituting equation (5) into equations (7)–(9) yields

where (∇∇φ)_ij=∂_ijφ and (∇φ∇φ)_ij=∂_iφ∂_jφ (it is worth noting that ). Under our assumed scale separation, the Hamiltonian (6) (or any large-scale integral) is unchanged by coarse-graining its integrand, thus we can rewrite equations (6) and (10)–(12):

where

With our thinness and slowly varying envelope assumptions, the modulus Y that appears in the first term of equation (14) is much larger than the other moduli in equation (14). It follows that , hence that the value of the local excess length Δ=a²|∇φ|² is nearly constrained locally (see Supplementary Note 1 for a detailed derivation). Under this constraint, the bending and gravity terms in equation (14) are minimized at the well-studied finite wavenumber , around which we can expand these terms.

The elastic/smectic energy functional

We now rearrange the energy density (14) and write it as a sum of two new components:

The first component in equation (15), which we denote the elastic term, takes the form

where we find an effective reference metric, , that is slightly shifted from the actual reference metric of the shell towards a flat Euclidean metric according to the stretchability parameter η. The elastic energy component (16) penalizes deviations of the coarse-grained metric from this effective reference metric.

To illustrate the nature of this energy component, it is helpful to consider a wrinkled configuration of constant amplitude a. The coarse-grained metric (5) of such a configuration is identical to that of a surface whose topographic map is given by the wrinkled pattern, where wrinkles represent level-sets spaced in the z direction. Therefore, equation (16) penalizes deviations of that smooth surface’s geometry from the target geometry . However, this geometry is the coarse-grained version of many different wrinkle patterns, for example, taking a topographic map of the same surface from different projection angles or with a different choice of level-set spacing (wrinkle amplitude).

This degeneracy is broken by the second component in equation (15), the smectic term, which takes the form

The first term, smectic compression, penalizes deviations of wrinkle spacing from the preferred wavenumber k₀. Its modulus comes from the bending-gravity interplay, but is also proportional to the local excess length Δ which may therefore be thought of as a (perhaps spatially varying) local smectic order parameter. The second term, smectic curvature, penalizes wrinkle bending (Fig. 2d). Interestingly it has nothing to do with elastic bending, and actually comes only from residual tangent-to-sheet elastic stretching. A geometric intuition can be given to this term—unlike straight wrinkles, bent wrinkles have nontrivial Gaussian curvature which oscillates at half the wrinkle wavelength. If the target geometry is smooth at the wrinkle scale, that implies strain which results in elastic stretching, which is the source for the second term in (equation 17).

A very similar analysis can be performed when the substrate is soft elastic solid rather than liquid; the wrinkle wavenumber k₀ is then to be replaced with , where E_s is the substrate’s Young’s modulus and the substrate’s Poisson’s ratio is assumed 0.5 (ref. 8). The relevant dimensionless stretchability parameter η is , which does not depend on the thickness (as is expected from dimensional analysis). The following results hold in either physical scenario.

We next explore some of the outcomes of this analogy between wrinkles and smectics. to test our predictions, we employ two-dimensional (2D) finite-element simulations (Abaqus) of an elastic spherical shell with bulk gravitational force. The elastic and gravitational moduli, as well as the sheet geometry, can be tweaked to give us full control over the relevant dimensionless parameters (η, γ). Fluid surface tension exists in the experiment and is neglected in the simulation, however, in the relevant parameter regime its effect is negligible (Supplementary Fig. 2). Figure 2c,d demonstrates the good agreement in equilibrium configurations between the simulation and an experiment of a PDMS shell on water, performed with the same physical parameters. An extensive experimental and numerical research of such systems is in progress.

Domain walls

From the two smectic moduli in equation (17), one can extract a smectic penetration length

where λ₀ is the wrinkle wavelength. The penetration length λ_pen is a typical length over which smectic layers—or wrinkles—may bend smoothly. Hence, the wall between two smectic domains, with an angle difference between their wrinkle directions, will be of typical width 2λ_pen/tan ω (ref. 19). From equation (18) one learns that domain wall width is proportional to , hence proportional to the amplitude. In systems with a non-Euclidean underlying geometry, the amplitude may vary across the system (since ). Consequently domain walls may vary in width, possibly even along one continuous wall, as can be well seen in Fig. 3a. Remarkably, equation (18), based entirely on the theory of smectic liquid crystals, provides a precise prediction for the relation between the width of domain walls and local amplitude in wrinkled patterns. This prediction is in excellent agreement with measurements taken on our wrinkled elastic sheet simulations (Fig. 3c).

As Fig. 3 suggests, the above prediction appears to be valid only above certain values of the parameter η. Nonetheless, its breakdown is also explained by smectic theory; indeed, many systems that exhibit layered patterns, ranging from diblock copolymer melts^15,20,21 to Rayleigh–Bénard convection^22,23, show a transition from ‘Chevron’-shaped domain walls at small wall angles to Ω-shaped walls at large angles. This bifurcation is achieved via the formation of a pair of disclinations at each layer, to better maintain layer spacing in the system. A similar transition was observed, though to the best of our knowledge not explicitly discussed, in wrinkled herringbone patterns^24,25, seemingly via the formation of d-cones (ref. 26). The critical domain boundary angle for this transition is made wider when the thickness of the sheet (or equivalently η) increases, thus not allowing focusing of the curvature into small d-cones cores. Since the width of the Ω-shaped domain wall scales as the wrinkle wavelength λ₀, and not as the penetration length, one expects domain walls to abruptly change their scaling with decreasing η, as seen in Fig. 3. For wider domain wall angles, the Chevron-type walls and consequently the width ∝λ_pen relation are expected to persist onto narrower sheets.

Undulation instability

In smectic systems, certain boundary conditions may induce an undulation instability in the smectic structure^27,28. At its onset, the wavelength of this undulation is

where the d is a system scale at which boundary conditions are introduced. To test if similar behaviour exists in wrinkled systems, we designed a numerical experiment (Fig. 4) that follows the footsteps of ref. 27 in their smectic-A experiment. In our simulations, we hold a thin strip at its boundaries at a wavelength λ_B which is slightly mismatched with the preferred wrinkle spacing λ₀=0.95λ_B. One resolution for such mismatch may be an undulation, which allows layer spacing in the bulk to be smaller than on the boundary, while paying some cost for layer bending. The predicted optimum for this exchange is given by equation (19). We examine a sequence of systems with fixed geometry and wrinkle wavelength, however, with penetration/undulation lengths varying between simulations (achieved via changing both fluid density and sheet thickness so that λ₀ remains constant, however, η changes).

Indeed, an undulation in the wrinkling pattern appears (Fig. 4a). The observed undulation wavelength scales as the square root of the penetration length (Fig. 4b), as accurately predicted by equation (19). When using the ribbon width as a surrogate for the system scale d, we get a factor of ≈1.7 between the observed wavelength and the predicted one, asserting that ribbon width may not be the relevant system scale.

We further observe that wrinkle undulations spontaneously appear in certain regions of larger systems with free boundary conditions as a result of the global geometric incompatibility. We show examples for this observation in Fig. 4c and in Supplementary Fig. 1. We conjecture that this instability sets the scale for the herringbone patterns seen in biaxially compressed sheets on elastic substrates^25,29,30, as shown in Fig. 1d. The scaling implied by this conjecture agrees with the stability analysis scaling presented in ref. 31.

Domain size scaling

Last, our smectic analogy approach can further give a coarse prediction for the typical size and number of domains in systems that wrinkle as a result of geometric incompatibility. The coarse-grained metric g^(cg) of a domain-type patterns is approximately that of a piecewise-linear configuration that approximates the curved target metric . The strain that results from the difference between the two induces the elastic energy (16), therefore preferring smaller domains (that is, a better piecewise-linear approximation). However, the smectic energy term (17) is concentrated into domain walls, thus penalizing many-domain solutions. The interplay between these two contributions sets the typical domain size.

We consider a spherical shell with radius of curvature R. For domains of typical size L, the maximal (average) elastic strain scales as ∝ L²R⁻² (R⁻² is the Gaussian curvature mismatch). Integrated over the entire sheet of typical size D, this will result in total elastic energy . The energy per-unit-length of a domain wall is proportional to the smectic curvature modulus divided by the smectic penetration length¹⁹, that is, scales as ∝Δ^3/2η^1/2λ₀∝Δ^3/2t. Upon introducing the typical excess length implied by the global geometric incompatibility, , and the total length of the walls ∝D²/L, we get the total smectic contribution . Minimization of the total energy (15) suggests the scaling

where γ is the Föppl–von Kármán number. This general scaling relation should be taken with caution, since the number of domains is discrete and is very much affected by boundary conditions. A quantitative experimental or numerical verification of equation (20) must therefore involve very large systems (so that boundary effects are irrelevant), with typical domain size extracted via statistical analyses of the pattern. Another experimental complication involves the very small exponent by which the domain size in equation (20) depends on physical parameters. Preliminary numerical observations may support relation (20) (Fig. 5), however, a robust numerical verification using a finite-element thin-sheet simulation is computationally inconceivable as a result of the extreme multi-scale nature of the problem. Nevertheless, we hope that a more efficient simulation based on our phase field φ approach may in the future allow verification of this scaling.