††thanks: These two authors contributed equally††thanks: These two authors
contributed equally††thanks: Corresponding author<EMAIL_ADDRESS>
# Geometric localization in supported elastic struts

T.C.T. Michaels  Paulson School of Engineering and Applied Sciences, Harvard
University, Cambridge, MA 02138, USA  R. Kusters  University Paris Descartes,
Center for Research and Interdisciplinarity (CRI), 10 Rue Charles V, Paris
Department of Applied Physics, Eindhoven University of Technology,
Eindhoven, The Netherlands  A.J. Dear  Paulson School of Engineering and
Applied Sciences, Harvard University, Cambridge, MA 02138, USA  C. Storm
Department of Applied Physics, Eindhoven University of Technology,
Eindhoven, The Netherlands  Institute for Complex Molecular Systems, Eindhoven
University of Technology, Eindhoven, The Netherlands  J.C. Weaver  Paulson
School of Engineering and Applied Sciences, Harvard University, Cambridge, MA
02138, USA  L. Mahadevan  Paulson School of Engineering and Applied Sciences,
Harvard University, Cambridge, MA 02138, USA

###### Abstract

Localized deformation patterns are a common motif in morphogenesis and are
increasingly finding widespread applications in materials science, for
instance as memory devices. Here we describe the emergence of spatially
localized deformations in a minimal mechanical system by exploring the impact
of growth and shear on the conformation of a semi-flexible filament connected
to a pliable shearable substrate. We combine numerical simulations of a
discrete rod model with theoretical analysis of the differential equations
recovered in the continuum limit to quantify (in the form of scaling laws) how
geometry, mechanics, and growth act together to give rise to such localized
structures in this system. We find that spatially localized deformations along
the filament emerge for intermediate shear modulus and increasing growth.
Finally, we use experiments on a 3D printed multi-material model system to
demonstrate that external control of the amount of shear and growth regulates
the spatial extent of the localized strain texture.

††preprint: APS/123-QED

The deformation of spatially extended elastic structures, such as filaments,
plates and shells is often elastically constrained by the surrounding medium.
The simplest formulation of this elastic constraint is due to Winkler
Winkler1867  who proposed a local linear elastic model for the medium. Since
then, there have been many variants that account for both nonlocality and
nonlinearity, particularly in the context of localized deformations in such
systems that take the form of creases, localized wrinkles etc. Cao2012 ;
Jin2011 ; Takei2014 ; Hunt1998 ; Hunt2000 ; Pocivavsek2008 ; Yoo2008 ;
Diamant2011 ; Brau2013 ; Auguste2014 ; Jin2015 . These localized patterns in
spatially extended dynamical systems are of intrinsic interest due to their
potential role as mediators of turbulence in fluids, optics and beyond
Knobloch2015 ; Pomeau_1986 ; Coullet_2000 ; Coullet_2003 . In the context of
elasticity, a recent proposal suggests the use of localized dimples as
programmable elastic bits (e-bits) to store memory on a featureless elastic
shell that requires minimal substrate infrastructure Chung_2018 , and raises
the question of whether there are alternatives that take advantage of
substrate elasticity to confine deformations and thus generalize these ideas
to a much broader class of systems.

In this paper, we address this question by studying the appearance of
localized structures in a minimal model system of an extensible growing
adherent filament that is attached to a rigid substrate by a set of
stretchable and shearable springs. Using a combination of theory, simulations,
and experiments, we show that this minimal system can give rise to robust
spatial localization as a function of growth and shear. More specifically, we
derive explicit scaling relationships describing the extent and amplitude of
localized deformations, as well as a phase diagram describing the range of
mechanical parameters that admit localized structures. We then demonstrate
these rules in practice using a 3D printed multi-material model system. Our
results identify the physical parameters that control the extent of
localization, suggesting a strategy for easy and robust programming mechanical
deformations at multiple scales.

Figure 1: Schematic representation of the discrete rod model of geometric
localization. A flexible filament (with stretching stiffness $\bar{S}$ and
bending stiffness $\bar{B}$) is attached to a rigid substrate through an array
of springs that can extend vertically (with spring stiffness $\bar{K}$) and
shear horizontally (with shear modulus $\bar{G}$). The the flexible filament
undergoes uniform growth, such that its rest length $L$ increases to
$L(1+\epsilon)$, where $\epsilon$ is the growth strain. Depending on shear
modulus $\bar{G}$, filament growth $\epsilon$ can result in the formation of
localized deformation patterns. Understanding the conditions for the emergence
of such localized deformations as a function of mechanical parameters is the
subject of this paper. Figure 2: Computer simulations of discrete rod model.
(a-b) Depending on the shear modulus of the springs, filament growth results
in uniform buckling for high shear modulus $\bar{G}=0.01$ (a), and localized
deformations for low shear modulus $\bar{G}=0.0002$ (b). In both cases
$N=100$, $n_{s}=316$, $b=1$, $h_{0}=1$, $\bar{K}=0.1$, $\bar{S}=100$,
$\bar{B}=0.4$, and $\epsilon=3.8\%$. (a) and (b) also show the curvature
($W^{\prime\prime}\times 10^{5}$, top purple line) and strain ($U^{\prime}$,
bottom blue line) fields along the filament. (c-d) Filament deformations for
increasing growth strain $\epsilon$ (c), and decreasing shear modulus for the
springs $\bar{G}$ (d). In (c) growth strain is $\epsilon=0.0004,$ $0.012$,
$0.059$, and $0.075$ (bottom to top) at constant shear rate $\bar{G}=0.002$.
In (d) the shear modulus is $\bar{G}=0.05,$ $0.01$, $0.005$, and $0.002$
(bottom to top) for $\epsilon=0.075$. In both (c) and (d) $N=50$.

