# Unfolding the Sulcus

Sulci are usually seen in combination with complementary protrusions known as gyri on the surface of the primate brain, but are also seen on the palm of our hand, in our elbows and knees, in swollen cellular foams (such as bread) and gels, and in geological strata; a few representative examples are shown in Fig. 1. While observations of sulci are ancient, their systematic study is fairly recent; an early reference is the reticulation patterns in photographic gelatin [1], and there has been a small interest in these objects both experimentally [2–6] and theoretically [3,7–9], starting with the pioneering work of Biot [10] over the past 50 years. Despite this, there is no careful analysis of the fundamental instability and bifurcation that leads to sulci. Here we study the formation of a sulcus in a bent slab of soft elastomer, e.g., PDMS: as the slab is bent strongly, it pops while forming a sulcus that is visible in the lower right panel of Fig. 1; releasing the bend causes the sulcus to vanish continuously, in sharp contrast with familiar hysteretic instabilities that pop in both directions. We find that sulcification is a fundamentally new kind of nonlinear subcritical surface instability with no scale and a strongly topological character, yet has no energetic barrier relative to an entire manifold of linearly stable solutions. We also argue that sulcification instabilities are relevant to the stability of soft interfaces generally, and provide one of the first physical examples of the consequences of violating the complementing condition [11] (during loading) and quasiconvexity at the boundary [12] (during unloading), keystones in the mathematical theory of elliptic partial differential equations and the calculus of variations.

To understand the unusual nature of sulcification, we first recall Biot’s calculation for the linear instability of the surface of an infinite half-space of an incompressible elastomer that is uniformly stressed laterally. Because the free surface is the softest part of the system, and there is no characteristic length scale in the equations of elasticity or in the boundary conditions, instability first arises when the Rayleigh surface wave speed vanishes. For incompressible rubber, Biot [10] showed that this threshold is reached when the compressive strain exceeds 45.3% [10], at which value, all surface modes are unstable while the fastest growing one has an infinite wave number. Since every free surface looks like a half-space locally, Biot’s instability lurks at every free boundary. Finite geometries typically break this infinite degeneracy and lead to a hierarchy of ordinary buckling instabilities that preempt the surface instability [13]. However, since Biot’s calculation was limited to a linear analysis of the problem, it could not address the question of whether the instability was supercritical or subcritical or its ultimate nonlinear saturation. Since the basic problem is scale-free and translation invariant (the sulcification instability can arise anywhere along the surface), the nonlinear problem is numerically intractable without explicit regularization and a careful limiting process requiring that we unfold the sulcus literally and figuratively.

Therefore we consider the bent strip geometry shown in Fig. 1(d), and break scale invariance by assuming that a thin skin of a stiff material is attached to the surface of the bent slab. Furthermore, the curvature maximum at the bottom of the horseshoe where the highest strains are achieved naturally breaks translation invariance. For planar deformations, the total energy of the system is given by

where $\mathbf{x}(\mathbf{X})$ is the deformation of a strip occupying a rectangular material volume $\Omega \subset {\mathbb{R}}^{2}$ and subject to the incompressibility constraint $\mathrm{det}(\nabla \mathbf{x})=1$, ${\mu}_{0}$ is the shear modulus of the incompressible elastomer, ${B}_{0}$ is the stiffness of a semiflexible skin and $s$ is the arc length parameter of the upper surface $\Gamma \subset \partial \Omega $. Our model corresponds to having a simple neo-Hookean elastomer free energy for the bulk and a Bernoulli-Euler curvature energy for the skin. To understand the onset of the sulcification instability, we extremized the energy above using a custom-built finite element method with continuous strains and a hierarchical mesh (see supplementary information (SI), Ref. [14]). We enforced incompressibility and self-contact using pressure fields and by assuming left-right symmetry about the sulcus. This model has three relevant length scales: a regularization length ${l}_{r}=\sqrt[3]{\frac{{B}_{0}}{{\mu}_{0}}}$, the length of the strip ${L}_{s}$, and its thickness ${W}_{s}$. We use ${\mu}_{0}$ and ${L}_{s}$ to scale all quantities so that $B=\frac{{B}_{0}}{{\mu}_{0}{L}_{s}^{3}}$ and the aspect ratio of the strip ${L}_{s}/{W}_{s}$ are the only dimensionless parameters.

