Surface sulci in squeezed soft solids

The squeezing of soft solids, the constrained growth of biological tissues, and the swelling of soft elastic solids such as gels can generate large compressive stresses at their surfaces. This causes the otherwise smooth surface of such a solid to becomes unstable when its stress exceeds a critical value. Previous analyses of the surface instability have assumed two-dimensional plane-strain conditions, but in experiments isotropic stresses often lead to complex three-dimensional sulcification patterns. Here we show how such diverse morphologies arise by numerically modeling the lateral compression of a rigidly clamped elastic layer. For incompressible solids, close to the instability threshold, sulci appear as I-shaped lines aligned orthogonally with their neighbors; at higher compressions they are Y-shaped and prefer a hexagonal arrangement. In contrast, highly compressible solids when squeezed show only one sulcified phase characterized by a hexagonal sulcus network.

The squeezing of soft solids, the constrained growth of biological tissues, and the swelling of soft elastic solids such as gels can generate large compressive stresses at their surfaces. This causes the otherwise smooth surface of such a solid to becomes unstable when its stress exceeds a critical value. Previous analyses of the surface instability have assumed two-dimensional plane-strain conditions, but in experiments isotropic stresses often lead to complex three-dimensional sulcification patterns. Here we show how such diverse morphologies arise by numerically modeling the lateral compression of a rigidly clamped elastic layer. For incompressible solids, close to the instability threshold, sulci appear as I-shaped lines aligned orthogonally with their neighbors; at higher compressions they are Y-shaped and prefer a hexagonal arrangement. In contrast, highly compressible solids when squeezed show only one sulcified phase characterized by a hexagonal sulcus network. Complex patterns often arise from simple causes, in such instances as the fractal structures in physical aggregation phenomena [1] or the labyrinthine structures, spots and stripes in chemical systems [2]. In purely mechanical systems, there are two basic instabilities associated with the buckling of a slender filament or sheet [3,4] and the cracking of a bulk solid [5,6]. The first instability arises because of a competition between compression and bending and is embodied in the ratio of two length scales, while the second arises because of a competition between bulk and surface effects and embodied in the ratio of two energy scales. Here we show that an even simpler system -a thick isotropic, homogeneous elastic layer which is subject to planar compression and whose top surface is free -is susceptible to the formation of various sulcal patterns (cusped folds) as a function of the applied compressive strain. This process is based on a surface instability with no length scale, in sharp contrast to the thin film based pattern formation [7] where the film thickness is the intrinsic length scale.
The elastic instability of a compressed surface was first studied by Biot [8] who showed that a half-space of incompressible neo-Hookean material becomes unstable to any smooth perturbation when the surface stretch ratios λ x and λ y (defined so that an uncompressed surface has λ x = λ y = 1) reach a critical value λ x λ y ≈ 0.544. More recent studies have found that a subcritical instability in the form of a sulcus is energetically favorable when λ x λ y < ∼ 0.647 [9], i.e. a nonlinearly unstable sulcified state with no nucleation threshold can exist at lower compression than predicted by Biot's linear analysis. Apart from some early and preliminary studies in gels [10,11], theoretical and numerical studies of sulcification [9,[12][13][14][15][16] assume two-dimensional plane-strain conditions, and hence exclude all realistic sulcal morphologies that involve three-dimensional deformations seen in most phys-ical experiments [17,18]. These three-dimensional morphologies are also seen in biological tissues, including particularly prominent sulci in the primate brain [19], but also in tumors [12] and other organs, where there is evidence of mechanical forces driving folding [20]. In addition to being a new paradigm for mechanical pattern formation, this system also serves as a model for unusual thermodynamical phase transitions without a barrier [9], classical nucleation [21], interface instabilities etc.
Experimental observations of sulcal patterns can be realized by isotropic compression induced, for example, by swelling gels [17,18,22] and typically lead to threedimensional patterns that are hard to control. Numerical simulations provide means for controlled studies, but simulating sulcification in three dimensions is challenging as it involves finite deformations, tracking the free boundaries of the unknown sulci, and accounting for selfcontact of the free surface. To overcome these difficulties, we use a finite element method to approximate the solid with a dense rectangular mesh of tetrahedrons [23], based on a discretized finite strain elasticity theory with a neo-Hookean energy density where F is the deformation gradient, J = det(F), and µ and K are shear and bulk modulus, respectively. To relax the integrated energy density (1) towards minima, we use damped Newtonian dynamics of the nodal degrees of freedom. Self-avoidance is implemented by short range repulsive contacts between the edges and faces that make up the surface. The numerical model is well suited for sulcification studies as Newtonian dynamics naturally allows the system to change rapidly when there are large configuration changes while the quadratic elastic potentials allow for stable integration of the equations of motion with a fixed time step. To confirm that our results are independent of the mesh geometry, we also used a hexagonal prism mesh and checked that our results are robust [23].
To study sulcification under isotropic compression, we simulate a layer with thickness h in the stress-free state, and lateral dimensions corresponding to a square of side L = 10h; the lattice has ∼ 320 × 320 × 40 nodes with ≈ 2 × 10 7 elements. The base of the layer is clamped and periodic boundary conditions are applied along the lateral edges to prevent edge effects from constraining sulcal morpologies. Since many soft materials are approximately incompressible, we assume that the bulk modulus K = 30µ. Our simulations start from an isotropic compressed reference state with stretch λ = λ x = λ y = 0.54 in both planar directions, corresponding to strain = λ − 1 = −0.46 and λ x λ y = 0.40 well beyond the Biotpoint (λ x λ y ≈ 0.544 [8]) and T-point (λ x λ y ≈ 0.647 [9]). To trigger sulcification the featureless flat surface is perturbed by small random vertical displacements of the nodes (maximum amplitude 10 −3 × lattice constant), after which the system is allowed to relax to a sulcified state while keeping the strain constant. In Fig. 1a we see the appearance of a densely sulcified state characterized by isolated Y-shaped triple-junctions of sulci lying on an approximately triangular lattice. The three-fold symmetry of the junctions is consistent with their angle being ≈ 120 • ; occasionally some triple-junctions share an arm with their neighbor, although usually they are isolated. Once the layer is fully relaxed in its equi-biaxially compressed state, we quasistatically decompress it while preserving the lateral stress isotropy (see movie [23]). Decompression of the layer transforms the Y-shaped sulci, one by one, to I-shaped sulci. In Fig. 1b we show the layer at λ = 0.61 where Y's and I's coexist. At λ = 0.67 (Fig. 1c), at the Biot point, all Y's have transformed to I's, but their number remains approximately constant. As the layer is allowed to relax still more, the I's shorten and eventually unfold at λ ≈ 0.73 as the T-point is approached (Fig. 1d). We use unloading rather than loading to probe the patterns to circumvent that the absence of non-smooth perturbations prevents sulcus formation during loading until the Biot-point [8] is reached, in contrast with any physical experiment where the T-point determines sulcus formation [9].
To understand the patterns qualitatively, we note that the formation of a sulcus relaxes stress primarily in the direction perpendicular to it (similar to a crack), so that it is unfavorable for adjacent sulci to be parallel. The arrangement of the I-shaped sulci in a square lattice with alternating orientations is a natural solution. Indeed, in Fig. 1 we see such a pattern although the system may get trapped in metastable states that break this order at times. For Y-shaped sulci, on the other hand, the simplest plane-filling symmetric pattern is based on a hexagonal lattice, although again we see imperfections at high compression in our unfolding simulations; we note that these symmetries are not those of the underlying lattice [23].
Perfectly ordered sulcus patterns can be constructed numerically by using a "mold" to imprint the desired pattern onto a precompressed layer. Given the hysteretic nature of the sulficification transition [9], when the mold is removed, the pattern persists. In Fig. 2a  conditions, simulating a unit cell with two horizontal and two vertical I's on a lattice of dimensions ∼ 125×125×40 with a simulation domain of unknown side length a s . For each simulation we fix λ and vary a s to find the energy minimum. In Fig. 2b, we show the hexagonal Y-pattern and again find the optimal hexagonal spacing a h for a given fixed compression. The surface profile of an I-sulcus shown in Fig. 2a reveals that it is symmetric while the surface profile of a Y-sulcus (Fig. 2b) reveals that it has a sharp asymmetric center. A comparison of the energies of these patterns with those of freely relaxed layers in Fig. 2c shows that the perfectly ordered patterns have lower energy. Furthermore, we find that I-patterns are stable for λ > ∼ 0.61 and Y-patterns are stable for λ < ∼ 0.63, while in the narrow regime 0.61 ≤ λ ≤ 0.63 their energies are nearly equal, and explains the coexistence of Y's and I's seen in Fig. 1b. The optimal spacing of these structures is given by the relations a s ≈ 3h and a h ≈ 2h but it increases weakly with compression in both cases. The square and hexagonal symmetries of the sulcus patterns are similar to those seen in other pattern forming systems such as fluid convection [24] and elastic fracture [25]; however, there are fundamental differences from a mechanistic perspective, and we will not pursue the mathematical analogies further here.
Sulcification is an energetic consequence of the exchange of stability between a uniformly deformed state and a set of localized states. To quantify this, in Figs. 3a,b we show the depth-averaged energy distributions of the I-and Y-patterns, obtained by integrating energy density over the thickness of the layer in material coordinates. They reveal that Y's span triangular areas on the surface, whereas I's span elliptical areas. Although sulcified states are favored in terms of total energy, between sulci the energy density increases with respect to the reference state. Similarly, energy distribution in the thickness direction (Fig. 3c) shows that, while relaxing compression relieves energy near thesurface, the material at some depth gains energy as it conforms to the buckling surface.
Having considered isotropic compression, we now look at anisotropic compression, characterized by the strain λ x λ y ≈ 0.65 λ x λ y = 0.54 ratio y / x < 1 ( x = λ x − 1, y = λ y − 1). We perform these simulations by starting with a compressed reference state and then unloading the layer quasistatically keeping y / x constant. The lattice size and simulation domain are as in the simulation of Fig. 1. In the perfectly anisotropic case ( y / x = 0, see Fig. 4) as expected, we see stripes of sulci in the direction perpendicular to the direction of compression. By including compression in the y-direction we find that when y / x ≈ 0.5 stripes begin to break up. When the strains are set almost equal ( y / x = 0.85 in Fig. 4), Y-shaped sulci appear, but anisotropy is still apparent from the pattern. As in the isotropic case ( y / x = 1), Y's transform to I's with decompression, but now all the I's are oriented perpendicular to the direction of highest compression. We observe an almost identical unfolding threshold λ x λ y ≈ 0.62 for all y / x ; this deviates from previous numerical results for the plane strain case [9] where sulci unfold at λ x λ y ≈ 0.647 because of the weak dependence of the critical strain on the finite mesh size; here sulci vanish when their size become comparable to the mesh spacing. The near threshold behavior of these sulci calls for a more careful analysis, but the transversely isotropic and planestrain cases have similar hysteresis effects associated with the presence of two critical points [9].
Having characterized the patterns on the surface of in-  compressible hyperelastic solids, we turn briefly to consider the other limit of soft surfaces of highly compressible materials, such as solid foams. Here, the bulk modulus K ∼ µ enters as a relevant parameter, and the Poisson effect that couples the transverse directions is weaker than in the nearly incompressible materials considered above. When such a solid is isotropically compressed, we find that sulci form a connected hexagonal network, typically with some imperfections, as shown in Fig. 5a. For K < ∼ 2µ the hexagons persist all the way to the Tpoint upon unloading, which itself shifts to higher strain with decreasing K. These findings can be summarized in a simple phase diagram (Fig. 5b) of sulcus morphologies as a function of bulk modulus and compression. We note that the phase boundaries in the diagram are only qualitative guides since the simulations indicate regions of coexistence of the different morphological states.
Our simulated patterns are able to capture the range of experimental observations of sulcification in swollen gel layers [17] shown in Figs. 1e and 1f. Indeed, prior observations show the domains of ordered I patterns and YI-mixtures, reminiscent of our freely relaxed layers, although the transition from I's to Y's has not been previously attributed directly to increasing compression. Furthermore, the spacing between sulci in our simulations agrees well with the experimentally observed spacing [17] with a similar weak strain dependence as observed for uniaxially compressed sulci [26]. Sulcus patterns in compressible hydrogels (which are poroelastic and thus compressible over long time scales) have been observed to relax with time into honeycomb structures [10,11] that are similar to our compressible hexagon patterns. In a biological setting, several organs, including the cerebral cortex and cerebellum in the brain [19], have sulcified surfaces, with the cerebellum showing striped patterns while the cortex showing triple-junctions that are similar to those seen in our simulations. Recent experiments [20] show the presence of residual strains in these tissues, consistent with the hypothesis that sulicfication might be a simple consequence of relative growth. From a technological perspective, since sulci form so easily on the surface of soft solids, they should be easy to manipulate as well. Efforts to design and control smart surfaces using temperature-responsive gels [27], voltage-responsive elastomers [28] and mechanical strain are just beginning, and point the way to functional patterning via sulcification.
We thank E. Hohlfeld and R. C. Hayward for discussions. TT acknowledges the Academy of Finland for funding. The computational resources were provided by the Finnish IT Center for Science (CSC).
Supplementary material for 'Surface sulci in squeezed soft solids'