## Theory

Our physical model of geometric localization is a flexible filament that can
swell or grow, but is constrained by an elastic substrate to which it is
attached. The filament can bend transversely and stretch longitudinally, while
the substrate resists lateral shear and transverse stretch induced by the
filament. We consider two limits of this problem: (i) a discrete version,
where the filament is assumed to be made of a set of rigid links with soft
connectors that resist bending, connected to a substrate via a set of discrete
springs that resist changes in their natural length and orientation relative
to the rigid rods (Fig. 1), and (ii) a continuum version, which maps onto the
growth of a filament embedded in an elastic medium. In either case, we seek
the solution with minimal total elastic energy for a filament with initial
rest length $L$ which has grown/swollen to a total length of $L(1+\epsilon)$,
where $\epsilon$ measures the relative growth (i.e. strain) of the filament.

### Discrete rod model

The discrete rod model consists of a flexible filament made of $n_{s}$
discrete segments of rest length $\ell_{0}$ attached to a solid substrate
through a set of $N$ Hookean springs that can extend vertically and shear
horizontally (Fig. 1). The total energy of the discrete rod model is

$\displaystyle\mathcal{E}$
$\displaystyle=\frac{\bar{S}}{2\ell_{0}}\sum_{n=1}^{n_{s}}\left(\ell_{n}-\ell_{0}\right)^{2}+\frac{\bar{B}}{\ell_{0}}\sum_{n=1}^{n_{s}-1}\left(1-\cos\theta_{n,n+1}\right)$
(1)
$\displaystyle\quad+\frac{\bar{K}}{2b}\sum_{i=1}^{N}\left(h_{i}-h_{0}\right)^{2}+\frac{\bar{G}}{2}\sum_{i=1}^{N}\alpha_{i}^{2},$

where $\bar{S}$ denotes the filament stretching stiffness, $\bar{B}$ is the
bending rigidity of the filament, $\bar{K}$ is the spring stiffness, and
$\bar{G}$ is the shear modulus of the springs. Moreover, $\ell_{n}$ is the
length of the $n$-th filament segment, $b=L/N$ is the spacing between the
adherent springs, where $L=n_{s}\ell_{0}$ is the total length of the
supporting filament. In (1), the first two terms correspond to the stretching
and bending energies of the filament, while the last two terms correspond to
the stretching and shear energies of the springs connecting the filament to
the rigid substrate. The bending energy is consistent with the semi-flexible
approximation valid for weak bending, where
$\cos\theta_{i,i+1}=\hat{t}_{i}\cdot\hat{t}_{i+1}$, $i=1,\dots,n_{s}$, where
$\hat{t}_{i}$ is the unit tangent vector of the $i$-th filament segment.
$\alpha_{i}$ denotes the angle that the $i$-th spring makes with the vertical,
hence measuring the amount of shear; $h_{i}$ is the vertical extension of the
$i$-th spring and $h_{0}$ is the rest length of the springs.

