Surface creasing of soft elastic continua as a Kosterlitz-Thouless transition

TL;DR

This paper models surface creasing in soft elastic materials as a Kosterlitz-Thouless transition, identifying ghost fibers as quasi-particles that emerge at a critical strain, aligning with experimental observations.

Contribution

It introduces a novel theoretical framework linking surface creasing to topological phase transitions, providing a unified explanation for crease formation and profiles.

Findings

Creases arise from unbinding of ghost fiber dipoles at a critical strain.

Predicted surface profiles match experimental data.

Proposes further experimental validation of the model.

Abstract

Harnessing a model from composite materials science, we show how point-like cusped surface features arise as quasi-particle excitations, termed "ghost fibers", on the surface of a homogeneous soft elastic material. These deformations appear above a critical compressive strain at which ghost fiber dipoles unbind, analogous to vortices in the Kosterlitz-Thouless transition. Finite-length creases can be described in the same framework. Our predictions for crease surface profiles and onset strain agree with previous experiments and simulations, and further experimental tests are proposed.

Figures

c in Trujillo *et al.* ([2008](#bib.bib5)) exhibits both motifs.) Within our framework, such patterns again have a natural interpretation as charge-crystallized shear lag quasi-particles. In the limit where the spatial extent of quasi-particles is very small compared with their lattice spacing (i.e., they are ghost fiber-like), the creasing pattern is predicted to be a 2d Coulomb crystal with hexagonal symmetry (e.g. Bonitz *et al.* ([2008](#bib.bib42))). In fact the hexagonal dimple crystal observed in Brojan *et al.* ([2015](#bib.bib25)); Jiménez *et al.* ([2016](#bib.bib26)) is suggestive o

Equations

$$\oint \vec{\gamma }\cdot d\mathbf{A}=\frac{f_0}{G}.$$

$${\nabla }_r^2{u}_z(z,\mathbf{r})=2\pi R\tan (\theta (z)){\delta }^2(\mathbf{r}),$$

$$u_z(z,\mathbf{r})=-2\pi R\tan (\theta (z))C(\mathbf{r}),$$

$$U_{\text{gf}}=E_{\text{core}}+\pi Gh{R}^2\overline{{\tan }^2\theta }\ln (L/R).$$

$$k_B{T}_{\text{eff}}\sim G⟨{ϵ}_d^2⟩{\lambda }_d^2{l}_d,$$

$$⟨{ϵ}_d^2{⟩}_c=\frac{\pi }{2}\frac{R^2}{{\lambda }_d^2}\overline{{\tan }^2\theta }.$$

$$H=\sum _i{E}_{\text{core},i}+4{\pi }^{2}Gh\sum _{i<j}{R}_i{R}_j\overline{\tan {\theta }_i}\overline{\tan {\theta }_j}C({\mathbf{r}}_i-{\mathbf{r}}_j).$$

$$Z=\sum _{\{n_i\}}\int \prod _i{d}^2{r}_i{y}_0^{{\sum }_i{n}_i^2}{e}^{4\pi \ln y_0{\sum }_{i<j}{n}_i{n}_{j}C({\mathbf{r}}_i-{\mathbf{r}}_j)},$$

$${\nabla }_r^2{u}_z=2\pi s\sum _{\text{innies}}{\delta }^2({\mathbf{r}-\mathbf{r}}_i)-\alpha ,$$

$$u_z(\mathbf{r})$$

$$(\nabla u_z{)}^2$$

Full text

Cusped inward folds known as creases form on compressed surfaces of a variety of soft elastic materials Jin et al. (2014), including natural rubber Southern and Thomas (1965); Gent and Cho (1999), polymer gels Tanaka et al. (1987); Trujillo et al. (2008); Yoon et al. (2010), silicone elastomers Cai et al. (2012); Chen et al. (2012, 2014); Liang et al. (2016); Tallinen et al. (2014, 2016), starchy foods Hong et al. (2009); Cai et al. (2010), and the developing mammalian brain Tallinen et al. (2013, 2014, 2016); Fernández et al. (2016). In the latter context, creases are called “sulci”. Unlike the long-wavelength buckling of a compressed beam, or the smooth sinusoidal wrinkles observed on the skin of drying fruit or a tensioned elastic sheet Cerda and Mahadevan (2003); Huang et al. (2007); Paulsen et al. (2016), creases are sharply localized in both their elastic deformation and stresses, thereby defying a linear perturbation analysis Gent and Cho (1999); Hong et al. (2009); Hohlfeld and Mahadevan (2011, 2012); Cai et al. (2012). Owing to this difficulty, numerical minimization of a nonlinear neo-Hookean energy functional has become the standard theoretical tool for investigating the onset of creases Hong et al. (2009); Amar and Ciarletta (2010); Hohlfeld and Mahadevan (2011, 2012); Cai et al. (2012); Tallinen et al. (2013, 2014, 2016). A central claim in much of this work is that creasing is a fundamentally new, nonlinear instability with no scale Hohlfeld and Mahadevan (2011, 2012). Experimental work has also studied the growth of pre-existing long creases, describing these behaviors in analogy to crack propagation Chen et al. (2012, 2014).

Here we develop a new quasi-particle framework for shear stress focusing in surface-compressed solids, assuming planar geometry and neglecting surface tension. We apply our theory to the creasing instability, obtaining a markedly different picture than Hohlfeld and Mahadevan (2011, 2012). We find evidence that (i) creasing onset maps to the Kosterlitz-Thouless (KT) transition Kosterlitz and Thouless (1973), (ii) nonlinear elasticity is needed only within a small region analogous to a vortex core, and (iii) compression-induced shear strain fluctuations set the fundamental, microscopic lengthscale in the problem. Our theory makes contact with experimental and simulation results on critical strain, surface profiles, and crease patterns. In particular, we obtain a universal critical compressive plane strain above which creases emerge, in reasonable agreement with the measured value of $35\%$ Trujillo et al. (2008); Hong et al. (2009); Hohlfeld and Mahadevan (2012). Finally, the theory points to a set of minimal physical ingredients for creasing, and suggests a possible unification with ridging (formation of localized surface protrusions) Takei et al. (2014), and dimple crystallization Brojan et al. (2015); Jiménez et al. (2016).

Our point of departure from prior work is to consider a distinct regime of zero-length creases, qualitatively similar to those observed in Yoon et al. (2010); Cai et al. (2012); Chen et al. (2012), immediately upon nucleation, and those in Tallinen et al. (2013), as the critical point is approached from above. Foundational to our theory is the observation that zero-length creases and zero-length anticreases (hereafter, innies and outies) also appear in a very different continuum elastic context, namely the shear lag model of composite materials science and engineering Cox (1952); Hull and Clyne (1996). In this model, one assumes that shear coupling is supported at the interface between a low-dimensional reinforcing phase (i.e., 1d fibers or 2d slabs) and a surrounding 3d matrix phase. Next, an approximation is made that the transfer of axial loads between the two components is accomplished entirely via tension or compression in the reinforcing phase, and pure shear in the matrix. Axial loads refer to external or internal forces (such as those arising from differential growth of the two components) acting parallel to a long axis of the reinforcing phase. In the case of a fiber-matrix composite, the model predicts that matrix shear stress and strain fall off as $1/r$, where $r$ is the perpendicular distance from a fiber. Thus, the matrix deformation and hence the surface profile scales as $\ln r$ (see Figure 1).

Mapping shear lag to 2d electrodynamics —

Let us take all forces in the shear lag model along $z$. In the matrix phase, the only non-negligible components of the strain tensor $\epsilon_{ij}=(\partial_{i}u_{j}+\partial_{j}u_{i})/2$ have one index equal to $z$ and the other not equal to $z$. Defining a 2d vector of shear strains $\vec{\gamma}=2(\epsilon_{xz},\epsilon_{yz})$, force balance on a volume element of matrix takes the form of Gauss’s Law

[TABLE]

Here $G$ is the matrix shear modulus while $f_{0}$ is the net force supported at the interface and enclosed by the free body diagram (Gaussian surface). In terms of the wetted perimeter $p_{0}(z)$ and interfacial shear stress $\tau_{0}(z)$, the interfacial force is $f_{0}=\int dz\,p_{0}\tau_{0}$. Our neglect of $\epsilon_{zz}$ is justified by requiring $p_{0}\tau_{0}$ vary slowly with $z$.

Equivalently, one has a Poisson’s equation for the (scalar) deformation field $u_{z}$. For the important special case of a thin fiber source at the origin,

[TABLE]

where $R$ is the fiber radius and $\tan\theta=\tau_{0}/G$ is the interfacial shear strain, as indicated in Figure 1. $2\pi R\tan\theta$ is the “charge” per unit length of fiber. (Note that in the case of a viscoelastic matrix with storage and loss moduli given by $G^{\prime}$ and $G^{\prime\prime}$, respectively, a frequency-dependent “dielectric function” appears as $1+iG^{\prime\prime}/G^{\prime}$.) The solution of Equation 2 is

[TABLE]

where $C(\mathbf{r})=-\ln(|\mathbf{r}|/R)/2\pi$ is the 2d Coulomb potential (Green’s function).

Guided by this mapping, we ask whether $u_{z}$ might be associated with a quasi-charge excitation in which the fiber of the conventional shear lag model is an abstraction. For system size $L$ and vertical thickness $h$, the elastic strain energy (electrostatic energy) $\frac{G}{2}\int dV\,\gamma^{2}$ required to create an isolated such “ghost fiber” is

[TABLE]

The overbar denotes an average over $z$, and the core contribution is from deformations within $r<R$. Notably, this result resembles the energy of an isolated vortex Kardar (2007), and $\frac{G}{2}\gamma^{2}$ is the leading order term of the energy density used in the aforementioned neo-Hookean simulations imposing an incompressibility constraint. Motivated thus, we now turn to the statistics of ghost fibers.

Effective thermodynamics of shear lag quasiparticles —

Creasing in a flat geometry requires some spatial inhomogeneity to break translation invariance and nucleate the creases. Typically there is disorder at both the mesoscale and microscale, due to dust, surface imperfections, and the polymer network itself. In experiments, nucleation sites have also been introduced in a controlled manner via surface-embedded microspheres Yoon et al. (2010). Creasing simulations rely upon defect seeding Cai et al. (2012), inhomogeneous surface normal forces Hohlfeld and Mahadevan (2012), random vertical displacements Tallinen et al. (2013), or random variations in growth Tallinen et al. (2014). In all such cases, these defects pattern the surface with stress concentrators that give rise to strain fluctuations within a thin subsurface layer, when the system is macroscopically compressed.

The central postulate of this work is that disorder and drive-induced shear strain fluctuations (from one or more sources) can play the role of thermal fluctuations, i.e. there is an “effective temperature” that increases during a slow compression of a disordered elastomer or gel, enabling a strongly athermal system of ghost fibers to access many microstates (position and charge configurations). Similar approaches have been used to describe the statistics of granular materials, pinned vortex lattices, artificial spin ice, and other athermal systems Bi et al. (2015); Cugliandolo (2011); Casas-Vázquez and Jou (2003); Nisoli et al. (2010); Kolton et al. (2002). Here, our effective temperature postulate is of the form

[TABLE]

where $\langle\epsilon_{d}^{2}\rangle$ is the mean squared amplitude of $\epsilon_{xz}$ and $\epsilon_{yz}$ fluctuations, and $\lambda_{d}$, $l_{d}$ are the characteristic wavelength and skin depth of these fluctuations. Equation 5 is reminiscent of the Lindemann criterion for bulk melting of a harmonic solid Lindemann (1910). In the following, we will identify the effective system thickness $h$ with $l_{d}$, below which the system is cold and inactive.

Working within the microcanonical ensemble, the configurational entropy of a system containing a single ghost fiber is $S=2k_{B}\ln(L/R)$. The free energy cost to create the ghost fiber is $F=U_{\textrm{gf}}-T_{\textrm{eff}}S$. In the thermodynamic limit where the finite core energy is dominated by the logarithmically divergent term, $F<0$ for mean square strain fluctuations greater than the critical value

[TABLE]

Because the quantity that maps to electric charge is an odd function of $\theta$, a charge dipole corresponds to an innie-outie pair of surface deformations, equivalent to a tension-compression pair of ghost fibers (see Figure 1). The energy required to create a pure dipole is finite, in contrast to Equation 4, and thus 3d ghost fibers are analogous to 2d vortices: for $\langle\epsilon_{d}^{2}\rangle<\langle\epsilon_{d}^{2}\rangle_{c}$, the system contains tension-compression bound pairs of ghost fibers, and at $\langle\epsilon_{d}^{2}\rangle_{c}$ there is an unbinding transition (KT transition).

We now consider the grand canonical ensemble. A charge neutral system of ghost fibers has Hamiltonian

[TABLE]

In the regime where elastic deformations within the core are linear, or weakly nonlinear, or there are strong nonlinearities but they are confined to $r\ll R$, simple scaling arguments indicate $E_{\textrm{core}}\sim GV_{\textrm{core}}\epsilon_{\textrm{core}}^{2}\sim Gh\pi R^{2}(\overline{\tan\theta})^{2}$. This quantity would appear to vary from one quasi-particle to another because the charges $\sim R_{i}\overline{\tan\theta_{i}}$ are here continuous degrees of freedom. However, we can exploit the arbitrariness of the $R_{i}$ in order to take the core energy as a meaningful chemical potential $\mu$. The appropriate choice is $R_{i}=s|\overline{\tan\theta_{i}}|^{-1}$, where $s\sim\sqrt{\mu/(Gh\pi)}$. This brings the partition function into the Coulomb gas form

[TABLE]

where $n_{i}=\pm 1$ and $y_{0}=\exp[-\mu/k_{B}T_{\textrm{eff}}]$ is the ghost fiber fugacity. The price paid for replacing continuous charges with discrete ones is that $R$ is now an ambiguous “lattice constant”. However, we have made available a small and well-defined lengthscale $s$; this can presumably replace $R$ as the short distance cutoff.

The fugacity $y_{0}$ and “coupling constant” $K=-(\ln y_{0})/\pi$ are related because we are considering a specific physical system (e.g. Kardar (2007)). Intersection of the line $y_{0}=e^{-\pi K}$ and the line of fixed points $y_{0}=-\pi^{-2}(K^{-1}-\pi/2)$ determines the critical inverse coupling $K_{c}^{-1}=1.06$, and hence the critical mean square strain fluctuation $\langle\epsilon_{d}^{2}\rangle_{c}=(s^{2}/\lambda_{d}^{2})K_{c}^{-1}$. Note $K_{c}^{-1}$ is depressed from the mean field value $\pi/2$ obtained earlier.

Below $\langle\epsilon_{d}^{2}\rangle_{c}$, large-scale surface deformation would not be seen because the tension-compression pairs are tightly bound. The appearance of spatially separated, cusped surface deformations at a critical point that has no explicit dependence on system thickness or shear modulus is consistent with creasing experiments Trujillo et al. (2008). So too, we argue from data in Cai et al. (2012); Chen et al. (2012); Yoon et al. (2010); Tallinen et al. (2013), is the notion of zero-length creases at the critical point. (But see Chen et al. (2014) for a different interpretation of crease lengths.) In light of the apparently universal 35% onset strain Trujillo et al. (2008); Hong et al. (2009); Hohlfeld and Mahadevan (2012), one is tempted to identify $\lambda_{d}$ with the lattice constant $s$ (possibly scaled by a numerical prefactor). In fact, the simulations of Tallinen, et al. use random vertical displacements of mesh surface nodes, consistent with this picture Tallinen et al. (2013). Their method of introducing fluctuations suggests a way to estimate the critical compressive plane strain $\epsilon_{c}$, given $\langle\epsilon_{d}^{2}\rangle_{c}$, but first we need a brief digression.

A consequence of all the stresses in the problem being shear stresses is that the area of the free surface is approximately conserved (within linear elasticity). One can see this by considering a 2d finite element in simple shear, and noting that its perimeter change is a second order effect. Here and in one place to follow, we shall make use of this approximate “area conservation principle”.

Returning to the problem at hand, we approximate the shear strain fluctuations by a square wave with amplitude $\sqrt{\langle\epsilon_{d}^{2}\rangle}$, such that the corresponding fluctuations in vertical displacement are a triangular wave with amplitude $\frac{1}{2}\lambda_{d}\sqrt{\langle\epsilon_{d}^{2}\rangle}$. Setting $s=\lambda_{d}/2$, as suggested by the method in Tallinen et al. (2013), one then finds $\epsilon_{c}=1-(1+K_{c}^{-1})^{-1/2}=38\%$ (mean field theory) and $=30\%$ (renormalization group), which bracket 35%.

Post-KT dynamics of ghost fibers —

Two features of creasing experiments and simulations remain to be explained by our quasi-particle theory: (i) that only innies and not outies appear to be seen, and (ii) that innies smoothly become finite-length creases. In this section we consider (i), and in the next section we will consider (ii).

The KT transition does not involve (or at least, does not require) self-contact in the core region. Yet self-contact is generically observed Hong et al. (2009); Hohlfeld and Mahadevan (2012); Chen et al. (2012); Cai et al. (2012); Tallinen et al. (2013). We propose that self-contact ensues at strain $\epsilon_{sc}>\epsilon_{c}$, and point out that it can only be available to innies, because a self-contacting outie is an unphysical concept. Appealing to area conservation (now only qualitatively useful since the core deformations are strongly nonlinear), outies should incur a higher energy penalty than innies, in the regime $\epsilon>\epsilon_{sc}$, because they are not as effective at sequestering surface area. The system cannot exactly get rid of its outies, however. Doing so would generate a nonzero net charge, causing $u_{z}$ to grow with system size as $(\sum_{i}n_{i})\ln L$, clearly inconsistent with creasing experiments. What happens, we propose, is that an outie’s $R$ increases while its $|\overline{\tan\theta}|$ decreases, in such a way that its charge $-2\pi R|\overline{\tan\theta}|$ stays fixed. In other words, the negative point charges get smeared out into a negative background charge (cf. the non-neutral Coulomb gas Minnhagen (1987)). The Poisson equation for this situation reads

[TABLE]

where $-\alpha$ is a uniform negative charge density, interpreted as the surface curvature the system would have, if the innies were removed.

One might ask if this innie-outie shape asymmetry could be present even during the KT transition. We suggest the answer is no, because the ghost fiber bound pairs that exist below $\epsilon_{c}$ must have an essentially net zero surface deformation in order to be consistent with the observed flat surface. The shape asymmetry is something that arises in connection with the energy penalty for unable-to-self-contact outies. One might also ask if there exist conditions in which a system of outies in a neutralizing positive background is realized. Tentative support for this idea comes from recently observed ridges Takei et al. (2014), which are somewhat reminiscent of anticreases.

Ghost slabs —

Having considered ghost fibers, we now consider a thin “ghost slab” of length $2\ell$ and height $h$ characterized by charge density $\rho\sim\delta(x)[H(y+\ell)-H(y-\ell)]$. (The case $\alpha\neq 0$ will be treated momentarily). A straightforward application of Green’s method yields material deformation

[TABLE]

and strain energy density

[TABLE]

Equations 10 and 11 are plotted in Figure 2 for $\ell=1/2$. In the thermodynamic limit, the rotational contribution to the ghost slab entropy is insignificant, and the only important contribution to the elastic energy comes from the monopole term. Similar arguments to those used above lead to a critical strain fluctuation for slabs $\langle\epsilon_{d}^{2}\rangle_{c,\textrm{ slab}}\sim(\ell^{2}/\lambda_{d}^{2})\,\overline{\tan^{2}\theta}$. This result reveals ghost fibers to be a limiting case of ghost slabs (i.e., as $\ell\to R$).

Experiments and simulations employing uniaxial strain tend to generate straight, parallel creases; when the strain is large, the creases can span the system size, forming a 1d array Cai et al. (2012); Chen et al. (2012); Tallinen et al. (2013). Within our electrostatics framework, for $\alpha\neq 0$, surface profiles of such systems are predicted to have the same form as the electrostatic potential of a 1d Coulomb crystal (a periodic stack of infinite, charged sheets embedded in a charge-neutralizing, uniform background). Figure 3a shows a test of this prediction, with no adjustable fit parameters apart from an overall prefactor. The possibility of the background being “polarizable” is investigated via the $d$-dimensional charge screening equation $(\nabla^{2}_{d}-\lambda^{-2})\phi=-4\pi Q\delta^{(d)}$, where $Q$ is charge and $\lambda$ is the screening length Johnson-McDaniel and Owen (2012). The $d$=1 solution, $\phi(x)\sim e^{-x/\lambda}$, is fit to the profile of an isolated long crease in Figure 3a inset. The qualitative shape of the crease array in Figure 3a (parabolic crests between sharp cusps) occurs in other settings such as a mouse cerebellum and a winter squash (Figure 3b,c). Prior work has modeled these as elastic materials Lejeune et al. (2016); Hu et al. (2011), suggesting the same mechanism may be at play, in spite of the different (curved versus planar) geometry.

Experiments and simulations employing equibiaxial strain tend to generate a square lattice of short, straight creases with each nearest neighbor pair having relative orientation of 90 degrees; at high strains, a hexagonal lattice of 3-fold symmetric, Y-shaped creases is also seen Tallinen et al. (2013). (Figure 3c in Trujillo et al. (2008) exhibits both motifs.) Within our framework, such patterns again have a natural interpretation as charge-crystallized shear lag quasi-particles. In the limit where the spatial extent of quasi-particles is very small compared with their lattice spacing (i.e., they are ghost fiber-like), the creasing pattern is predicted to be a 2d Coulomb crystal with hexagonal symmetry (e.g. Bonitz et al. (2008)). In fact the hexagonal dimple crystal observed in Brojan et al. (2015); Jiménez et al. (2016) is suggestive of this, although its discoverers give a wrinkling, not creasing, interpretation. Outside the limit of fiber-like quasiparticles, the ground state crystal structure is an interesting topic for future work.

Discussion —

Upon building a composite materials-inspired quasi-particle framework, we have uncovered a creasing scenario that involves at least three distinct regimes. For in-plane compressive strain $\epsilon<\epsilon_{c}$, the system contains tightly bound pairs of ghost fibers whose deformation fields largely cancel. For $\epsilon_{c}<\epsilon<\epsilon_{sc}$, the pairs are unbound, giving rise to spatially separated, cusped surface deformations. Also in this regime, innies (outies) smoothly become finite length creases (anticreases), as ghost fibers smoothly change dimensionality into ghost slabs. For $\epsilon>\epsilon_{sc}$, anticreases smear out into a charge-compensating background, while repulsive interactions between creases causes them to organize into a Coulomb crystal. The KT transition is predicted to occur at a universal critical strain $\epsilon_{c}=30\%$ (renormalization group analysis) and self-contact is expected to commence at a slightly higher strain $\epsilon_{sc}$, which we speculate may be the previously measured onset strain, $35$%. Both transition points should be observable, by virtue of the characteristics of the regimes they delineate (e.g. surface profiles), providing a means for experimental tests of our theory.

Acknowledgements.

The authors are grateful to Ryan Hayward and Joseph Paulsen for many helpful discussions, and acknowledge financial support from NSF-DMR-CMMT Award Number 1507938.