Our simulations start with an initially flat rubber strip that is bent and quasistatically compressed between parallel, rigid plates separated by a distance $\Delta $. As the control parameter $\Delta $ is decreased, compressive strain on the inner surface of the strip increases and ultimately drives sulcification. To ensure that the scale of the furrow is not numerically under-resolved, we use a recursively refined finite element mesh near the incipient sulcus to keep the mesh scale roughly an order of magnitude smaller than ${l}_{r}$ (see SI, Sec. B). Using a novel continuation method for variational inequalities (see SI, Sec. C), we computed both stable and unstable extrema of $E(\mathbf{x})$, and explored the limit $B\to 0$. This yields the central result of our study, the family of bifurcation diagrams shown in Fig. 2(a), where we plot the minimum height of the slab $h$ as a function of $\Delta $.

Each $h-\Delta $ curve is a bifurcation diagram for a different value of $B$: solid lines represent linearly stable solutions, while the dotted lines represent linearly unstable solutions. Each curve follows the characteristic $S$ shape of a hysteretic transition, associated with a sudden change in $h$ and formation or relaxation of a finite size sulcus when a critical value of $\Delta $ is passed in loading or unloading. For every value of $B$, extrema on the top branch have smooth surfaces, those on the middle branch have a pendant of size ${l}_{r}$ [see inset a of Fig. 2(a)], and those on the bottom branch have a self-contacting sulcus (insets b and c). As $B$ is decreased, ( $B\in [{10}^{-7}-{10}^{-12}]$, ${l}_{r}/{L}_{s}\in [4.6\times {10}^{-3}-{10}^{-4}]$), the hysteresis in a typical loading cycle becomes atypically one-sided as the branch of unstable extrema (dotted lines) swings up toward the top stable branch and the $S$-shaped bifurcation diagrams converge to a master $T$-shaped diagram traced by the thick gray line with two critical points.

For fixed $\Delta $ and $B$ an unstable solution represents a saddle point in the energy landscape; its energy relative to the configuration with a flat surface is $\Delta E={E}_{u}(\Delta )-{E}_{s}(\Delta )$ where ${E}_{u}(\Delta )$ and ${E}_{s}(\Delta )$ are the energies of the unstable and top, stable branches at $\Delta $ respectively, and is an upper bound on the height of the barrier to nucleating a sulcus. Figure 2(c) shows $\Delta E$ as a function of $B$ and confirms the convergence of the family of bifurcation diagrams toward the limiting $T$-shaped bifurcation diagram, as well as the existence of a nonlinear surface instability with no energetic barrier over an extended range of strains. We see that the instability thus differs significantly from traditional first order phase transitions in that the deformation is continuous, occurs in simple elastic continua, and is well defined in the limit of vanishing surface energy. The presence of a metastable region bracketed by a pair of critical strains along which the stable and unstable solutions coincide as the skin becomes vanishingly thin naturally explains the discrepancy between Biot’s prediction and a large number of experiments on creasing and sulcification [4] (and references therein) over the past half century. Unfolding the instability without breaking translation symmetry at the surface, e.g., in a swollen, adhered layer of gel, then naturally leads to extreme sensitivity to imperfections, and a hierarchy of complex subcritical instabilities connecting Biot’s instability and buckling (see SI, sec. E), and the ability to control sulcification [6].