Figure 3: Theoretical analysis of geometric localization. (a) Vertical
displacement field $W(x)$ obtained by numerical integration of (3) (solid
line) and comparison to the analytical perturbation solution, (8) (black
dotted line). Parameters: $G=5\times 10^{-5}$, $S=10$, $B=0.2$, $K=0.01$,
$\epsilon=9.5\times 10^{-3}$, $L=23\times\lambda$. (b) Bifurcation diagram for
(3) as a function of growth strain $\epsilon$. For $\epsilon<\epsilon^{\star}$
there is no deformation, while for $\epsilon>\epsilon^{\star}$ either uniform
or localized deformations emerge. For $\epsilon>\epsilon^{\star}$, the
amplitude of deformations increases in proportion to
$\sqrt{\epsilon-\epsilon^{\star}}$, as predicted by (8) (solid line). The
figure also shows the width $w$ and wavelength $\lambda$ of localized
deformations as a function of $\epsilon$, with solid lines indicating the
theoretical predictions. Data points are from numerical simulations of the
discrete rod model for the following parameters: $N=80$, $\bar{S}=100$, $b=1$,
$n_{s}=316$, $\bar{K}=0.1$, $\bar{B}=0.4$ and $\bar{G}=0.001$. (c) Double
logarithmic plots of amplitude $A$ and width $w$ of localized solutions as
function of shear modulus $G$. The solid lines, which have slopes of 1/4,
respectively, -1/4, are the theoretical predictions of (10) and (9), $A\propto
G^{1/4}$, $w\propto G^{-1/4}$. Simulation parameters: $S=10$, $B=1$, $K=1$,
$\epsilon=0.21$, $L=20\times\lambda$. (d-e) Uniform or localized deformations
are selected depending on their relative energies, as described by
$(\mathcal{E}_{\rm{loc}}-\mathcal{E}_{\rm{uni}})/\mathcal{E}_{\rm{uni}}$,
where $\mathcal{E}_{\rm{loc}}$ is the energy of the localized state and
$\mathcal{E}_{\rm{uni}}$ is the energy of the uniform state. Panel (d) shows
that, at constant $G$, there is a critical value for $\epsilon$ above which
the localized solution has lower energy compared to the uniform buckled
solution. The critical $\epsilon$ increases with increasing $G$. Panel (e)
shows that, at constant $\epsilon$, there is a critical shear modulus $G$,
below which the localized solution is lower in energy than the uniform
deformation. The critical $G$ increases with increasing $\epsilon$. Simulation
parameters for (d-e) are: $N=50$, $\bar{S}=100$, $b=1$, $n_{s}=197$,
$\bar{K}=0.1$, $\bar{B}=0.4$, $\bar{G}=0.001,0.0025$, and $\epsilon=0.1,0.13$
as indicated on the graph.

The procedure we use to “grow” the filament is as follows: we initiate a
configuration where the springs are in their rest configuration and vertical.
We then increase the rest length of the filament segments $\ell_{0}$ by a
small amount $\epsilon\ell_{0}$, i.e. $\ell_{0}\to\ell_{0}(1+\epsilon)$. Since
thin rods are much easier to bend than to stretch, we focus our analysis on
filaments with $\bar{S}\gg\bar{K}$; some extensibility is retained for
increased numerical stability. To minimize the total energy $\mathcal{E}$ of
the system we either use Newton’s method, present in the FindMinimum routine,
or a global energy search present in the NMinimize routine. Both routines are
implemented in Wolfram Mathematica (see SI Sec. S2).

### Continuous model

To determine how localized filament deformations emerge during growth, we map
the discrete rod model onto a set of coupled differential equations that
describe the growth of the adhering filament in the continuum limit, valid for
small deformations. To this end, we consider the elastic energy (1) in the
limits $b,\ell_{0}\to 0$ and $N,n_{s}\to\infty$ with $Nb=n_{s}\ell_{0}=L$
constant. In this limit, the total elastic energy ${\cal E}$ of the system is
written as a function continuous horizontal $U(x)$ and vertical $W(x)$
displacements, where $x$ is the arclength of the filament, as (see SI Sec.
S1):

$\displaystyle{\cal E}$
$\displaystyle=\int_{L(1+\epsilon)}\left(\frac{S}{2}\left[U^{\prime}(x)+\frac{W^{\prime}(x)^{2}}{2}-\epsilon\right]^{2}\right.$
$\displaystyle\qquad\left.+\frac{B}{2}W^{\prime\prime}(x)^{2}+\frac{K}{2}W(x)^{2}+\frac{G}{2}U(x)^{2}\right)dx.$
(2)

The elastic constants in the continuum model relate to those in the discrete
model (denoted with bar) through $G=\bar{G}/(h_{0}^{2}b)$, $K=\bar{K}/b^{2}$,
$S=\bar{S}b/{\ell_{0}}$ and $B=\bar{B}/h_{0}$ (see SI Sec. S1). The various
energy terms in (2) have a straight-forward physical interpretation: the first
term is the stretching energy of the filament and is proportional to the
square of the elastic strain
$\epsilon_{\rm{el}}=\epsilon_{\rm{tot}}-\epsilon$, where $\epsilon_{\rm{tot}}$
is the total strain and $\epsilon$ is the growth strain. The second term in
(2), proportional to the square of the local curvature $W^{\prime\prime}(x)$,
describes the bending energy of the filament. The last two terms in (2)
correspond to the energy contributions associated with the stretching and the
shear of the substrate; these contributions are proportional to the squares of
the vertical and horizontal displacement fields, respectively.

To find the system configurations of minimal total energy, we consider the
Euler-Lagrange equations associated with (2), $\delta{\cal E}/\delta
W(x)=\delta{\cal E}/\delta U(x)=0$, which read (see SI Sec. S1.6):

$\displaystyle
BW^{\prime\prime\prime\prime}(x)+S\left[\epsilon-U^{\prime}(x)-\frac{W^{\prime}(x)^{2}}{2}\right]W^{\prime\prime}(x)$
$\displaystyle+$ (3a) $\displaystyle+KW(x)-GU(x)W^{\prime}(x)$
$\displaystyle=$ $\displaystyle 0$ $\displaystyle
S\left[U^{\prime\prime}(x)+W^{\prime}(x)W^{\prime\prime}(x)\right]-GU(x)$
$\displaystyle=$ $\displaystyle 0.$ (3b)

