Fluid-supported elastic sheet under compression: Multifold solutions
Leonardo Gordillo Edgar Knobloch





leonardo.gordillo@usach.cl


Departamento de Física, Universidad de Santiago de Chile, 

 Av. Ecuador 3493, Estación Central, Santiago, Chile


knobloch@berkeley.edu


Department of Physics, University of California at Berkeley, Berkeley,
California 94720, USA

The properties of a hinged floating elastic sheet of finite length
under compression are considered. Numerical continuation is used to
compute spatially localized buckled states with many spatially localized
folds. Both symmetric and antisymmetric states are computed and the
corresponding bifurcation diagrams determined. Weakly nonlinear analysis
is used to analyze the transition from periodic wrinkles to single
fold and multifold states and to compute their energy. States with
the same number of folds have energies that barely differ from each
other and the energy gap decreases exponentially as localization increases.
The stability properties of the different competing states are established
and provide a basis for the study of fold interactions. 


68.08.-p: Liquid-solid interfaces, 46.32.+x: Static buckling and
instability, 47.54.-r: Pattern selection; pattern formation



Fluid-supported elastic sheet under compression: Multifold solutions
    Edgar Knobloch
    December 30, 2023
====================================================================





§ INTRODUCTION


Recent experiments <cit.> have shown that an elastic
sheet on the surface of a liquid undergoes a series of transitions
when compressed in the longitudinal direction: for small compression
the sheet deforms into a periodic array of wrinkles; with increasing
compression the deformation of the sheet becomes more and more nonlinear
and the deformation spontaneously localizes. This process continues
until the sheet makes contact with itself.

These experiments have stimulated several efforts at modeling this
sequence of transitions <cit.>.
The equilibrium configurations are described by a fourth-order equation
for the deformation angle as a function of arclength. This equation,
as well as others describing post-buckling phenomena of elastic media
under compression <cit.> and torsion <cit.>,
is of Swift-Hohenberg type and its solutions bear a number of similarities
to existing analysis of steady states of the Swift-Hohenberg equation
with a quadratic nonlinearity <cit.>, or equivalently
to the small amplitude behavior of the quadratic-cubic Swift-Hohenberg
equation (SH23) as described, for instance, in Refs. <cit.>.
Specifically, on an infinite domain one expects that the primary instability
to the wrinkled state is accompanied by the simultaneous appearance
of two branches of spatially localized solutions, one characterized
by deformation that peaks in the center of the domain and the other
that dips in the center. The theory tells us that on an infinite domain
these states are described by the normal form for the reversible Hopf
bifurcation (in space) with 1:1 spatial resonance <cit.>.
In particular, the branches of localized solutions are initially exponentially
close to one another <cit.>, but gradually
separate as the compression of the sheet increases and the solutions
begin to localize.

Of course, in the experiments the sheet is of finite length. In this
case existing theory shows that the 1:1 reversible Hopf bifurcation
breaks apart into a primary bifurcation to a wrinkled state followed,
at small amplitude, by a secondary bifurcation that creates a pair
of spatially modulated wrinkled states, one with maxima in the center
and the other with minima in the center. This bifurcation is an example
of an Eckhaus instability of the periodic wrinkled state and its dependence
on the domain length has been studied in detail <cit.>.
In particular, it is known that subsequent secondary bifurcations
lead to the so-called multipulse states, i.e., states with more than
one localized structure within the available domain <cit.>.
However, since the periodically wrinkled state does not exhibit a
saddle-node bifurcation the correspondence with SH23 does not extend
to larger amplitude states; in particular, homoclinic snaking <cit.>
is absent.

In this paper we use the similarity between the present problem and
the Swift-Hohenberg equation to identify new classes of localized
wrinkled structures on finite domains <cit.>. We focus
on the case of a thin elastic sheet allowing us to neglect its contribution
to the gravitational potential energy when its is deformed. In this
case the model equation is odd in the dependent variable, a symmetry
respected by the hinged boundary conditions we employ. As a consequence
the system admits four types of localized states instead of the two
mentioned above. two even parity states under spatial reflection and
two odd parity states. Of these the even states correspond to deformation
that either peaks or dips in the center of the domain while the odd
parity states consist of symmetry-related peak-trough combinations.
Consequently the system exhibits four branches of spatially
localized states, and so its properties in fact resemble more closely
those of the small amplitude states in the cubic-quintic Swift-Hohenberg
equation <cit.> rather than those of the
quadratic-cubic Swift-Hohenberg equation.

Each of these solutions is a single fold solution, with a fold either
in the center of the sheet or split between the two end points. However,
the problem also admits a vast array of multifold states consisting
of both identical and nonidentical folds at different locations along
the sheet. We compute here both single fold and multifold states using
numerical continuation techniques that have proved invaluable in the
studies of the Swift-Hohenberg equation, as well as weakly nonlinear
theory of the type first employed in the present context in <cit.>
and discuss the stability properties of these distinct states.

The paper is organized as follows. In the following subsection we
summarize the basic properties of the model we use. We then focus
on finite sheets, and discuss in turn periodic structures in sheets
with hinged ends and the wrinkle-to-fold transition. The word 'fold'
is here employed to mean a strongly spatially localized deformation.
We compute fold solutions with both positive and negative deformation
in the center, as well as multifold structures. We also compute the
free energy of these structures, and identify structures corresponding
to minima of the free energy with the stable states observed in experiments.
Our results can be considered to be an extension of the work of <cit.>
to finite domains and extend the analysis in <cit.>
to multifold states.



 §.§ Governing equations 


Consider an elastic sheet of length L and bending modulus B
lying on a semi-infinite fluid substrate of density ρ. In static
equilibrium under gravity, a compression Δ in the x direction
will deform the sheet in the y direction as shown in Fig. <ref>.
The shape of the sheet can be parametrized in terms of the arclength
s and the local angle ϕ(s) between the sheet and the horizontal
axis. Cartesian coordinates [x(s),y(s)] are related
to this parametrization through x(s)=∫_-L/2^scosϕ(s) s,  y(s)=∫_-L/2^ssinϕ(s) s.
Accordingly, the total energy stored (per unit of length in the z
direction) due to the deformation is 

    E[ϕ]=∫_-L/2^+L/2s (1/2Bϕ̇^2+1/2ρ gy^2cosϕ),

where the dot denotes derivative with respect to s and g is
the acceleration due to gravity. The first term is the bending energy
while the second represents the gravitational potential energy of
the deformed state. The sheet shape is linked to the total compression
through 

    Δ[ϕ]=∫_-L/2^+L/2s (1-cosϕ).

Thus the load P required to keep the system in equilibrium is equal
to dE/dΔ.

The natural lengthscale of the system, λ=[B/(ρ g)]^1/4,
can be used to nondimensionalize the quantities x,y,s,Δ and
L. Likewise, the energy E can be scaled by B/λ and
P by B/λ^2. This rescaling leads to a simpler version
of Eq. (<ref>) with B=ρ g=1.