As $B\to 0$, the sequence of saddle-node fold points encountered during loading converges to a limiting, infinitely sharp fold point when the surface strain at the lowest point on the inner surface of the horseshoe, $\mathbf{X}=\mathbf{0}$, reaches a critical value of 45.6% consistent with Biot’s classical result; in Fig. 2(b) we see the convergence of the critical compressive stain, $1-{\lambda}_{x}$ where ${\lambda}_{x}$ is the principal stretch of the deformation gradient $\nabla \mathbf{x}$ along the free surface, for finite $B$ to Biot’s predicted value at $B=0$. A numerical linearized spectral analysis of the loaded slab also confirmed Biot’s prediction that the Rayleigh surface wave speed vanishes at $\mathbf{X}=\mathbf{0}$ just as the critical strain is achieved, and corresponds to the failure of the complementing condition [11,15], wherein infinitesimal periodic solutions at the boundary grow at a rate that diverges as the inverse of the wavelength. Over-damped dynamical simulations—which trace steepest descent contours of the energy landscape—reveal that nonlinear effects reorganize these surface waves into a self-similar furrowing process; after a short transient, depending on $B$, the growth of a sulcus, which occurs via rolling, not snapping, is described by the self-similar form ${\mathbf{x}}_{s}(\mathbf{X},{\Delta}^{\star})+\sqrt{\lambda t}{\mathbf{v}}^{\star}(\mathbf{X}/\sqrt{\lambda t})$ where $\lambda $ is a dimensional constant and ${\mathbf{v}}^{\star}$ is the numerically computed scale-invariant form of the sulcus [15], and ${\mathbf{x}}_{s}(\mathbf{X},\Delta )$ is the branch of smooth solutions and ${\Delta}^{\star}$ is the value of $\Delta $ at Biot’s limiting fold point.

The complementary sequence of saddle-node fold points for decreasing $B$ encountered during unloading are actually “corners” associated with the loss or gain of a self-contacting sulcus as the unstable surface pendant just closes to form a cavity of fixed size ${l}_{r}$. As these corners converge to the limiting “ $T$ point,” the maximal surface strain outside the self-contacting region approaches a limiting value of 35.4% that is attained at a sequence of points converging to $\mathbf{X}=\mathbf{0}$. The convergence to this strain is traced by the lowest curve in Fig. 2(b) with the asymptote marked by the dashed line. The middle solid curve of Fig. 2(b) is another estimate of the critical strain computed by measuring the strain at $\mathbf{X}=\mathbf{0}$ for a sequence of extrema for corresponding values of $\Delta $ on the top branch. The $T$-point critical strain (like the Biot critical strain) is universal for free surfaces of incompressible materials, consistent with recent experimental observations [4,6]; however they both change with applied normal stress (i.e., indentation), for compressible materials [15] etc.

To understand why the $T$-point bifurcation and the entire unstable manifold are not captured by linearized analysis, we note that before the sulcus reaches the regularization scale ${l}_{r}$, it shrinks according to the form ${\mathbf{x}}_{s}(\mathbf{X},\Delta )+l(\Delta ){\mathbf{v}}_{T}\mathbf{(}\mathbf{X}/l(\Delta )\mathbf{)}$ where $l(\Delta )\ge 0$ vanishes at the $T$ point and ${\mathbf{v}}_{T}\approx {\mathbf{v}}^{\star}$. (See insets b and c in Fig. 2(a), and the relative scale factor of 6.5.) Since the elastic stress is determined by $\nabla \mathbf{x}$, this transformation shrinks the size of the sulcus without altering the local stress balance; therefore all the material and contact nonlinearities remain relevant even for vanishingly small sulci.

We tested our theory with experiments using a commercial Sylgard 180 Elastomer to form $36\times 26\times 4\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ slabs that were placed between parallel rigid plates attached to linear motors and compressed in small increments of $200\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ in a second, separated by 50 s to allow for the equilibration of the slab. We tracked sulcification optically by imaging the refracted image of a laser sheet that passed through the slab along its bending axis. When the sulcus formed it sharply refocused the laser sheet into an almond-shaped caustic pattern surrounding a dark shadow (SI, sec. D). Figure 3 (left) shows the evolution of a central raster scan of the caustic pattern during a loading cycle (vertical axis) (see SI, sec. D). Analogous ray traced light distributions for the simulation, using the measured laser profile and assuming left-right symmetry, the measured system geometry, and $B={10}^{-11}$ (physically ${l}_{r}\approx 18\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$), are also shown in Fig. 3 (right) for comparison; the numbered red dots correspond to the numbered red dots in Fig. 2(a). We see that with no free parameters we can capture the one-way hysteretic transition associated with sulcification.