To compare with analytic results, we have used a shooting algorithm to
numerically solve (3) subject to clamped boundary conditions
$U(-L/2)=U(L/2)=0,\ W(-L/2)=W(L/2)=0,\ W^{\prime}(-L/2)=W(L/2)=0$ (see SI Sec.
S2). Note that these equations decouple and are fully integrable in two
limits: (i) for $G\to\infty$, and (ii) for $G=0$. Both situations lead to a
uniform wrinkling solution with characteristic wavelength
$\lambda=(B/K)^{1/4}$ (see SI Secs. S3.2 and S3.3).

## Analysis

### Lowering shear modulus induces localization

We have simulated the discrete rod model under different conditions of growth
$\epsilon$ and shear modulus $\bar{G}$ for clamped and free filament edges.
For clamped edges, increasing $\epsilon$ above a critical threshold
$\epsilon^{\star}$ and for large shear modulus $\bar{G}$, we observe a
uniformly periodic deformation of the filament (Fig. 2(a)). This situation is
consistent with the classical Winkler foundation Winkler1867 , corresponding
to the $\bar{G}\to\infty$ limit where springs can extend only vertically. If
we decrease $\bar{G}$, allowing for the adherent springs to hinge, we find
that, upon uniform growth $\epsilon$, the vertical displacement field along
the filament is not uniformly periodic, but rather becomes localized (Fig.
2(b)). Since the filament itself is effectively inextensible
($\bar{S}/\bar{K}=10^{3}$), this localized strain solution is accompanied by a
shear field of the adherent springs (horizontal displacement). While the
buckling wavelength remains unaltered, we find that localization is favored by
increasing $\epsilon$ (Fig. 2(c)) or decreasing $\bar{G}$ (Fig. 2(d)).
Decreasing the shear modulus further leads to an increase of the localization
width. Eventually, for vanishing shear modulus $\bar{G}=0$, we recover the
uniformly deformed solution, when the width of the localized deformation
approaches the system size $L$.

Our simulations also point to an important role of the boundary conditions:
releasing the edges results in the occurrence of localized structures that
occur near the filament edge and upon further growth travel inward. In SI Sec.
S2 we discuss this aspect further.

### Localized deformations as solutions to the continuous equations

To rationalize the emergence of spatial localized deformations in terms of
specific combinations of the underlying physical parameters, we set out to
solve the continuous equations, (3), analytically in the limit of low shear
modulus $G$ by employing weakly nonlinear analysis Drazin_2002 . We construct
this solution on the half-axis $0\leq x\leq L/2$ by employing the mirror-
symmetry of the vertical displacement around the midpoint $x=0$ (see SI Sec.
S3).

Figure 4: Phase diagram as a function of shear modulus $G/K$ and growth strain
$\epsilon$ separating regions where uniform or localized buckling deformations
are minimum energy solutions. (a) In the purple region, bounded by the
theoretical lines $(G/K)^{\rm{crit}}$ and $(G/K)^{\star}$, localized
deformations exist and have lower energy than uniform solutions. Diagram
calculated for the same parameters as in Fig. 3(d). (b) Phase diagram in high
$G$ regime obtained from numerical simulations of discrete rod model for the
parameters of Fig. 3(d) and comparison with theoretical prediction (solid
line). Circles indicate localized deformations, triangles indicate uniform
deformations. At the crosses we detected bistability. (c) Phase diagram in low
$G$ regime obtained from numerical simulations of continuum model for the
parameters of Fig. S5(b) and comparison to the theoretical prediction (12)
(solid line, which has slope -1).

#### Linear theory and bifurcation analysis

As a first step, we study the linearized equations
$BW^{\prime\prime\prime\prime}(x)+S\epsilon W^{\prime\prime}(x)+KW(x)=0$,
$SU^{\prime\prime}(x)=GU(x)$ subject to the boundary conditions associated
with clamped edges (see SI Sec. S3.1). For $\epsilon\leq\epsilon^{\star}$,
where

$\epsilon^{\star}=\frac{2\sqrt{BK}}{S},$ (4)

we find that the only solution consistent with the boundary conditions is
$W(x)=U(x)\equiv 0$, i.e. no deformation. However, for
$\epsilon\geq\epsilon^{\star}$, we find a non-trivial, uniformly periodic
wrinkling solution for the vertical displacement; the wavelength of wrinkling
is $\lambda=(B/K)^{1/4}$, consistent with previous reports Landau ; Cerda_2003
. Thus, analysis of the linearized equations implies a bifurcation point when
$\epsilon=\epsilon^{\star}$ (Fig. 3(b)). The physical interpretation of
$\epsilon^{\star}$ is the maximal amount of strain that can be accommodated by
filament stretching (compression) without causing any bending and spring
deformation. Indeed, $\epsilon^{\star}$ is inversely proportional to $S$;
decreasing $S$ leads to a higher strain threshold $\epsilon^{\star}$. Above
the critical swelling/growth strain $\epsilon^{\star}$, the excess strain
$\epsilon-\epsilon^{*}$ is accommodated by filament buckling.