We study quasistatic deformations of an elastic sheet under the following
two conditions: (i) a fixed load P (dead loading), or (ii)
fixed compression Δ (rigid loading). An equation for
the shape of the sheet, ϕ(s), can be derived in both cases using
a Lagrangian approach as shown in <cit.>.
Whether minimizing the functional E[ϕ], the energy, for given
Δ (rigid loading) or, equivalently, minimizing the functional
G[ϕ]≡ E-PΔ, the free energy, for given
P (dead loading), one finds that ϕ(s) satisfies 

    ⃜ϕ+3/2ϕ̇^2ϕ̈+pϕ̈+sinϕ=0,

where p depends on the boundary conditions. We adopt here hinged
boundary conditions, i.e., the conditions .y|_± L/2=.ÿ|_± L/2=0
(or, equivalently, .ϕ̇|_± L/2=.⃛ϕ|_± L/2=0),
so p is related to the static load P by <cit.>

    P=pcosϕ(± L/2)-ϕ̈(± L/2)sinϕ(± L/2).

With these boundary conditions the system is invariant under the transformations
ϕ→-ϕ and s→-s, which are linked to
invariance with respect to x→-x and y→-y.

Equation (<ref>) conserves the following quantity along
the sheet: 

    H[ϕ]=⃛ϕ-1/2^2+3/8ϕ̇^4+1/2pϕ̇^2-cosϕ.

For hinged sheets H=.-|_± L/2-cos.ϕ|_± L/2.
Throughout this article both dead and rigid loading approaches will
be considered and the results will therefore be represented either
in terms of the quantities (P,G) or (Δ,E),
as appropriate. The difference between these two approaches is fundamental
importance for the stability analysis as discussed in Sec. <ref>.



 §.§ Exact solutions for an infinite sheet


Despite the presence of nonlinear terms, Eq. (<ref>) admits
explicit solutions for the case of an infinite sheet. Two families
of spatially localized solutions can be found <cit.>:
a family of symmetric solutions, 

    ϕ_s(s)=4arctan[κsin k(s-s_0)/kcoshκ(s-s_0)]

and a family of antisymmetric solutions: 

    ϕ_a(s)=4arctan[κk(s-s_0)/kcoshκ(s-s_0)].

Here s=s_0 corresponds to the symmetry point while k=1/2√(2+P)
and κ=1/2√(2-P), i.e., the solutions exist only
for P<2. Evaluation in the limit s→±∞ shows
that for this case P=p. The extra parameter s_0 is a result
of the invariance of the infinite sheet under transformations of the
form s→ s+ const. In the literature, localized solutions
are referred to as folds. These solutions are members of a
wider family of solutions, obtained by adding an extra phase to the
trigonometric function, i.e., 

    ϕ_f(s)=4arctan[κ[k(s-s_0)-φ_f]/kcoshκ(s-s_0)].

Thus the symmetric and antisymmetric solutions correspond to the particular
cases φ_f=π/2 and φ_f=0, respectively. The
asymmetric family (φ_f≠0,π/2) has been shown to describe
the shape of an elastic sheet when extracted from a fluid bath <cit.>.

On the other hand, periodic solutions on an infinite sheet can be
found in terms of Jacobi elliptic functions, 

    ϕ(s)=2arcsin[k_k(qs-ϑ_0)],

or, equivalently, 

    y(s)=2k/q_k(qs-ϑ_0),

and exist for 0<k<1. The quantity q is related to k and the
static load P by the implicit relation 

    P=q^2+(1-2k^2)q^-2.

Thus q^2-q^-2<P<q^2+q^-2. The solutions (<ref>)
have a period equal to 4(k)/q, where (k) is the
complete elliptic integral of the first kind. Periodic solutions are
usually referred to as wrinkles in the literature, cf. <cit.>
and references therein.



 §.§ Explicit periodic solutions for finite sheets


Exact periodic solutions of Eq. (<ref>) for finite sheets
can easily be found by matching properly both ϑ_0 and
the period of Eq. (<ref>) with the boundary conditions.
For hinged boundary conditions, it follows that 

    ϑ_0=ν_n(k),   q_n=2n(k)/L,   n=1,2,3...,

where ν_n=1 for n odd and ν_n=0 for n even (see
<cit.> for details). As expected, the finiteness of
the sheet discretizes the parameter q into discrete families q_n(k)
labeled by the index n (the number of half-periods of the solution).
Solutions of this type will be referred to as ϕ_n. Here k
quantifies the amplitude of the solution as measured, for example,
by the mean compression per unit length Δ≡Δ/L:

    Δ=2I(k).

Here I(k)≡1-(k)/(k) and (k) is the complete elliptic
integral of the second kind. Remarkably, this expression does not
depend on n.

The mean energy per unit length, E≡ E/L, is given
by 

    E=2(q_n^2+q_n^-2)[k^2-I(k)]-[a_k-b_kI(k)],

where a_k=8q^-2(k^4+k^2)/15 and b_k=16q^-2(k^4-k^2+1)/15.
The mean free energy G=E-PΔ
now follows from Eqs. (<ref>)–(<ref>).
For a given n and L, Eqs. (<ref>)–(<ref>)
determine (P,Δ,E,G)
in terms of the single continuous parameter k. Figure <ref>
depicts the resulting solution branches in the (Δ,E)
plane for L=20π and different n values.



Periodic solutions for finite sheets bifurcate from the trivial state
ϕ(s)=0 at k=0. Since E(0)=K(0)=π/2, each periodic branch
ϕ_n bifurcates from the P-axis at 

    P_n^*=q_n^2+q_n^-2

and Δ_n^*=0. Thus P_n^*≥2 with equality at n=L/π.
To determine which solution emerges spontaneously when the sheet is
compressed, we expand Eqs. (<ref>) and (<ref>)
around k=0, yielding E≈(q^2+q^-2)k^2
and Δ(k)≈ k^2. It follows that
E≈(q^2+q^-2)Δ.
The branch with minimal energy in the (Δ,E)
plane near the trivial state is that which minimizes the slope q_n^2+q_n^-2≈ P_n^*
of the energy. Hence, the branch bifurcating at the lowest P_n^*
value is the only one that is stable in the vicinity of the primary
bifurcation.

The branches ϕ_n and ϕ_n+1 bifurcate simultaneously
from the trivial state whenever L_(n,n+1)=π√(n(n+1)) and
do so at P_(n,n+1)^*=2+[n(n+1)]^-1. These conditions
define the range of L for which the nth solution is stable:
L_(n-1,n)<L<L_(n,n+1). As L→∞, the width
of this range converges to π and its center to π n.

Far from the primary bifurcation, corresponding to k=0, other branches
may display lower energies than the branch that minimizes P_n^*.
This is illustrated in Fig. <ref>, which shows that
the lower n branches eventually cross the n=20 branch as Δ
increases.



§ FOLDS IN FINITE SHEETS




  §.§.§ Numerical continuation of localized patterns