The emergence of the $T$ and Biot points, and intervening metastable region in the $B\to 0$, can be understood in terms of a *nonlinear* generalization of Biot’s half-space problem. All our simulations show that when a half-space of incompressible elastomer is compressed by 34.5% it has an infinite degeneracy of energy minimizers: the trivial flat configuration, and a continuous family of isolated sulci which are stable up to translation and rescaling (i.e., $\mathbf{v}(\mathbf{X})\to l\mathbf{v}(\mathbf{X}/l)$ for any number $l>0$), i.e., these symmetries are spontaneously broken. Sulcification exchanges compressive strain for rotation and shear which are ultimately limited by self-contact. Beyond the lower-critical strain, forming a sulcus of size $O(l)\ll 1$ releases energy over a region of size $O({l}^{2})$, equivalent to the failure of quasiconvexity at the boundary. The spatial variations of $\nabla \mathbf{x}$ near the compressive strain maximum at $\mathbf{X}=\mathbf{0}$ act as symmetry-breaking perturbations and determine the ultimate scale $l$ of the surface fold. When the compressive strain is not localized to a point, the size of the sulcus is not set by the local geometry of $\nabla \mathbf{x}$ and the domain, and the $T$ bifurcation is sensitive to details and potential interactions between multiple sulci resulting in reticulation [1,3], or in a combination of buckling and sulcification [16,17] (see SI, sec. E).

More generally, $T$ bifurcations might arise in elastic systems with internal interfaces and nucleationlike processes in elliptic systems where nonlinearities enter in a scale-free way, e.g., the formation of cavities, bubbles and cracks [18,19]. These processes are notoriously difficult to control, displaying extreme sensitivity to imperfections, and are associated with a discontinuous transition in the microscopic state characterized by a critical size nucleus; e.g, a bubble or crack will grow only once it has reached a threshold size. $T$ points should exist in these systems in the limit when the surface energy vanishes and the size of the “defect” also vanishes, but the ratio of the two which corresponds to a critical pressure or stress remains finite.

## References

- S. E. Sheppard and F. A. Elliot, Ind. Eng. Chem.
**10**, 727 (1918). - E. Southern and A. G. Thomas, J. Polym. Sci., Part A: Gen. Pap.
**3**, 641 (1965). - T. Tanaka
*et al.*, in*Molecular Conformation and Dynamics of Macromolecules in Condensed Systems*, edited by M. Nagasawa (Elsevier, New York, 1988). - A. N. Gent and I. S. Cho, Rubber Chem. Technol.
**72**, 253 (1999). - A. Ghatak and A. L. Das, Phys. Rev. Lett.
**99**, 076101 (2007). - V. Trujillo, J. Kim, and R. C. Hayward, Soft Matter
**4**, 564 (2008). - T. Hwa and M. Kardar, Phys. Rev. Lett.
**61**, 106 (1988). - A. Onuki, Phys. Rev. A
**39**, 5932 (1989). - S. A. Silling, J. Appl. Mech.
**58**, 70 (1991). - M. A. Biot,
*Mechanics of Incremental Deformations*(John Wiley and Sons, New York, 1965). - S. Agmon, A. Douglis, and L. Nirenberg, Commun. Pure Appl. Math.
**17**, 35 (1964). - J. M. Ball and J. E. Marsden, Arch. Ration. Mech. Anal.
**86**, 251 (1984). -
Thus the Biot instability is actually the progenitor of all structural or buckling instabilities.

- See supplemental material at http://link.aps.org/supplemental/10.1103/PhysRevLett.106.105702.
- Evan Hohlfeld, Ph.D. thesis, Harvard University, 2008.
- L. Pocivavsek, R. Dellsy, and A. Kern
*et al.*, Science**320**, 912 (2008). - E. Sultan and A. Boudaoud, J. Appl. Mech.
**75**, 051002 (2008). - A. N. Gent and P. B. Lindley, Proc. R. Soc. A
**249**, 195 (1959). - J. M. Ball, Phil. Trans. R. Soc. A
**306**, 557 (1982). - J. Bradbury, PLoS Biol.
**3**, e50 (2005).