#### Weakly nonlinear theory

The linear theory implies that localized deformations exist for (3) only for
$\epsilon>\epsilon^{\star}$; we thus consider the non-linear terms in (3) as a
perturbation of the linearized equations in excess strain
$\epsilon-\epsilon^{\star}$. To do so, we rescale variables as
$\bar{x}=x/\lambda$, $W=\sqrt{\epsilon-\epsilon^{\star}}\,\lambda\,\bar{W}$
and $U(x)=(\epsilon-\epsilon^{\star})\,\lambda\,\bar{U}$, solve the resulting
equation for $U(x)$ exploting $(G/S)\sqrt{B/K}\ll 1$ and insert this solution
back into (3). This procedure (detailed in SI Sec. S3.5) yields at leading
order:

$\displaystyle\bar{W}^{\prime\prime\prime\prime}(\bar{x})+$ $\displaystyle
2\bar{W}^{\prime\prime}(\bar{x})+\bar{W}(\bar{x})-\delta\frac{(\bar{x}-\bar{L}/2)^{2}}{2}\bar{W}^{\prime\prime}(x)=0,$
(5)

where the relevant perturbation parameter of the problem emerges as
$\delta=(G/K)(\epsilon-\epsilon^{\star})\ll 1$. In (5) we recognize the linear
part, which yields uniform buckling with wavelength $\lambda$; as expected, we
obtain a new term, proportional to $\delta$, which comes from the excess
strain $\epsilon-\epsilon^{\star}$ and varies slowly on distances of scale
$\lambda$. This problem is analogous to the WKB approximation to the
Schrödinger equation with a slowly varying potential. The form of (5) thus
suggests the following envelope-type Ansatz for the vertical displacement

$\bar{W}(x)=A(\delta^{1/4}\bar{x})e^{i\bar{x}}+\textrm{c.c.}$ (6)

where $A$ is a slow-varying amplitude (envelope) function that depends on the
slow variable $X=\delta^{1/4}\bar{x}$ and c.c. stands for complex conjugate.
Inserting this Ansatz into (5) and collecting terms at various orders in
$\delta^{1/4}$, we arrive at the following amplitude equation describing the
long-scale behavior of the solution (see SI Sec. S3.6):

$\frac{\partial^{2}A(X)}{\partial X^{2}}=\frac{(X-X_{e})^{2}}{8}A(X),$ (7)

where $X_{e}=\delta^{1/4}L/\lambda$. The amplitude equation is a particular
case of the Weber differential equation
$y^{\prime\prime}(z)+\left(\nu+1/2-z^{2}/4\right)y(z)=0$ with $\nu=-1/2$,
whose solution is expressed in terms of the parabolic cylinder function
$y(z)=D_{-1/2}(z)$ Abramowitz . Hence, the final perturbation solution for the
vertical displacement is given by:

$\displaystyle W(x)$
$\displaystyle\propto\sqrt{\epsilon-\epsilon^{\star}}\lambda\,\cos\left(\frac{x}{\lambda}\right)\,\times$
(8) $\displaystyle\quad
D_{-1/2}\left[(\epsilon-\epsilon^{\star})^{1/4}\left(\frac{G}{2B}\right)^{1/4}\left(x-\frac{L}{2}\right)\right].$

The solution has the form of the uniform wrinkled solution $\cos(x/\lambda)$
modulated by an envelope function $A(X)$, which depends on the slow variable
$X=(\epsilon-\epsilon^{\star})^{1/4}(G/B)^{1/4}x$. Note that, independently of
growth strain $\epsilon$ or shear modulus $G$, the solution always selects the
wavelength $\lambda$. The accuracy of (8) against numerical integration of (3)
is shown in Fig. 3(a) and in SI Figs. S7-S9.

Figure 5: Experimental realization of geometric localization using 3D
printing. (a) Schematic representation of the experimental realization of the
system using a multi-material 3D printed plastic model for the growth of a
filament attached to a shearable substrate. The ends of the flexible filament
are clamped. (b-c) Images of swelling 3D-printed samples placed in isopropyl
alcohol showing the development of localized deformations over time: the
number of vertical springs is $n_{s}=28$ (b), and $n_{s}=56$ (c). The
different spacing between “springs” allows us to control the mechanical
parameters $G,K,S$, hence the form of localized deformations. Scale bar is 1
cm. (d) Measured amplitude of deformations for the samples in (b) and (c) as a
function of growth strain $\epsilon-\epsilon^{\star}$ is compared to theory,
(10) (solid line).

#### Scaling laws for localized solution