Since no explicit exact localized solutions of Eq. (<ref>)
exist for elastic sheets of finite length and hinged boundary conditions,
we employ numerical continuation <cit.>
to follow different solution types of this problem through parameter
space. For this purpose, we implemented Eq. (<ref>) in the
AUTO^© code <cit.> with hinged boundary
conditions, with the independent variable s transformed so the
original half domain [-L/2,0] becomes [0,1].
Since the quantities (Δ,E,G,P)
are all functions of the parameters (p,L) of the system,
we plot the results as a function of P instead of p, for fixed
length L=lπ (l∈ℝ).



Numerical continuation requires an initial guess. Taking advantage
of the existence of explicit periodic solutions, we performed continuation
on a 128 point spatial grid starting from solutions (<ref>)
for a given n and a high value of k, i.e., ϕ_n far from
the trivial state. The code was tuned to detect bifurcation points
along this branch as P increases towards the onset value P^*.
Figure <ref>(a-b) displays the branches
ϕ_40 and ϕ_41 in black, both for L=40π. Each
periodic branch displays a set of bifurcation points with decreasing
separation as P increases towards P^*. For ϕ_40 these
bifurcation points are shown Fig. <ref>(b). These
bifurcation points will be referred to as M_n,m, where the integer
n specifies the primary periodic state and the integer m counts
the bifurcations from this state starting from the primary bifurcation
P_n^*. Thus M_40,1 denotes the secondary bifurcation closest
to P_40^*.

At each of the secondary bifurcation points M_n,m we switched
the branch direction and ran the numerical continuation for decreasing
P until P=0 was attained. The new branches are also displayed
in Fig. <ref> but now in color. Some of the solutions far
from the bifurcation points M_n,m are depicted in Fig. <ref>.
A simple inspection [see inset in Fig. <ref>(a)]
shows that the branch starting at M_40,1 corresponds to a localized
antisymmetric fold while that starting at M_41,1 is a symmetric
fold. The former bifurcation point is located quite far from P_40^*,
i.e., at finite amplitude, while the latter lies very close to P_41^*.
At this value the ϕ_41 state barely departs from the trivial
state. For L=41π the situation is the opposite: the symmetric
fold is the one that bifurcates at finite amplitude (from the ϕ_41
branch), while the antisymmetric fold bifurcates at a very small amplitude
(from the ϕ_42 branch). In any case, the symmetric and antisymmetric
branches approach each other very fast as they depart from their respective
branch point with decreasing P. 


Solutions ϕ_n,m with m>1 represent m-fold solutions,
i.e., arrays of m equispaced folds. Examples are shown in Fig. <ref>.
Although reported in <cit.>, these multifold solutions
have not been described in detail. Multifold solutions ϕ_n,m
for a fixed l do not exist for arbitrary n and m. For instance,
for l=40, one-fold solutions (m=1) are observed only for n=40
and 41 while two-fold solutions (m=2) are only found for n=39,40,41,42.
In general, 2m+2 solutions with m folds were detected for n
varying from l-m+1 to l+m.



 §.§ Weakly modulated solutions


Multiscale analysis describes successfully several features of this
system and, in particular, the wrinkle-to-fold transition <cit.>.
Herein, we use this approach to understand the emergence of the solutions
obtained in Sec. <ref>. The system is
assumed to be subject to hinged boundary conditions and finite. The
domain of s is mapped into s∈[0,L] for the sake
of simplicity in the application of the boundary conditions.

We seek multiscale solutions of Eq. (<ref>) in the form

    y(s)=∑_j=1^∞ϵ^ju_j(s,S=ϵ s),

where the functions u_j are of 𝒪(1) and
ϵ is a small parameter measuring the distance from some
primary instability threshold, ϵ^2≡ p^*-p. 
Since ϕ(s)=arcsin(∂_sy) with ∂_s
now replaced by ∂_s+ϵ∂_S, Eq. (<ref>)
generates a hierarchy of equations for the functions u_j(s,S).
The first three equations are 

    𝒪(ϵ^1):  ℒ_p^*∂_su_1    =     0,
    𝒪(ϵ^2):  ℒ_p^*∂_su_2    =     -ℒ'_p^*∂_su_1,
    𝒪(ϵ^3):  ℒ_p^*∂_su_3    =     -ℒ'_p^*∂_su_2-𝒩_p^*∂_su_1,

where the fourth-order linear operators ℒ_p^* and
ℒ'_p^* and the nonlinear operator 𝒩_p^*
are, respectively, defined as 

    ℒ_p^*w    ≡   (∂_s^4+p^*∂_s^2+1)w,
    ℒ'_p^*w    ≡   (5∂_s^4+3p^*∂_s^2+1)w,
    𝒩_p^*w    ≡   1/2⃜ww^2+10ẅ”+15/2ẅw”+3ẅ^2w+…
           +4⃛wẇw-ẅ+p^*[1/2ẅw^2+3w”+ẇw].

Here the overdots denote derivatives with respect to s while the
primes denote derivatives with respect to S. At the instability
threshold p^*=2 the linear operators reduce to ℒ_2≡(∂_s^2+1)^2
and ℒ_2'≡(∂_s^2+1)(5∂_s^2+1)
and have a common kernel spanned by the null eigenfunctions of the
second order linear operator ∂_s^2+1. In this case the
𝒪(ϵ^1) solution of (<ref>)
is u_1(s,S)=A(S)sin(s-θ),
i.e., the dominant term u_1(s,S) is a carrier wave with wave
number q=1 modulated by a slowly varying envelope A(ϵ s).
It now follows that the right side of the 𝒪(ϵ^2)
equation vanishes, and hence that u_2≡0. An equation for
A(S) is therefore obtained by imposing a solvability condition
at 𝒪(ϵ^3). A straightforward calculation
yields 

    A”-1/4A+1/8A^3=0.

Equation (<ref>) has two solutions given by Jacobi elliptic
functions, 

    ϵ A_(S)     =     4ϵκ__k(κ_S-ϑ_),
    ϵ A_(S)     =     4ϵκ_^1/2k_k(κ_S-ϑ_),

where ϑ is an arbitrary phase and 

    ϵ^2κ_^2    =    1/4(2-p)(2-k^2)^-1
    ϵ^2κ_^2    =    1/4(2-p)(2k^2-1)^-1.

Thus each solution branch is parametrized by k, just as for the
periodic case in Sec. <ref>. The final step is to multiply
the expressions for the modulation amplitude A by the carrier wave,
sin(qs-θ), and match the result to the hinged boundary conditions.
Note that both the modulation amplitude and the carrier wave are periodic
functions 0f s.



  §.§.§ Wave number and phase selection


Since the wave number q=1, the weakly modulated solutions satisfy
hinged boundary conditions only for a simple set of combinations of
θ,ϑ,L and ϵκ, displayed in Table <ref>.
This set can be identified on realizing that the chosen boundary conditions
are satisfied if and only if one of the two periodic functions, either
the trigonometric or the elliptic one, has a node (zero) at the boundary
while the other has an antinode (local maximum or minimum) at the
same location. 