MESH GEOMETRY
A layer with thickness h in the stress-free state is confined to a square domain of width W . The surface of the layer is assumed to be normal to the y-direction and its base is clamped. We assume periodic boundary conditions along the the edges in the x and z-directions. The layer is discretized into a rectangular mesh, and each rectangle is divided into five tetrahedrons as indicated in Fig. 6. The arrangement of the tetrahedrons in any two neighboring rectangles is reflected with respect to the face they share; this imposes mesh symmetry with respect to reflections in the x, y and z-directions. The number of nodes in the planar and vertical directions is adjusted so that each edge of the rectangle has length a.

TETRAHEDRON ELEMENT
In its stress-free configuration, each tetrahedron is defined by its natural state in terms of the matrix wherex 1 ,x 2 andx 3 are vectors describing the tetrahedron, see Fig. 7a. The growth, or swelling of the tetrahedral elements is characterized by the tensor where g x , for example, indicates the stress-free growth ratio in the initial x-direction. Then the deformed configuration (Fig. 7b) of the tetrahedron, including growth, is characterized by where x 1 , x 2 and x 3 are the deformed basis vectors. Here F is the elastic deformation gradient, and we have assumed a multiplicative decomposition of the total deformation gradient A into a form similar to that used in finite strain plasticity or growth processes. This allows us to obtain F from eq. (4) by using the relation To connect the kinematic relations linking deformation and stress, we assume that our elastomer may be modeled as a compressible neo-Hookean material with a strain energy density (in its undeformed configuration) where µ and K are the shear and bulk modulus, respectively, and J = det(F). The corresponding Cauchy stress, i.e., the force per unit area in the deformed configuration, is given by We note that this allows us to characterize each element in terms of its local constitutive behavior, in contrast with other structural models that use beams, plates or shells as the elementary units constituting the solid. Furthermore, for the purposes of computing the equilibrium configurations, we also assume that each element also has internal damping. Since the relative velocities of the nodes are defined by the rate of deformation is given by and the viscous stress by where the elemental viscosity is η Then, the traction on each deformed face (i = 1, 2, 3, 4) of the tetrahedron is given by where is the total stress and n i are normals with lengths proportional to the deformed areas of the faces, see Fig. 7b. Nodal forces are obtained by distributing the traction of each face equally on to its three vertices.