With a perturbative solution to the continuous equations (3) at hand, we can
now predict, from first principles, the scaling behavior of a number of key
observables, such as the width $w$, the wavelength $\lambda$, and the
amplitude $A$ of the localized deformations. The resulting scaling
relationships provide a useful means for designing localized deformations
across different scales. In particular, from (8), the width $w$ of the
localization scales as:

$\frac{w}{\lambda}\propto\delta^{-1/4}=(\epsilon-\epsilon^{\star})^{-1/4}\left(\frac{G}{K}\right)^{-1/4}.$
(9)

Hence, increasing growth strain $\epsilon$ beyond the bifurcation point
$\epsilon^{\star}$, or increasing $G$ favors localization (Fig. 3(b,c)). The
wavelength of the localized deformation $\lambda=\left(B/K\right)^{1/4}$ is
independent of growth strain $\epsilon$ and shear modulus $G$ (Fig. 3(b)).
Finally, the amplitude $A$ of the localization scales as

$\frac{A}{\lambda}\propto(\epsilon-\epsilon^{\star})^{1/2}\left(\frac{G}{K}\right)^{1/4}\left(\frac{L}{\lambda}\right)^{1/2}.$
(10)

We have tested these scaling predictions using numerical realizations of the
continuous equations obtained with a shooting algorithm. As shown in Fig.
3(b,c) and Fig. S10, the scaling relationships of (9) and (10) with $\epsilon$
and $G$ are verified by numerical analysis, thus confirming the power of our
theory to capture the key mechanical parameters that determine the onset and
the characteristics of localized solutions.

### Energetics of localized solution

We also studied the energetics of the single localized and uniform buckled
solutions (the energetics of higher-order solutions is discussed in SI Sec.
S2; also see SI Sec. S2 for details on energy calculation procedure). In Fig.
3(d), we plot the energy difference between localized and uniform buckling
deformations,
$(\mathcal{E}_{\rm{loc}}-\mathcal{E}_{\rm{uni}})/\mathcal{E}_{\rm{uni}}$, as a
function of excess strain $\epsilon-\epsilon^{\star}$ at fixed values of shear
modulus $G$. We find that there is a critical strain $\epsilon$ above which
the localized solution has lower energy than the uniform solution; this
critical strain increases with increasing $G$ or decreasing $K$. In fact, the
energy difference between the localized and the uniform states is set by an
interplay between shear energy and extensional energy of the substrate
springs; localizing the strain of the filament increases the amount of shear,
while decreasing the extensional energy of the substrate springs (see SI Sec.
S3.7).

Similarly we can also pick a solution for a given $\epsilon$ and study how
varying $G$ impacts its total energy. Since we do not perturb the solution it
will remain trapped in the local minimum. We find that below a critical $G$,
the energy difference between localized and uniform states becomes negative,
i.e. the localized state has a lower energy compared to the uniform solution
(Fig. 3(e)). The critical $G$ increases for increasing $\epsilon$, hence
expanding the window of values for shear modulus for which we can expect to
find localized states.

### Phase diagram for localized deformations

In summary, our analysis implies that localized deformations are stable
solutions for intermediate values of shear modulus $G$ (Fig. 4). By analyzing
the energy of uniform and localized solutions, we found that decreasing $G$
from the $G\to\infty$ limit below a critical value $G^{\rm{crit}}$ causes
localized deformations to have a lower total energy than the homogeneous
buckled solution. This critical value is given by (see SI Sec. S3.7)

$(G/K)^{\rm{crit}}\propto\epsilon-\epsilon^{\star};$ (11)

it increases with increasing $\epsilon$. Moreover, our perturbative solution,
(8), valid for low $G$, shows that, in this limit, decreasing $G$ causes the
width $w$ of localization to increase; eventually, localization disappears
when $w$ approaches system size $L$; using (9), we find that this occurs below

$(G/K)^{\star}\propto(\epsilon-\epsilon^{\star})^{-1}(\lambda/L)^{4}.$ (12)

These two effects give rise to the phase diagram in Fig. 4(a), where localized
deformations emerge for intermediate values of $G$. We have verified these
relationships using computer simulations both in the high $G$ (Fig. 4(b)) and
low $G$ regimes (Fig. 4(c)).

We note that our system displays bistability. If we prepare the homogeneous
state at $G=0.001$, which is well below the value at which it is the global
minimum (Fig. 4(b)), we find that random perturbations of the filament, up to
the buckling amplitude, relax back to the homogeneous state. We have similarly
verified the stability of the localized state at $G=0.005$ and found that this
state relaxes back for amplitudes of perturbations up to the buckling
amplitude. At higher values of $G$ we recover homogeneous buckling or less
pronounced localized buckling.

## Experiments