Since both elliptic functions localize as k→1, some of
the solutions from Table <ref> may display asymmetric folds
at the boundaries. This happens whenever the trigonometric function
has a node at the boundary. These types of solutions were also detected
in our numerical continuation. However, in the remainder of this article,
and for the sake of simplicity, we restrict our analysis to solutions
that do not display folds at the boundaries. The solutions that remain
are summarized in Table <ref>. Solutions related by the
symmetries x→-x, y→-y have been omitted. 








  §.§.§ Branches in parameter space


The weakly modulated solutions y(s)=A(ϵ s)sin(s-θ)
can be used to build expansions for the static load P, the mean
energy E, the mean compression Δ
and the mean free energy G for a sheet of finite length
L. For a static load P the condition p=2-ϵ^2 translates
into P=2-ϵ^2+𝒪(ϵ^4) using
Eq. (<ref>) and the solutions obtained in Sec. <ref>.
Accordingly, (E,Δ,G)=∑_j=1ϵ^j(E_j,Δ_j,G_j).
The expressions for (E_j,Δ_j,G_j)
can be simplified by averaging over one period of the fast oscillation,
e.g., ∫_0^Lds sin^2(s)A(ϵ s)=1/2∫_0^Lds A(ϵ s).
As a consequence, the first and third order terms all vanish and the
remaining quantities only depend on the envelope A(S). The following
expressions for the quantities P, G, Δ,
E to 𝒪(ϵ^4) are obtained:

    [                                P                                =                    2-ϵ^2+𝒪(ϵ^4),;                                G                                =     1/2ϵ^4[1/2A^2+2A'^2-1/8A^4],;                                Δ                                =  1/4ϵ^2A^2+1/4ϵ^4[3/16A^4+A'^2],;                                E                                = 1/2ϵ^2A^2+1/2ϵ^4[3A'^2+1/16A^4]. ]

The overbars indicate spatial averages. After replacing A by A_
or A_ from Eqs. (<ref>)–(<ref>), all the
averaged quantities can be expressed in terms of the complete elliptic
integrals K(k) and I(k) (see Appendix <ref>).
For our reduced set of solutions, ϵ=2mK(k)/(κ L), where
κ is either κ_ or κ_). Hence, all
solution branches are fully parametrized in terms of k. Inversion
of Eqs. (<ref>)–(<ref>) yields an expression
for the static load. For q=1 we obtain, through 𝒪(ϵ^4),

    [[ P_; P_ ]]     =     2-2^4m^2/L^2K^2(k)[[  2-k^2; 2k^2-1 ]],
    [[ G_; G_ ]]     =    2^8m^4/3L^4K^4(k)[[     3-k^2-2(2-k^2)I(k); 3k^4-k^2-2(2k^2-1)I(k) ]],
    [[ Δ_; Δ_ ]]     =    2^4m^2/L^2K^2(k)[[   1-I(k); k^2-I(k) ]]+
          2^6m^4/3L^4K^4(k)[[     9-2k^2-7(2-k^2)I(k); 9k^4-2k^2-7(2k^2-1)I(k) ]],
    [[ E_; E_ ]]     =    2^5m^2/L^2K^2(k)[[   1-I(k); k^2-I(k) ]]+
          2^7m^4/3L^4K^4(k)[[     3+2k^2-5(2-k^2)I(k); 3k^4+2k^2-5(2k^2-1)I(k) ]].

We use the above expressions to characterize the dnoidal and cnoidal
solutions in a bifurcation diagram, and the parameter k as a measure
of their amplitude. In particular, the condition k→0
determines the threshold for the appearance of each solution type.



 §.§ Wrinkle-to-fold transition


Experiments have shown that an elastic sheet over a fluid substrate
undergoes a wrinkle-to-fold transition when continuously compressed
from its original length <cit.>. The wrinkle-to-fold
transition is characterized by the emergence of a localized fold from
a periodic pattern. Figure <ref> shows the
critical loads P_n^* for the appearance of a periodic state
with n=38,39,...,42 wavelengths in the domain (solid lines), as
well as the secondary bifurcation thresholds M_n=l,m=1^
for the appearance of the lowest energy dnoidal states (m=1) obtained
from the k→0 limit of the dnoidal solutions of the previous
section (red filled circles). In contrast, the cnoidal solutions bifurcate
from the periodic states within an 𝒪(L^-3)
neighborhood of L=L_(n-1,n+1), the crossing points for like-parity
modes n-1 and n+1. The predicted bifurcation thresholds M_n=l,m=1^
(blue open circles) are also obtained from the k→0 limit
of the corresponding analytical solution. These bifurcation points
and the analytical expressions for the small amplitude solutions present
nearby were used as an input into a numerical continuation routine
and used to extend these thresholds to noninteger values of L/π
(M_n,m=1 in Fig. <ref>). Numerical continuation
reveals that the M_n,m=1 curves resemble parabolas with vertices
located in the vicinity of M_n=l,m=1^. These parabolas extend
between the codimension-two points L_(n-2,n) and L_(n,n+2)
(or approximately, between M_n=l-1,m=1^ and M_n=l+1,m=1^).
Also indicated are the codimension-two points L=L_(n,n+1) (open
red diamonds) corresponding to the crossing of periodic modes n
and n+1, i.e., the points at which the primary bifurcation to a
symmetric wrinkle state changes to an antisymmetric wrinkle state
as L changes. The curves M_n,m=1 and P_n^*, as well
as the predicted locations M_n=l,m=1^ and M_n=l,m=1^,
are all compressed towards the P=2 axis as L increases.



As shown in Fig. <ref> the multiscale analysis
generates a discrete set of bifurcation points instead of the continuous
curves of bifurcation points obtained by numerical continuation. The
solid curves for L∈(lπ,(l+1)π), l∈ℤ^+, can
only be obtained via numerical continuation. This limitation of the
multiscale analysis has its origin in the structure of the equations
analyzed in Sec. <ref>, which requires l∈ℤ^+
to satisfy simultaneously the boundary conditions and the cancellation
of the right-hand side of the 𝒪(ϵ^2)
equation in Eq. (<ref>). If we instead follow <cit.>
and minimize the functional G on the assumption that
y(s)=A(S)cos(qs) we obtain Eq. (<ref>) when q=1
and a similar equation to Eq. (<ref>) but with q-dependent
coefficients whenever q>1. The latter, however, yields inaccurate
solutions, which are an artifact of improper assumption on the s-dependence
of y(s) when q>1.



 §.§ The π/2-shifted localized branch


The exact nonlinear solutions in Eqs. (<ref>)–(<ref>)
show that the symmetric and antisymmetric folds on an infinite sheet
are equivalent from the energy point of view. Such states are therefore
energetically degenerate. However, this is no longer the case on a
finite sheet and the question arises whether whether a sheet whose
wrinkle-to-fold transition yields a symmetric (antisymmetric) fold
may also display an antisymmetric (symmetric) fold. Results from Sec. <ref>
provide hints about the dynamical features of these two states. For
l∈ℤ^+ symmetric and antisymmetric single fold states
(i.e., n=l,m=1) can be understood in terms of the carrier wave
phase θ. The idea is simple: since both Jacobi elliptic functions
of the approximate solutions in Table <ref> have a maximum
at the center of the sheet, the choice of whether θ=0 or θ=π/2
determines the symmetry of the central fold. For n=l even (odd),
the dnoidal (cnoidal) branch yields the antisymmetric (symmetric)
single fold state and the cnoidal (dnoidal) branch the symmetric (antisymmetric)
state. However, there is an important difference between these two
branches. The cnoidal branch, whose carrier wave is shifted by π/2
relative to the dnoidal branch, emerges from the trivial state through
a primary bifurcation to the right of the primary bifurcation to the
periodic state (and half as close to it as the dnoidal branch bifurcation
on the left). The resulting predictions of the theory of Sec. <ref>
are displayed in the inset of Fig. <ref>(a) and
compared there with the corresponding results from numerical continuation.
The accuracy of the approach is very good.

The theoretical predictions show that in the limit L→∞
both branches merge into a single one in terms of (Δ,E,G,P).
As L increases, the secondary bifurcation to dnoidal states is
pushed downward along the primary periodic branch while the cnoidal
bifurcation point is also pushed towards the primary bifurcation point.
These two bifurcation points coalesce in the limit L→∞.
The corresponding large amplitude dnoidal and cnoidal branches become
rapidly very close as P decreases, a prediction that is confirmed
by numerical continuation to values of ϵ far beyond ϵ=1,
i.e., P≪1.



 §.§ Multifold solutions




The modulation analysis carried out above also provides a clearer
picture of the emergence of multifold (m≥2) solutions for l≡(L/π)∈ℤ^+.
The condition q=1 requires that n=l, where the labels n and
m, introduced in Sec. <ref>, represent
the carrier and envelope wave numbers, respectively. The bifurcations
for m≥2 follow the same principle as the branches of localized
states when m=1, i.e., the dnoidal m-branch bifurcates from
a periodic state ϕ_n at finite amplitude while the corresponding
cnoidal m-branch bifurcates from the periodic state ϕ_n+m
at a very small amplitude, and to the right of the primary bifurcation
to ϕ_n and half as close to it as the bifurcation to the dnoidal
branch. For example, when m=2 the cnoidal branch in general bifurcates
from ϕ_42 and does so very close to the point P_42^*,
although it may bifurcate directly from the trivial state but only
at the crossing point of P_38^* and P_42^*. Likewise,
the bifurcation points to the m≥2 dnoidal branches spread out
from the primary bifurcation point to the nth periodic branch as
m^2. We also used numerical continuation to follow the M_n,m
bifurcation points as a function of L as done in Sec. <ref>
for M_n,1. Just as the single-fold cnoidal bifurcation point
lies in the vicinity of L_(n-1,n+1), the crossing
point of the primary bifurcations to the n-1 and n+1 periodic
branches, the m-fold cnoidal bifurcation point is located in the
vicinity of L_(n-m,n+m), the crossing point of the
primary bifurcations to the n-m and n+m periodic branches. Likewise,
the M_n,m bifurcation curves for noninteger L/π resemble
parabolas with vertices near M_l=n,m^ extending from L_(n-2m,n)
to L_(n,n+2m). The range of existence of the parabola clarifies
why the number of m-fold states observed is 2m+2 and which n
values are allowed for a given m and L.

The analysis also reveals another interesting feature. The choice
of l and m determines (a) if the solution consists of an array
of symmetric or antisymmetric folds, and (b) the local shape
of each fold. The rule for (a) is very simple: these states are allowed
if l is divisible by m. For the dnoidal solutions, if l/m
is even, the solutions will be an array of antisymmetric folds whereas
for l/m odd the solution will be an array of symmetric folds with
alternating orientation. If folds are allowed at both boundaries,
another possible solution, consisting of alternating antisymmetric
folds for l/m even, can be observed. For l=40, periodic multifold
states are observed for m=2,4,5,8,10,20, as displayed in Fig. <ref>
(m=20 is not shown). The figure also shows a state with folds at
either boundary, for m=8, to illustrate the complementary family
of periodic solutions.

When l is not divisible by m the folds forming the array become
distorted and resemble the family of arbitrary phase folds described
in <cit.> for an infinite length sheet. The rule
(b) for the local shape of each fold can be understood in
terms of the phase φ_f in Eq. (<ref>), which
is now fixed by the carrier wave. The jth fold of the dnoidal solutions
displays a local phase φ_f=-lπ(j-1/2)/m, whereas
for the cnoidal solutions φ_f=-π/2-lπ(j-1/2)/m.
Accordingly, the fundamental period of an array of folds in a dnoidal
solution is given by 2L/(l,2m) while the fundamental
period for cnoidal solutions is L/(l,m) provided
l is divisible by m and l/(l,m) is odd,
or 2L/(l,m) in all other cases. Examples of
distorted folds are displayed in Fig. <ref>
for L=40π. Solutions which do not satisfy q=1, e.g., ϕ_n=39,m
or ϕ_n=41,m in a L=40π domain, obey a similar rule for
the local phase as that obtained by swapping the integer l for
n, i.e., replacing l in φ_f by n to obtain ϕ_n≠ l,m.
Further analysis for this general case is required to understand this
feature.



 §.§ Discussion




  §.§.§ Energy gap


The energy of the dnoidal and cnoidal branches for l∈ℤ^+
and n=l far from the bifurcation can be found by expanding the
expressions for P, G, Δ and E
obtained in Sec. <ref> around k=1.
If we write k=1-ξ, ξ≪1, we obtain 

    g_m≡1/mG_m^±(P)     =    8/3(2-P)^3/2(1∓3/2ξ),
    
    e_m≡1/mE_m^±(Δ)     =     2δ_m-1/48δ_m^3(1∓3ξ),

where the - (+) corresponds to the dnoidal (cnoidal) branch
and δ_m≡Δ/m is a single-fold average compression.
Notice that the results are no longer expressed in terms of mean values.
At leading order both branches have exactly the same energy. Equation (<ref>)
shows that, for fixed P (dead loading), the G spectrum is a
set of equispaced free energies, where the quantity g_m≈8/3(2-P)^3/2
can be identified with the free energy of a single fold. On the other
hand, for fixed compression Δ (rigid loading), the energy
of a single fold is e_m≈2δ-δ^3/48. The energies
display maxima at δ_max=2^5/2, with e_max=2^9/2/3.
Numerical continuation shows that this point is related to self-contact
of the sheet. A comparison between the numerical results and the analytical
approach at leading order is shown in Fig. <ref> for
q=1. The overlap between the numerical continuation results and
the analytical approximation is remarkable. 




Figure <ref> shows the difference between the m=1
dnoidal and cnoidal branches obtained from numerical continuation
and compares the result with a theoretical prediction obtained by
inverting the series expansions for either P or Δ to obtain
ξ in terms P or Δ. The final result for the energy
gap between the dnoidal and cnoidal branches is 

    δ g_m   ≡    g_m^+-g_m^-=-64(2-P)^3/2e^-1/2l_m(2-P)^1/2,
    δ e_m   ≡    e_m^+-e_m^-=-δ_m^3e^-1/8l_mδ_m,

where l_m≡ L/m can be identified as the spacing between
the folds. The predicted exponential decay in both (P,G)
and (Δ,E) planes matches the numerical results
almost perfectly. The exponentially small residual difference is a
consequence of the finite sheet size and can be understood as the
interaction between the exponential tails of localized solutions <cit.>.

The general case q≠1 far from the bifurcation is more complicated
and is not well described by the weak modulation approach. Numerical
continuation shows that, at leading order, the energies for q≠1
are in fact exactly the same as those for the q=1 branches: the
solutions appear to shift their q towards q=1 as k→1,
a feature that is beyond the scope of the approach in Sec. <ref>.



  §.§.§ Stability analysis


The analysis carried out in the previous sections provides a framework
for understanding how single fold and multifold solutions bifurcate
from periodic states. However, this type of analysis does not address
the stability of the solutions, although we expect solutions with
lowest energy such as the single fold solution ϕ_40,1 for
L=40π to be stable. For example, it may happen that several solutions
are stable simultaneously. Thus a stability analysis provides a link
between the theoretical predictions and observations in experiments.





To explore the stability properties of the solutions identified above,
we carried out a numerical study based on a variational approach.
Since the equations of a fluid-supported elastic sheet are static
(time-independent), we consider a solution as stable if and only if
it corresponds to a local minimum of energy in the space of functions
satisfying all the constraints, as well as the boundary and symmetry
conditions imposed by the problem formulation. The stability of a
solution is thus guaranteed if the energy functional is positive definite
in its vicinity, i.e., if and only if E[ϕ_0+δϕ]>E[ϕ_0], ∀ δϕ ∈ 𝔽_ϕ,
where E is defined as in Eq. (<ref>), ϕ_0
is the solution under study, δϕ≠0 is the variation of
ϕ_0 and 𝔽_ϕ is the space of (continuous)
functions such that ϕ_0+δϕ satisfies the prescribed
constraints (e.g. Δ[δϕ]=0 under rigid loading),
and the boundary and symmetry conditions (hinged, symmetric or antisymmetric).

We define the inner product (f,g)≡∫_-L/2^L/2f· g ds=∫_-L/2^L/2f^†g ds.
It can be shown (see Appendix <ref>) that the functionals
E,G evaluated at y=y_0+δ y can be written as 

    E[y]=     E_0+(δ E/δ y,δ y)+(δ y,δ^2E/δ yδ y)+𝒪(δ y^3),
    
    G[y]=     G_0+(δ G/δ y,δ y)+(δ y,δ^2G/δ yδ y)+𝒪(δ y^3),

where E_0=E[ϕ_0], G_0=G[ϕ_0]
and all the functional derivatives are evaluated at ϕ_0 (equivalently,
y=y_0). The variation δ y is related to the corresponding
variation δϕ by the geometrical relation ∂_s(y+δ y)=sin(ϕ+δϕ)
(see Appendix <ref> for further details). Similar
relations can be written for P and Δ. By construction,
the terms in E or G linear in δ y vanish when y=y_0
is a solution of the original variational problem. Hence, under dead
loading (P fixed) the solution is stable provided 

    (δ y,δ^2G/δ y^2δ y)=(δ y,[δ^2E/δ y^2-P_0δ^2Δ/δ y^2]δ y)>0.

Here P_0=P[y_0]. In the rigid loading case, P is no longer
fixed and the stability condition becomes 

    (δ y,[δ^2E/δ y^2-P_0δ^2Δ/δ y^2]δ y)-δ P(δΔ/δ y,δ y)>0

with an extra constraint imposed by the fixed compression: 

    (δΔ/δ y,δ y)=0.

Since the operator δ^2G/δ y^2 in the square brackets
is a real symmetric operator, the inner products behave as a quadratic
form in function space. Thus, the spectral decomposition of the operator
yields a set of eigenvalues and the corresponding eigenfunctions containing
information about the stability of the solution.

In practice, we first recovered ϕ_0 and the corresponding
y_0 from the AUTO^© calculations and resampled them
into an N vector. The derivative operators from Eqs. (<ref>)-(<ref>)
were likewise discretized into N× N arrays using Fourier derivative
matrices. Using these results, we built a discrete N× N matrix
of the free energy functional δ^2G/δ y^2, 𝐆_N× N.
The stability of solutions in the dead loading case simply depends
on the spectral decomposition of 𝐆_N× N, that is,
on the eigenvalues of the standard eigenvalue problem 

    𝐆[δ y]=λ[δ y].


The rigid loading case requires further work because of the extra
constraint. In this case we discretized the compression condition
by introducing a matrix Δ_N×1' as the discrete version
of δΔ/δ y, and the (N+1)×(N+1)
bordered Hessian matrix 𝐇 given by 

    𝐇≡[[        0 -Δ_1× N';  -Δ_N×1'   𝐆_N× N ]].

As shown in <cit.>, the stability problem
can be studied via the generalized eigenvalue problem 

    𝐇[[ δ P; δ y ]]=λ[[      0 0_1× N;  0_N×1 I_N× N ]][[ δ P; δ y ]].

Thus the stability of a solution under dead or rigid loading can be
tested by checking the signs of the eigenvalues of the standard (Eq. <ref>)
or generalized eigenvalue (Eq. <ref>) problems, respectively.
If all the eigenvalues are positive, linear stability is guaranteed.



Within this framework, we now proceed to describe the stability properties
of the solutions for a L=40π sheet. Under dead loading, all the
solutions described so far are unstable to amplitude modes, which
is as expected as ∂ G/∂ P<0 for all the branches
studied. The rigid loading case is much more interesting, however.
In the following section, we describe the stability properties of
the periodic solutions ϕ_40 and ϕ_41, the symmetric
and asymmetric single fold solutions ϕ_40,1 and ϕ_41,1,
and the two-fold solutions ϕ_39,2, ϕ_40,2, ϕ_41,2
and ϕ_42,2. For the sake of completeness, we also examined
the stability of solutions displaying half-folds at the boundaries,
which we denote by the superscript b. For instance, the solution
ϕ_40,1^(b) has half-folds at s=± L/2, i.e.,
1/2+1/2=1 complete fold. 






  §.§.§ Stability of periodic solutions


Figure <ref> shows the six lowest eigenvalues
λ as a function of the load P for the periodic solutions
ϕ_40 (solid lines) and ϕ_41 (dashed lines). The corresponding
eigenfunctions evaluated at P=1.9 are shown in the right insets
following the same color scheme as used for the eigenvalues. The figure
indicates that the periodic solution ϕ_40 becomes unstable
at P=1.995, when its lowest eigenvalue becomes negative. The corresponding
eigenfunction represents a buckling mode that leads to either a centered
single fold state or a pair of half-folds on the boundaries. We will
refer to these eigenfunctions as fold eigenmodes. At P≈1.980,
a second eigenvalue becomes negative, which buckles the solution into
a pair of folds (or a 1/2+1+1/2-fold state);
at P≈1.956, a third instability sets in, buckling the solution
into three folds or a 1/2+1+1+1/2-fold state,
etc. The points where these solutions cross λ=0 are consistent
with the branch points identified using AUTO^© as leading
to solutions that subsequently evolve into fold solutions. Similar
behavior is observed for the periodic solution ϕ_41, with
fold eigenvalues crossing λ=0 sequentially as the number
of folds increases. Besides the fold eigenvalues all of which decrease
monotonically as P decreases, a second family of eigenvalues is
also present. However, these never cross λ=0 even though
they decrease monotonically as P decreases. The associated eigenfunctions
do not buckle the solution into localized folds but instead introduce
spatial modulation on top of the periodic pattern. We conclude that
the periodic solutions are stable with respect to spatial modulation,
at least when the wavelength of the periodic state is close to the
natural wavelength (no Eckhaus instability).

Figure <ref>(a) shows the construction of
a n=40 pseudo-solution arising from the second unstable fold mode.
The pseudo-solution is obtained by adding a multiple of the eigenfunction
(here the multiple is ±1) to the periodic state ϕ_40.
Panel (a) shows that when the multiple is +1 the instability
suppresses the oscillations in the center and near the boundaries,
leaving a 1+1 state. When the multiple is -1 the result is a
state resembling a 1/2+1+1/2 state.



  §.§.§ Stability of single fold solutions


The corresponding diagram for single fold solutions is shown in Fig. <ref>.
The four single fold solutions have positive eigenvalues and hence
are all stable. The first two solutions ϕ_40,1 and ϕ_40,1^(b),
have identical eigenvalues, a consequence of the fact that ϕ_40,1^(b)
can be obtained through a half-domain translation of ϕ_40,1.
This translation respects the boundary conditions at s=± L/2.
The lowest eigenvalue of ϕ_40,1 is inherited from the periodic
branch ϕ_40 and represents the amplitude mode of ϕ_40
[Fig. <ref>(b)]. 
The second lowest eigenvalue is of particular interest, because its
eigenmode is related to the splitting of the central fold into two
[Fig. <ref>(c)]. These eigenmodes,
which will be referred to as splitting eigenmodes, represent modulation
of a carrier wave that shifts continuously with respect to the base
state, with a π-phase shift at the center of the fold. With increasing
amplitude this mode leads to growing separation of the central fold
into two adjacent folds, in contrast to the usual fold eigenmodes
corresponding to a 0 or π shift. The third eigenvalue is linked
to a fold eigenmode which buckles the solution into a 1/2+1+1/2
array. Larger eigenvalues show similar behavior involving larger numbers
of folds. However, none of these modes is unstable, and so no finite
amplitude solutions of this type are present.

The stability behavior the ϕ_41,1 state is similar to that
of ϕ_40,1 although there are some important differences. First,
the solution ϕ_41,1^(b) no longer corresponds
to a simple half-domain shift of ϕ_41,1: to satisfy symmetry
conditions, a π/2 phase shift of the carrier wave plus
an inversion of the solution in one half of the domain is also required.
As a result the eigenvalues of ϕ_41,1 and ϕ_41,1^(b)
are no longer identical. The slight difference between them (≲10^-4)
is due to the presence of a defect in the center of the domain in
the latter case. The lowest eigenvalue mode is again related to the
periodic solution, ϕ_41,0. The next lowest eigenvalue mode
is now a 1+1 eigenmode while the third is a splitting eigenmode.
Notice also that the second mode has the lowest eigenvalue for P⪅1.982
but this eigenvalue, like the others, remains positive and no instability
is triggered.





  §.§.§ Stability of two-fold solutions


We also studied the stability properties of the two-fold solutions.
The eigenvalues of the ϕ_39,2, ϕ_40,2, ϕ_41,2
and ϕ_42,2 solutions are shown in Fig. <ref>
(solutions with folds on the boundaries are also shown). All the two-fold
solutions are unstable displaying a single negative eigenvalue. At
the bifurcation point to a two-fold solution (e.g. at P≈1.98
for ϕ_39,2) the amplitude eigenvalue inherited from the periodic
state vanishes but becomes positive as P decreases just as in the
single fold case. However, this eigenvalue is passed by that of another
mode in the vicinity of the bifurcation point (e.g. at P⪅1.975
for ϕ_39,2). The corresponding eigenmode is non-periodic and
remarkably interesting. For the ϕ_41,2 and ϕ_42,2
solutions, these two lowest eigenvalue modes display the same modulation
as the periodic solution but with a ±π/2 shift in the
carrier wave. Because of the localization of the folds, this feature
is a reminiscent of a Goldstone translation mode. The eigenmode acts
on the solutions by translating the left fold to the right (+π/2
shift) and the right fold to the left (-π/2 shift), or
vice versa, without changing their shape. Consequently we refer to
this mode as an translational attraction (or repulsion) eigenmode.
The other lowest eigenvalue mode is slightly different: on a one side
of a fold the mode displays a π shift while on the other side
there is no phase shift. These eigenmodes also generate attraction
(or repulsion) between folds but also change their shape. We refer
to them as rolling attraction (or repulsion) eigenmodes. While in
the translation mode both the carrier wave and the envelope translate,
in the rolling case only the envelope shifts, resulting in a local
deformation of the folds. The situation for ϕ_39,2 and ϕ_40,2
is similar with the difference that the translation and rolling eigenmodes
change their order. Figures <ref>(d) and
(e) show pseudo-solutions that indicate the tendencies represented
by these modes. We conclude that the two-fold states are in all cases
unstable and that the instability takes the form of a gradual approach
of the two folds, ultimately resulting in a stable single fold state.
The manifestation of the instability may be very slow, however, as
it depends on the energy difference between these two states. This
conclusion is supported by the fact that the unstable eigenvalue of
the two-fold state approaches zero from below, asymptotically exponentially,
as P decreases. This comes as no surprise since the two folds can
only interact via the overlapping tails of their profiles and these
are exponentially small when the folds are sufficiently far apart,
a fact that suggests that as the solutions increasingly localize (i.e.,
P decreases), the time required for the system to reach the single
fold state becomes exponentially long. This decay may arise through
attraction via either the translation or rolling eigenmodes. In contrast,
two-fold solutions with folds on the boundary decay into a single
fold state much more rapidly as the folds at the boundary decay leaving
only the fold at the center.





§ CONCLUSIONS


In this article, we reported on the existence of multifold solutions
in a thin floating elastic sheet of finite length with hinged boundary
conditions. Starting from known large amplitude periodic solutions
<cit.>, we used numerical continuation to uncover a
series of solution branches that bifurcate from the periodic (wrinkle)
branch as a function of the load P or the associated compression
Δ. Each branch consists of an array of localized (fold) solutions.
Symmetric and antisymmetric single fold solutions, which are well
known for the infinite sheet case <cit.>, are found
to bifurcate from branches of distinct periodic states and inherit
their symmetry properties. The multifold states with m=2,3,…
folds in the domain bifurcate from the periodic states in subsequent
bifurcations and are, we believe, new.

Based on our numerical continuation results, we performed a weakly
nonlinear analysis of the multifold solutions. These states form as
a result of a slow spatial modulation of the periodic wrinkle solutions.
The spatial modulation satisfies an amplitude equation whose solutions
can be written down in terms of Jacobi elliptic functions. The results
determine the branches of multifold solutions in the parameter space
and match the numerical continuation results remarkably well. In particular,
we determined the energy of multifold states and showed that the energies
of like-number multifold states are almost equal, with energy gaps
that become exponentially small as the sheet length increases (see
Eq. (<ref>)).

Finally, we also studied the stability of both the periodic and multifold
solutions within a variational framework. We showed numerically that
the symmetric and antisymmetric single fold solutions are both stable,
but that only one of them bifurcates from a stable periodic branch
and so inherits its stability. We also showed that all the two-fold
solutions are unstable, but with a common instability mode characterized
by the mutual attraction/repulsion of folds (see Sec. <ref>).
However, the relevant eigenvalue comes very close to zero as the localization
increases, implying that the transition to a single fold state requires
essentially exponentially long times. This result is in agreement
with the numerical simulations in <cit.> (with clamped
boundary conditions) in which a three-fold state is observed to decay
into a single fold state but did so so slowly that it was thought
to be stable.

These results have been obtained for the special case of hinged boundary
conditions. However, some numerical experiments employ clamped boundary
conditions, see e.g. <cit.>, and in this case the
route to localization appears to be different. Thus an extension of
the type of analysis performed here to the case of clamped boundary
conditions ϕ=ϕ”=0 at s=± L/2 would also be valuable.

Our results can be used to provide a rule of thumb for experiments
employing hinged boundary conditions where localization of wrinkles
into folds is crucial. Multifold solutions may emerge naturally in
certain geometries (for instance, in a floating elastic annulus or
an indented circular sheet <cit.>).
The prediction of the energy gap between multifold and single fold
states obtained here can be used to determine which states can emerge
deterministically or stochastically in the presence of thermal noise.
This is of particular importance in microscopic-scale experiments
<cit.> and in experiments in which the
floating raft is granular <cit.>. Although in
most experiments folds are allowed to form spontaneously, our calculations
also provide a basis for predicting the outcome of experiments with
manually induced folds, which may be of particular interest in technological
applications. 

This work was partially funded by the Berkeley-Chile Fund for collaborative
research and CONICYT-USA PII20150011. L.G. wishes to acknowledge support
by Conicyt/Becas Chile de Postdoctorado 74150032 and Conicyt PAI/IAC
79160140. We are grateful to Enrique Cerda, Marcel Clerc and Punit
Gandhi for fruitful discussions. 







§ JACOBI ELLIPTIC INTEGRALS


The quantities related to A in Sec. <ref>
can be written in terms of fundamental Jacobi elliptic integrals,

    (A _^2,A _^2)     =     2^4κ^2(D_2,C_2),
    (A_^'2,A_^'2)     =     2^4κ^4(D_2^',C_2^'),
    (A _^4,A _^4)     =     2^8κ^4(D_4,C_4),

where the coefficients D_2,C_2,D_2^',C_2^',D_4,C_4
are given, respectively, by the following mean values of the Jacobi
elliptic functions: 

    [                             D_2                               ≡                            _k^2                               =                          1-I(k);                           D_2^'                               ≡                           _k^'2                               =            1/3[k^2-(2-k^2)I(k)];                             D_4                               ≡                            _k^4                               =       1/3[(3-k^2)-2(2-k^2)I(k)];                             C_2                               ≡                         k^2_k^2                               =                        k^2-I(k);                           C_2^'                               ≡                        k^2_k^'2                               =           1/3[k^2-(2k^2-1)I(k)];                             C_4                               ≡                         k^4_k^4                               = 1/3[k^2(3k^2-1)-2(2k^2-1)I(k)]. ]




§ EXPANSION OF FUNCTIONALS


Consider a small perturbation of a solution y_0 satisfying the
variational problem, y=y_0+δ y, where δ y is assumed
to be of order ε, i.e. δ y=ε y_1≠0
and y_1∼𝒪(1): 

    y     =     y_0+ε y_1.

The function ϕ(s) can be expanded in terms of ε,

    ϕ=ϕ_0+εϕ_1+ε^2ϕ_2+…

and likewise for the load P, 

    P     =     P_0+ε P_1+…

and the compression and energy functionals, 

    Δ    =    Δ_0+εΔ_1+ε^2Δ_2+…
    
    E     =     E_0+ε E_1+ε^2E_2+…
    
    G     =     G_0+ε G_1+ε^2G_2+…

Each term in these expansions can be obtained from the definitions
of these relevant quantities (Sec. <ref>).
For instance, the function ϕ is related to y through ẏ=sinϕ,
and it can be shown that the first and second order terms in the expansion
of ϕ are, respectively, 

    ϕ_1    =    ϕ_0 ∂_sy_1,
    ϕ_2    =    1/2tanϕ_0^2ϕ_0(∂_sy_1)^2.

Analogously, the functionals Δ and E can be written as

    Δ    =Δ_0+ε(δΔ/δ y,y_1)+ε^2(y_1,δ^2Δ/δ yy_1)
    
    E     =E_0+ε(δ E/δ y,y_1)+ε^2(y_1,δ^2E/δ yy_1),

where Δ_0 and E_0 and their functional derivatives
are evaluated at the solution y_0. For the sake of simplicity,
we only consider here the terms that are relevant for the stability
analysis. For the displacement, we obtain: 

    δΔ/δ y    =tanϕ_0D_s,
    δ^2Δ/δ y^2    =D_s^†Σ_s^ΔD_s,

while for the energy, 

    δ^2E/δ y^2    =D_ss^†Σ_ss^ED_ss-D_s^†Σ_s^ED_s+Σ_0^E,

where D_s and D_ss are, respectively, the first and second
order linear differential operators and D_s^† and D_ss^†
are the corresponding adjoint operators. The scalar operators Σ_s^D,
Σ_ss^E, Σ_s^E and Σ_0^E are given
by 

    Σ_s^Δ≡   1/2^3ϕ_0,
    Σ_0^E≡   1/2^2ϕ_0(y_0ϕ̇_0+cosϕ_0),
    Σ_s^E≡   1/4y_0^2^3ϕ_0+…
       [ϕ̇_̇0̇^2(1/2+tan^2ϕ_0)+ϕ̈_0tanϕ_0]^2ϕ_0,
    Σ_ss^E≡   1/2^3ϕ_0.

A straightforward calculation shows that the free-energy second-order
term G_2 is given by 

    G_2=(y_1,δ^2G/δ yy_1)-P_1(δΔ/δ y,y_1)

, where the second order functional derivative is simply 

    δ^2G/δ y^2=δ^2E/δ y^2-P_0δ^2Δ/δ y^2.

Notice that the functional derivatives have been introduced in symmetrical
form, which simplifies the numerical calculations of Sec. <ref>.

 apsrev