SELF-AVOIDANCE OF THE SURFACE
Since our computations involve self-contact at the sulci, we must ensure that this is taken care of correctly. As the element faces at the surface of the layer form a lattice of triangles, self-avoidance is accommodated by processing 1) vertex-triangle contacts and 2) edge-edge contacts. If a separation d between a vertex and triangle is less than the contact range h = a/3 it is considered a contact. Contacts are penalized by an energy Ka 2 h−d

DAMPED DYNAMICS FOR ENERGY MINIMIZATION
We use damped second order dynamics for energy minimization. After the nodal forces are determined, Newton's equations of motion are solved for the nodes by an explicit scheme, x(t + ∆t) =x(t) + v(t + ∆t)∆t.
Here ∆t = 0.1a/ √ K is the time step, m = a 3 mass of a node, and γ = m viscous damping. Vectors f , v and x are force, velocity and position of a node, respectively.

ALTERNATIVE DISCRETIZATION: HEXAGONAL MESH
To explore possible discretization artifacts on the sulcal patterns, we also implemented our simulations on a hexagonal prism mesh as shown in Fig. 8. Each prism is divided to three tetrahedrons whose arrengement in any two neighboring prisms is reflected. The reason for this is that division of a prism into tetrahedrons breaks its symmetry, but by the altering arrengement the symmetry is recovered at the level of the whole mesh. We simulate spontaneous sulcification from a compressed reference state followed by quasistatic decompression. We obtain similar patterns , including the transition between Y-and I-phases, as with the rectangular mesh used in the simulations of the main text. The lattice only affects the orientation of the patterns, which simply reflects the fact that in a true continuum system the orientation is arbitrary, while the lattice effects can break this degeneracy to determine the global orientation of the pattern but not its symmetry or morphology. We also simulated full loading/unloading cycles with both meshes to confirm that their energies and thresholds are similar. Figure 9 indicates that under loading the surfaces remain flat, with energy equal to the reference state, until Biot-threshold is reached. Under unloading sulci persisted until λ ≈ 0.729 on a rectangular lattice and λ ≈ 0.722 on a hexagonal lattice. The slightly higher unfolding strain of the hexagonal mesh is most likely due to the mismatch between the mesh geometry and preferred orthogonal arrangement of sulci near the threshold.