Using a multi-material 3D printer (Stratasys Connex500), we constructed a
simple experimental realization of a growing filament bound to a shearable
substrate to asses the occurrence of localized deformations (Fig. 5(a)). The
3D printed samples consisted of three parts: (i) a rigid, non-swellable
substrate, (ii) a flexible and swellable “filament”, and (iii) a series of
flexible and non-swellable “springs” connecting the filament to the substrate
Guiducci_2015 . During the fabrication process, a photosensitive liquid
precursor (the 3D printer ink) is deposited in a voxel-by-voxel fashion.
Several precursors are used to print multiple materials with different
properties and the resulting modulus can be tuned by varying the concentration
of photo-initiator. A UV light cross-links the liquid precursors in a layer-
by-layer fashion and this process is repeated until the full 3D model is
built. Depending on the liquid precursor composition and the degree of cross-
linking, a broad range of mechanical properties can be achieved from stiff
thermoplastic-like to soft rubber-like materials. The degree of cross-linking
also directly influences the swelling capacity of the polymer, and for the
experiments described here, the component formulations were guided by previous
formulations developed by Guiducci, et al. Guiducci_2015 . Growth of the
filament is induced through swelling by submerging the sample in isopropyl
alcohol, and the sample deformations were filmed over 1 hour with a digital
camera. The structures were elevated off the bottom of the container to ensure
that no flexible component was in direct contact with the walls, which would
cause friction effects that could limit the motion of the filament or the
springs. The swelling experiments were conducted on two different geometries
with the same filament length $L$, but consisting of $n_{s}=28$ and $n_{s}=56$
vertical springs (Figs. 5(b,c)), i.e. the spacing $b$ between springs is
changed by a factor of 2 between the two experimental realizations. In both
cases, $b$ was sufficiently large to avoid differences in solvent diffusion
between the two samples. Since the material properties as well as thickness of
all components were identical for both realizations, changing $b$ allows us to
control the stretching modulus of the filament, $S=\bar{S}b$, as well as the
stretching and shear moduli of the springs, $K=\bar{K}/b^{2}$ and
$G=\bar{G}/(h_{0}^{2}b)$, while keeping the other parameters unchanged.
According to our theory, the width, wavelength and amplitude of localized
deformations scale as $w\propto G^{-1/4}$, $\lambda\propto K^{-1/4}$ and
$A\propto(G/K)^{1/4}$; hence, we predict that reducing $b$ results in
localized deformations with overall smaller width, wavelength and amplitude.
These predictions agree with the experimental observations (Figs. 5(b,c)).
Moreover, the critical strain for the onset of deformations scales as
$\epsilon^{\star}=\sqrt{BK}/S$; hence decreasing the spacing $b$ between
springs is expected to increase $\epsilon^{\star}$. This prediction is
consistent with the experimental observation that sample (c) requires a longer
time for the onset of localized deformations compared to (b). For both
geometries, we measured the amplitude of deformations as a function of excess
growth strain $\epsilon-\epsilon^{\star}$ and verified the square-root scaling
of the bifurcation predicted by (10) (Fig. 5(d)).

## Discussion

We have described and studied a minimal realization of spatially localized
deformations in a growing filament connected to a shearable substrate using a
discrete rod model, a continuum theory, and implemented the results using
experiments using a multi-material 3d printing framework. Our main results
take the form of scaling relationships for key parameters relating to
localized deformations, such as amplitude, wavelength, and width, and a phase
diagram for uniform and localized deformations. These can be particularly
useful as guiding principles for designing controlled deformations with
potential applications to mechanical memories. Indeed, our work shows that the
ability to achieve localized deformations in filaments bound to shearable
substrates depends on specific combinations of the mechanical parameters and
the size of the system. Upon appropriate rescaling using the scaling
relationships derived in this study, localized deformations could be
appropriately designed both at the macro- and micro-scales.

Our work complements previous extensive mathematical literature on localized
deformations in a variety of non-equilibrium settings Pomeau_1986 ;
Coullet_2000 ; Coullet_2003 ; Swift1977 ; Cross1993 ; diBernardo2008 ;
Devaney1977 ; Lee1993 , by providing a simple physical realization of such
localized deformations in equilibrium mechanical systems in a one-dimensional
setting, complementing earlier work by us in two-dimensional settings
Chung_2018 . Our work also provides an example of how to employ multi-material
3D printing as a reproducible, rapid and easily realizable method for studying
complex problems related to mechanically constrained growth at multiple
scales. Due to the ability to control the deposition of material with very
high resolution to create complex 3D structures, this approach could thus
provide a practical research platform for investigating mechanical feedback on
growth under different conditions, e.g. in the presence of gradients in
material properties. By incorporating extensions to 2D systems, our approach
could move us one step closer to understanding the two way-feedback between
mechanics and growth/swelling kinetics, which is a defining feature in the
growth of spatially extended structures in both materials science and biology.

## Acknowledgments

We acknowledge financial support from the Swiss National Science foundation
(TCTM), the Netherlands Organization for Scientific Research (NWO-FOM) within
the program “Barriers in the Brain: the Molecular Physics of Learning and
Memory” (No.FOM-E1012M, RK and CS), and the NSF BioMatter DMR 33985 (LM) and
MRSEC DMR 14–20570 (LM). We thank Nicholas W. Bartlett for technical support
during the swelling experiments.

## References

  * (1) E. Winkler. Die Lehre von der Elastizität und Festigkeit. (Dominicus, Prague, 1867).
  * (2) H. Diamant and T. A. Witten. Compression induced folding of a sheet: An integrable system. Phys. Rev. Lett. 107(16):164302 (2011).
  * (3) P. J. Yoo and H. H. Lee. Complex pattern formation by adhesion-controlled anisotropic wrinkling. Langmuir 24(13):6897-6902 (2008).
  * (4) F. Brau, P. Damman, H. Diamant, and T. A. Witten. Wrinkle to fold transition: influence of the substrate response. Soft Matter 9(34):8177-8186 (2013).
  * (5) L. Pocivavsek, R. Dellsy, A. Kern, S. Johnson, B. Lin, K. Y. C. Lee, and E. Cerda. Stress and fold localization in thin elastic membranes. Science 320(5878):912-916 (2008).
  * (6) G. W. Hunt, M. A. Peletier, A. R. Champneys, P. D. Woods, M. A. Wadee, C. J. Budd, and G. Lord. Cellular buckling in long structures. Nonlin. Dynam. 21(1):3-29 (2000).
  * (7) G. W. Hunt and M. A. Wadee. Localization and mode interaction in sandwich structures. Proc. R. Soc. Lon. A: 454, 1197-1216 (1998).
  * (8) A. Takei, L. Jin, J. W. Hutchinson, and H. Fujita. Ridge localizations and networks in thin films compressed by the incremental release of a large equibiaxial pre-stretch in the substrate. Adv. Mater. 26(24):4061-4067 (2014).
  * (9) Y. Cao and J. W. Hutchinson. Wrinkling phenomena in neo-hookean film/substrate bilayers. J. Appl. Mech. 79(3):031019 (2012).
  * (10) L. Jin, S. Cai, and Z. Suo. Creases in soft tissues generated by growth. Europhys Lett 95(6):64002 (2011).
  * (11) L. Jin, A. Takei, and J.W. Hutchinson. Mechanics of wrinkle/ridge transitions in thin film/substrate systems. J. Mech. Phys. Solids 81:22-40 (2015).
  * (12) A. Auguste, L. Jin, Z. Suo, and R.C. Hayward. The role of substrate pre-stretch in post-wrinkling bifurcations. Soft Matter 10:6520-6529 (2014).
  * (13) E. Knobloch. Spatial localization in dissipative systems. Annu. Rev. Cond. Matt. Phys. 6(1):325-359 (2015).
  * (14) Y. Pomeau. Front motion, metastability and subcritical bifurcations in hydrodynamics. Physica D 23:3-11 (1986).
  * (15) P. Coullet, C. Riera, and C. Tresser. Stable static localized structures in one dimension Phys. Rev. Lett. 84:3069-3072 (2000).
  * (16) P. Coullet, C. Riera, and C. Tresser. A new approach to data storage using localized structures., Chaos 14:193-198 (2004).
  * (17) J.Y. Chung, A. Vaziri, and L. Mahadevan. Reprogrammable Braille on an elastic shell. Proc. Natl. Acad. Sci. USA 115:7509-14 (2018).
  * (18) P.G. Drazin, Introduction to Hydrodynamic Stability (Cambridge University Press, 2002).
  * (19) E. Cerda and L. Mahadevan. Geometry and physics of wrinkling. Phys. Rev. Lett. 90(7):074302 (2003).
  * (20) L.D. Abramowitz and E.M. Lifshitz, Theory of elasticity (Pergamon Press, Oxford, 1970).
  * (21) M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (Dover, New York, 1965).
  * (22) A. Bergeon, J. Burke, E. Knobloch, and I. Mercader. Eckhaus instability and homoclinic snaking. Phys. Rev. E 78(4):046201 (2008).
  * (23) J. Burke and E. Knobloch. Snakes and ladders: localized states in the Swift-Hohenberg equation. Phys. Lett. A 360(6):681-688 (2007).
  * (24) L. Guiducci, J.C. Weaver, Y.J.M. Bréchet, P. Fratzl, and J.W.C. Dunlop. The geometric design and fabrication of actuating cellular structures. Adv Mater Interf 2, 1500011 (2015).
  * (25) J. Swift, P.C. Hohenberg (1977) Hydrodynamic fluctuations at the convective instability. Phys. Rev. A. 15:319-328.
  * (26) M.C. Cross, P.C. Hohenberg (1993) Pattern formation outside of equilibrium. Rev. Mod. Phys. 65:851.
  * (27) di Bernardo M, Budd CJ, Champneys AR, Kowalczyk P (2008) Qualitative theory of non-smooth dynamical systems. Appl Math Sci 163:47-119.
  * (28) Devaney RL (1977) Blue sky catastrophes in reversible and Hamiltonian systems. Indiana Univ Math J 26:247-263.
  * (29) Lee KJ, McCormick WD, Ouyang Q, Swinney HL (1993) Pattern formation by interacting chemical fronts. Science 261:192-194.