Eur. Phys. J. E 18, 343–358 (2005)
DOI: 10.1140/epje/e2005-00038-5

THE EUROPEAN
PHYSICAL JOURNAL E

Buckling of spherical shells adhering onto a rigid substrate
S. Komuraa, K. Tamura, and T. Kato
Department of Chemistry, Faculty of Science, Tokyo Metropolitan University, Tokyo 192-0397, Japan
Received 28 July 2005 / Received in ﬁnal form 2 October 2005
c EDP Sciences, Società Italiana di Fisica, Springer-Verlag 2005
Published online 15 November 2005 – 
Abstract. Deformation of a spherical shell adhering onto a rigid substrate due to van der Waals attractive
interaction is investigated by means of numerical minimization (conjugate gradient method) of the sum of
the elastic and adhesion energies. The conformation of the deformed shell is governed by two dimensionlessparameters, i.e., Cs / and Cb / where Cs and Cb are respectively the stretching and the bending constants,and  is the depth of the van der Waals potential between the shell and substrate. Four diﬀerent regimesof deformation are characterized as these parameters are systematically varied: (i) small deformation
regime, (ii) disk formation regime, (iii) isotropic buckling regime, and (iv) anisotropic buckling regime. Bymeasuring the various quantities of the deformed shells, we ﬁnd that both discontinuous and continuous
bucking transitions occur for large and small Cs /, respectively. This behavior of the buckling transition
is analogous to van der Waals liquids or gels, and we have numerically determined the associated criticalpoint. Scaling arguments are employed to explain the adhesion induced buckling transition, i.e., fromthe disk formation regime to the isotropic buckling regime. We show that the buckling transition takes
place when the indentation length exceeds the eﬀective shell thickness which is determined from theelastic constants. This prediction is in good agreement with our numerical results. Moreover, the ratio
between the indentation length and its thickness at the transition point provides a constant number (2–3)
independent of the shell size. This universal number is observed in various experimental systems ranging
from nanoscale to macroscale. In particular, our results agree well with the recent compression experimentusing microcapsules.
PACS. 46.32.+x Static buckling and instability – 68.35.Np Adhesion – 81.05.Tp Fullerenes and relatedmaterials

1 Introduction
1.1 Elastic sheets
Investigations of structures and properties of thin elastic
sheets are important for both practical and industrial reasons. Their applications in our daily life are such as cans,houses, domes, bridges, ships, planes, etc. When elastic
sheets are subjected to a large external force, they loose
their shape and buckle at a critical force as we often experience. In engineering of safety structures, it is desired
to increase the buckling threshold as large as possible. As
discussed below, one of the characteristic features of thinelastic sheets is that the energy required for stretching isvery large compared to that for bending. Hence pure bending deformations without any stretching are preferred ingeneral.
In recent years, considerable attentions have been paid
to thin materials which exhibit elastic behaviors at either microscopic or mesoscopic scales. For example, theconformation of two-dimensional sheets of graphitic oxideae-mail: komura@comp.metro-u.ac.jp

was investigated by electron microscopy and other techniques [1,2]. A sheet of graphitic oxide is ﬂat on averagedue to its ﬁnite in-plane shear elasticity in spite of outof-plane thermal ﬂuctuations [3]. For a two-dimensional
polymer network (polymerized silane monolayer) at the
air-water interface, a buckling phenomenon was observedby X-ray scattering [4]. The other group found a continuous buckling transition for a solid Langmuir monolayercomposed of phospholipid molecules deposited on the surface of formamide [5,6]. A more complex elastic sheet canbe found in a cell membrane skeleton called “cytoskeleton” which is a two-dimensional triangulated network consisting mainly of actin and spectrin molecules. Such kindof biological membranes take locally rough but globallyﬂat conﬁguration even in the presence of thermal ﬂuctuations [7]. Deﬂection of a cell membrane under a localized
force or torque was discussed theoretically [8].
From the theoretical point of view, the properties ofstretching ridges in a crumpled elastic sheet has been studied intensively during the past decade. The crumpling ofa thin sheet can be understood as the condensation ofelastic energy into a network of ridges [9,10]. The problem of ridges was initiated by Witten and Li who realized

344

The European Physical Journal Ethat the elastic energy of a ridge scales as X 1/3 , where
X is the length of the ridge [11,12]. Lobkovsky re-derived
the same ridge scaling relation by performing a boundary
layer analysis of the Föpple-von Kármán (FvK) equations
for plates [13,14]. Note that the term “plate” means athin sheet of elastic material of constant thickness whichis plane in its rest state. More recently, DiDonna andWitten investigated the compression and buckling of elastic ridges, and showed that the energy required to buckle
it is nearly a constant fraction (20%) of the total ridge
elastic energy [15,16]. In fact, this explains why crumpled sheets are qualitatively stronger than smoothly bent
sheets [17].
In general, the highly nonlinear FvK equations are very
diﬃcult to solve because they involve two types of deformation (stretching and bending) with energies of diﬀerent
orders of magnitude [18]. Nevertheless, the FvK equations
have been the subject of renewed interest in the context
of developable cone (d-cone) singularities. Surprisingly, ageometry of a d-cone is one of the solutions to the complicated FvK equations [19]. Some scaling relations for thecore size of the d-cone singularity was obtained in references [19,20]. Being inspired by these theoretical predictions, several macroscopic experiments were performed toinvestigate the shape, response, and stability of the d-conesingularities [21–23].
Other than the d-cone geometry, Audoly et al. predicted various buckling modes of a long rectangular elastic
plate subjected to the applied longitudinal and transverse
compressions [24–26]. Reference [27] reports the case inwhich a plate was initially bent in one direction into acylindrical arch, and then deformed in the other direction. Later the core energy of the d-cone singularity wasmeasured by piercing the plate around the singularity [28].
More recently, wrinkling of an elastic sheet was discussed
by Cerda and Mahadevan [29,30].

Conversely, the stretching energy becomes dominant assoon as the displacement ζ is larger than the thickness h.
For a spherical shell of radius R, on the other hand, therespective energies are given byEs ∼ Y hζ 2 /R2 ,Eb ∼ Y h3 ζ 2 /R4 .

(2)

Then the similar ratio for a shell becomes Es /Eb ∼(R/h)2 , which is typically very large. Hence the bending
deformation inevitably accompanies the stretching deformation for a shell. In contrast to the case of a plate, thisfact does not dependent on the ratio ζ/h of the shell.
As an example of this interplay between the stretching
and bending modes, the shape ﬂuctuations and the stability of a cylindrical shell (polymerized vesicle) was studied
before [31]. It was shown that the intrinsic curvature of
the shell leads to an enhanced coupling between the twoelastic modes, and act to suppress the shape ﬂuctuations
on large scales. A similar analysis for a spherical shell revealed that such a suppression eﬀect is more pronounced
when shells are closed [32–34]. The enhanced stability ofshells is crucial for constructing large structures such as
domes or bridges.
There are some experimental works which deal withthe deformation of spherical shells. By using actin-coatedvesicles [35], a buckling instability was observed whena large localized force is applied [36]. At macroscopic
level, on the other hand, contact and compression problem of a ping-pong ball was investigated by Pauchard
and Rica [37,38]. The same author reported that buckling instability occurs during the drying of sessile dropsof polymer solution [39]. Similar phenomena were foundalso by using droplets of colloidal suspensions [40,41].
Rather recently, elastic properties of polyelectrolyte capsules [42,43] are studied by AFM and reﬂection interference contrast microscopy [44–46]. The details of these experimental works will be discussed in Section 6 in order
to compare with our results.

1.2 Shells
It should be stressed, however, that most of the above
works are concerned with elastic sheets which are ﬂat in
the undeformed state, while less attention has been paid
for the properties of initially curved elastic sheets, i.e.,shells. Moreover, the FvK equations are valid only fora plane geometry, and a covariant generalization for any
curved shell is still lacking. In particular, a shell exhibits apeculiar elastic feature because the strain tensor is proportional to the ﬁrst power of the out-of-plane displacement(see Sect. 2 for the details), and the shell cannot be bent
without being stretched [18]. For a plate, typical stretching and bending energies per unit area scale asEs ∼ Y hζ 4 /4 ,Eb ∼ Y h3 ζ 2 /4 ,(1)

respectively, where Y is the Young’s modulus, h the thickness of the elastic sheet, ζ the magnitude of the out-ofplane displacement, and  is a typical length scale. Since
the ratio between the two energies is Es /Eb ∼ (ζ/h)2 , thestretching energy can be neglected in the limit of ζ  h.

1.3 Adhesion
Among various types of deformation, adhesion onto a substrate due to van der Waals interaction plays an importantrole especially in the ﬁeld of nanotechnology. For example,
the electric transport through carbon nanotubes is studied after their deposition on a substrate with which theyinteract each other. Unfortunately, it is known that theresistivity of the nanotube is aﬀected by its elastic deformation. Since there is little control over the alignment andthe shape of adsorbed nanotubes, it is crucial to know how
they deform on the substrate. The deformations of multiwalled nanotubes on a rigid substrate was observed and investigated using atomic force microscopy (AFM) [47] andmolecular-mechanics simulations [48]. Later a collapse of
a nanotube section due to the surface interaction was observed by using AFM [49]. More recently, a systematic numerical study on the deformation of an elastic nanotube
adhering onto a substrate was reported by the present authors [50]. However, we stress here that the adhesion and

S. Komura et al.: Buckling of spherical shells adhering onto a rigid substratecontact problem of a spherically closed shell has not yet
been investigated in detail. One exception is reference [51]in which some scaling arguments for the deformation andmechanical stability of fullerene-like hollow nanoparticles
were given by employing the shell theory. It was shownthat van der Waals interactions between a substrate andadhering nanoparticles can cause considerable deformations. Such an eﬀect is important for tribological applications of fullerene-like nanoparticles [52].
1.4 Present work
In this paper, we investigate both numerically and theoretically the deformation and the stability of a sphericallyclosed shell adhering onto a rigid substrate due to vander Waals attractive interaction. Some of the results havebeen published elsewhere [53]. To consider the shell adhesion, we propose a discretized model in which the equilibrium conﬁguration of the shell is determined according
to the competition among three energies, i.e., stretching,bending, and van der Waals energies. The total energy
is numerically minimized by using the conjugate gradient method. The eﬀects of thermal ﬂuctuations are nottaken into account in the present model. One of our mainﬁndings is the fact that the adhesion causes a bucklingtransition of a spherical shell as either the elastic properties and/or the strength of adhesion is varied. More interestingly, our systematic study revealed that the buckling transition can be either continuous or discontinuous
depending on the elastic properties of the shell such as
the stretching or bending constants. We identify a special point which is analogous to the critical point for vander Waals ﬂuids which exhibits liquid-gas coexistence. We
have also performed the scaling analysis to explain theadhesion induced buckling of spherical shells. Once thebuckling occurs, a polygonal structure consisting of ridges
and d-cones is created when the adhesion is strong enough.
Our work can be regarded in part as a contact problem of spherical shells. For two elastic bodies, their contact
problem under an applied load was solved by Hertz longtime ago [18,54]. After ninety years, the Hertz’s solution
was extended to take into account the inﬂuence of adhesion energy [55,56]. In the presence of the adhesion energy,
the apparent load acting between the two elastic bodies
is larger than the applied load. In contrast to the contact
problem of elastic bodies, there are few works which deal
with the corresponding problem of elastic shells.
This paper is constructed as follows. In the next section, we brieﬂy review the framework of the continuum
elasticity theory for shells. In Section 3, we describe ourmodel for shells adhering onto a rigid substrate. We also
explain the numerical method to calculate the equilibrium
structure. Then we present the obtained results together
with various quantitative analyses of the shell structuresin Section 4. In Section 5, we provide some scaling arguments concerning the deformation of shells, and compare them with our numerical results. Finally, the paper
is closed with discussions in Section 6 where we compare
our results with several previous experiments.

345

2 Shell theory
In this section, we describe the continuum version of theshell theory [18]. We collect some formulas from diﬀerential geometry which is the most appropriate formalism forthe classical theory of elastic shells. See reference [57] fora further treatment.
One can, in general, parameterize a two-dimensional
thin sheet in three-dimensional space by two real innercoordinates s = (s1 , s2 ). The shape of the sheet is thendescribed by a three-dimensional vector r = r(s). At eachpoint on the sheet, there are two tangent vectors ri =∂r/∂si with i = 1, 2. The outward unit normal vector n̂
is perpendicular to these tangent vectors, i.e.,n̂ =r1 × r2·|r1 × r2 |(3)

All properties related to the intrinsic geometry of the sheet
are expressed in terms of the metric tensor (or the fundamental tensor) deﬁned by
gij = ri · rj .

(4)

Two important quantities are the determinant and theinverse of the metric tensor which will be denoted byandg = det(gij ),(5)

g ij = (gij )−1 ,(6)respectively. In addition, one has to consider the (extrinsic) curvature tensor (or the second fundamental tensor)
given byhij = n̂ · ∂j ri = −n̂i · rj ,(7)

with ∂j ri = ∂ 2 r/∂si ∂sj . Note that a surface is uniquely
characterized by its metric tensor gij and the curvaturetensor hij .
Let us deﬁne the (undeformed) reference state as r =R. Here and below, we shall use capital letters in order to distinguish quantities in the reference state fromthe corresponding quantities in the deformed state. Thus,
Ri , N̂i , Gij , Hij represent the tangent and normal vectors, the metric and the curvature tensors in the reference
state, respectively. If the sheet is stretched, the distance
between two neighboring points in the sheet is changed.
This change can be expressed in terms of the strain tensor
uij deﬁned by [58]uij =1(gij − Gij ).
2(8)

The mixed strain tensor is obtained by raising one of theindices according to
ui j = uik g kj .

(9)Here and below, we use Einstein’s summation convention
and sum over all indices which appear twice. Likewise, themixed bending tensor
bij = hij − Hij

(10)

346

The European Physical Journal E

is taken as a measure for the bending deformation. Notethat this choice is not unique and alternative deﬁnitions
are possible. We will discuss this point later.
Consider a deformation of a shell which can beparametrized by
r = R + ui Ri + ζ N̂.

(11)

The contravariant vector ui (i = 1, 2) represents lateral(in-plane) displacement ﬁeld and ζ represents the transverse (out-of-plane) displacement ﬁeld. Both strain tensor
and bending tensor can be expressed in terms of the components of these displacement ﬁelds. Up to ﬁrst order inthe displacement r − R, the strain tensor turns out to beuij ≈1(Di uj + Dj ui ) − ζHij ,2(12)

where Y and ν are the (three-dimensional) Young’s modulus and the Poisson’s ratio, respectively. Then the elastic
energy density is given by 1/2A
1 h/2 ij
σ ηij
dz,
(18)
fe =
2 −h/2G
where h is the thickness of the shell and
A=

(19)

According to the shell theory, fe can be expressed in termsof uij and bij by constructing the possible invariants withrespect to coordinate transformations. Within the linear
elasticity, the deformation energy of an isotropic shell isgiven by [57]fe =where the covariant components of the lateral displacement ﬁelds are given by uj = Ai gij , and the covariant
derivative Di is deﬁned by
Di uj = ∂i uj − Γijk uk ,1[det(Gij − Hij z)]2 .
G

(13)

Yh
[νui i uj j + (1 − ν)ui j uj i ]2(1 − ν 2 )Y h3[νbi i bj j + (1 − ν)bi j bj i ].
+24(1 − ν 2 )(20)

From the principle of virtual work, we havewith the Christoﬀel symbol
j
Γik
=

1 jl
G (∂i Gkl + ∂k Gil − ∂l Gik ),
2

δW =

(14)

and ∂i = ∂/∂si . In a similar manner, the linear approximation for bij leads tobij ≈ Di Dj ζ − ζHi k Hjk + Hik (Dj uk )+ Hjk (Di uk ) + uk (Dj Hik ). (15)

where

∂W
∂W
δuij +
δbij ,
∂uij
∂bij

∂W
= N ij ,∂uij∂W= M ij
∂bij

(21)

(22)

are eﬀective membrane stress tensor and eﬀective moment
tensor, respectively. Then we can derive the constitutive
equations for shells:

In the Donnel-Mushtari-Vlasov approximation [57], thebending tensor is simpliﬁed to

N ij =Yh[νg ij uk k + (1 − ν)uij ],(1 − ν 2 )(23)

bij ≈ Di Dj ζ.

M ij =Y h3[νg ij bk k + (1 − ν)bij ].
24(1 − ν 2 )

(24)

(16)

As mentioned, we have deﬁned the bending tensor bij
in terms of the diﬀerence between the covariant components of the curvature tensors in the deformed andthe initial state (see Eq. (10)). If we had selected, forinstance, the diﬀerence between the mixed components,b̃i j = hi j − Hi j , or the contravariant components, b̄ij =hij − H ij , the results would be diﬀerent. This fact is ofgreat importance to the shell theory, although the term
Di Dj ζ in equation (15) is not aﬀected by the diﬀerent
measures of bending. From a purely formal point of view,
we may take any measure of strain and bending fromwhich the original measures can be recovered. There is,in fact, no physical ground for selecting one of the bending tensors.
Let σ ij be the components of three-dimensional stressin normal coordinates and η ij corresponding componentsof three-dimensional strain. For a bulk material that obeys
Hooke’s law, the stress tensor is given in terms of the strain
tensor by [57]


Y
ν
ij
ij
ij
k
σ =
g ηk
,
(17)
η +
1+ν
1 − 2ν

Using the identity det(bi j ) = 12 (bi i bj j − bi j bj i ), one
can rewrite the last two terms in equation (20) asfe =Yh[νui i uj j + (1 − ν)ui j uj i ]2(1 − ν 2 )κ+ bi i bj j + κ̄det(bi j ),
2

with
κ=

Y h3
,
12(1 − ν 2 )

(25)

(26)

Y h3
·
(27)
12(1 + ν)
Since bij has been deﬁned by the diﬀerence between two
curvature tensors, it has to satisfy certain compatibility
conditions in order to be itself a curvature tensor. In suchcase, κ and κ̄ are called bending rigidity and Gaussian
curvature modulus, respectively. Hence, in general, theGauss-Bonnet theorem will not apply to equation (25).
For a planar reference state, however, the Gauss-Bonnet
theorem applies.
κ̄ = −

S. Komura et al.: Buckling of spherical shells adhering onto a rigid substrate347

into account by using the model of polymerized membranes with a ﬁnite bending constant [9,62–65];
Eb =

1
αβ

Fig. 1. Beads and springs model for an elastic spherical shell
adhering onto a substrate. n̂α(β) is the unit normal vector of
the triangle α(β).

3 ModelConsider an elastic spherical shell interacting with a rigid
substrate as shown in Figure 1. The normal directions tothe substrate is taken as the z-axis, whereas the substratespans the xy-plane. In the discretized model, the conﬁguration of the shell is represented by a triangular mesh as asimplest approximation for a two-dimensional elastic material. In the absence of adhesion, the initial conﬁgurationof the shell is taken to be spherical. The initial conﬁguration is constructed by the Delaunay triangulation of thespherical surface [59,60]. Starting from an icosahedron as
the original network, we add new points on each triangle followed by a subsequent rescaling of all bonds to thedesired length [61,62]. Although there are always 12 grid
points which have ﬁve neighbors, this procedure ensuresthat most of grid points have six neighbors and each bondhas approximately the same length. The singularity associated with the ﬁvefold symmetry will be discuss later in
Section 6. In the present work, we studied shells consisting of N = 10 × 3k + 2 grid points with k = 1, . . . , 5. The
number of triangles is f = 2N − 4 while the number ofbonds is f = 3N −6. These quantities certainly satisfy the
Euler’s theorem; N + f − e = 2. In the next section, we
mainly present the results for k = 4, i.e., N = 812. The
size eﬀect will be separately discussed there. Hereafter weassociate all the grid points and bonds with beads andsprings, respectively.
To describe the deformation of an elastic shell, boththe stretching and the bending energies should be taken
into account [18,57]. Following the model of membranes
with crystalline order [63], or crushed elastic manifolds [9],the discretized stretching energy is given by the sum over
Hooke’s law of each spring:
Es =

1
n

2


Cs

Ln − L0
L0

2
·

(28)

Here Cs is the stretching (spring) constant, Ln is thelength of spring n, and L0 is the natural length of thespring (or the lattice constant) taken here as a constant.
On the other hand, the discretized bending energy is taken2

Cb |n̂α − n̂β |2 ,(29)

where Cb is the bending constant, n̂α(β) is the unit normal vector of triangle α(β), and the sum is taken over eachpair of triangles which share a common edge. The bending constant Cb plays the role of a Heisenberg exchange
coupling between neighboring normals. We note here thatboth Cs and Cb have the dimension of energy. The
comparison between the continuum elasticity theory andthese discretized elastic energies will be discussed later in
Section 5.
To calculate the adhesion energy of the shell, we consider a generalized Lennard-Jones type interaction that is
acting between each of the bead and the substrate: 28/3  σ 12  σ 3
W =
−
,
(30)
3
zi
zi
i
where zi is the height of bead i from the substrate, andthe sum is taken over all the beads. When the adhesionenergy of bead i is plotted against the distance zi , thedepth of the energy minimum is given by , and the distance corresponding to this minimum is 22/9 σ. The ﬁrst
repulsive term in equation (30) is responsible for the excluded volume interaction which prevents the beads from
penetrating into the substrate. The power of this repulsive
potential should not necessarily be 12, and a lower powersuch as 9 can be used as well. For our numerical calculations, it is more suitable to employ a stronger repulsive
potential.
The second term represents the long-ranged attractive
interaction between the beads and the substrate. The inverse cubic dependence of the above potential is brieﬂy
explained below [56,66]. The van der Waals attractive interaction between two atoms is generally given by the formCv(r) = − 6 ,r(31)

where C is a constant depending on the physical originof the attraction, and r is the distance between the two
atoms. The simplest approach to obtain the interaction
between an atom and a macroscopic body such as a substrate is to sum up the interactions between all pairs of
atoms [56,66]. We consider a case where a single atom
is placed at a distance D from a semi-inﬁnite medium of
density ρ. Then the total interaction energy is given by
w(D) = −

2πρC
,
12D3

(32)

which gives rise to the inverse cubic dependence of thepotential.
The total energy
Etot = Es + Eb + W

(33)

348

The European Physical Journal E

Table 1. The numerically obtained values of the critical point (Cs /)c , (Cb /)c , (H/R)c for the sizes N = 92, 272, 812, and2432. We also list the values of the input parameters L0 /σ, R/σ, and R/L0 . The values of (h/R)c at the critical point arecalculated by using equation (46). (H/h)c is the ratio between (H/R)c and (h/R)c . These numbers are almost independentof N .
N
92
272
812
2432

L0 /σ
0.1
0.1
0.1
0.1

R/σ
0.25
0.43
0.75
1.3

R/L0
2.5
4.3
7.5
13

(Cs /)c
5
80
210
500

is numerically minimized using the conjugate gradientmethod [67]. Like most methods of multidimensional minimization, it is performed as a series of one-dimensional
minimizations. For this purpose, a series of noninterfering,conjugate directions are constructed. As a result, minimization along one direction does not disturb the minimization in the other conjugate directions. Hereafter all
the energies and the lengths are respectively scaled byand σ which characterize the shape of adhesion interaction
in equation (30). There are three independent dimensionless parameters in the model, i.e., Cs / , Cb / , and L0 /σ.
In the present study, we have mainly varied Cs / and Cb / ,whereas the other parameters are ﬁxed. The scaled natural length of each bond L0 /σ is chosen such that the initial spherical conﬁguration of the shell does not store any
stretching energy Es , and its value is roughly L0 /σ ≈ 0.1.
(Notice that not all the bonds have exactly the same natural length because of the singularity associated with theﬁvefold symmetry.) For N = 812, the scaled radius of the
undeformed shell is R/σ ≈ 0.75 (see Tab. 1). We note thatthe bending energy is inherent even in the undeformed
spherical shell since the spontaneous curvature is not included in the present calculation. This assumption is justiﬁed such as for fullerene balls. (See also discussion in
Sect. 6 concerning the spontaneous curvature of shells.)
In our model, the eﬀect of thermal ﬂuctuation is notincluded. Hence our calculation corresponds to the zerotemperature numerical simulation. Since the excluded volume eﬀect of the surface is not included, we are dealing
with “phantom” shells. As we shall see later, the neglectof self-avoidance eﬀect is justiﬁed for most of the moderate deformations even when the buckling takes place.
Self-avoidance can be crucial such as when the shell collapses due to a large negative pressure [51,62,68].

4 Results
In this section, we collect and present our numerical results which are analyzed by various quantitative methods.
Some of them have been already published elsewhere [53].
We mainly discuss the results from the size N = 812. The
size dependence is discussed in Section 4.6.
4.1 Configurations
By looking at various equilibrium conﬁgurations of the deformed shells, we ﬁnd that there are basically four qualitatively distinct types of deformation as the combination(Cb /)c
0.35
1.5
6.1
22

(H/R)c
1
0.26
0.18
0.13

(h/R)c
0.29
0.089
0.064
0.046

(H/h)c
3.4
2.9
2.8
2.8

Fig. 2. Top, side and bottom views of the equilibrated conﬁgurations of adhering spherical shells when the sets of the scaled
elastic constants (Cs /, Cb /) are (a) (1000, 1000), (b) (150, 9),(c) (150, 2), and (d) (100, 1).

of Cs / and Cb / is systematically varied. Typical examples for these cases are shown in Figure 2 from (a) to(d). For a given parameter set, the deformed shell is seenfrom top, side, and bottom with respect to the substrate.
Figure 2a (Cs / = 1000, Cb / = 1000) corresponds to thesituation when both of the elastic constants are very large
compared to the adhesion energy. Here the shell hardly
deforms in spite of the adhesion, and keeps its sphericalshape (“small deformation regime”). This means that theeﬀect of adhesion is practically irrelevant.
When both of the elastic constants are simultaneously
decreased, we observe the case Figure 2b (Cs / = 150,

S. Komura et al.: Buckling of spherical shells adhering onto a rigid substrate

Cb / = 9). Here a ﬂat contact disk develops at the bottomof the shell as can be observed from the side view (“diskformation regime”). The area of the contact disk increases
as the adhesion energy becomes larger. The formation of
a ﬂat contact disk is reminiscent of the ﬂattening of elastic
tubes along the contact region [50].
Keeping the value of Cs / while decreasing the bending constant Cb / results in the buckling of the shell as
illustrated in Figure 2c (Cs / = 150, Cb / = 2). For sucha buckled conﬁguration, both the stretching and the bending energies are localized at a narrow “bending strip” ofcontact. The competition between the two energies determines its width, which will be discussed later in Section 5
within the scaling argument. The bending strip is formed
in a circular shape, and the whole conﬁguration of the shellis almost isotropic in the xy-direction (“isotropic buckling
regime”). The buckled region of the shell is bent inward,
but it does not violate the excluded volume eﬀect.
In the case of Figure 2d (Cs / = 100, Cb / = 1), it becomes energetically favorable for the buckled region to create a polygonal structure composed of a number of ridges
jointed by the d-cones. In contrast to the isotropic buckling in Figure 2c, the shell buckles in an anisotropic manner (“anisotropic buckling regime”). To characterize sucha buckled shape is related to the problem of post-buckling.
The number of ridges is dependent on the van der Waals
adhesive energy, which will be argued below. We remind
here that neither isotropic nor anisotropic buckling has
never been observed for elastic tubes [50].
In order to see the sequence of the systematic deformation more clearly, we ﬁxed the stretching constant to
Cs / = 100 and varied the bending constant Cb / from 50
to 1.1. The top, side and bottom views of the equilibrium
conﬁgurations are arranged in Figure 3. Figure 3a belongs
to the small deformation regime, (b) to the disk formation
regime, (c) and (d) to the isotropic buckling regime, and(e) and (f) to the anisotropic buckling regime, respectively.
It is remarkable that, in the anisotropic buckling regime,
only a slight change in the value of Cb / causes a big
diﬀerence in the ﬁnal conﬁguration, namely, pentagonal,
square, or triangular polygonal ridges (see also Fig. 2d).
It is likely that the pentagonal shape appears due to thepresence of the ﬁve-handed bead that is ﬁrst attaching
to the substrate in our simulation. The inﬂuence of isolated beads of ﬁvefold symmetry on a sphere [69,70] will
be discussed separately in Section 6.
4.2 Asphericity
Next we characterize the shape of the deformed shells
more quantitatively. For this purpose, we ﬁrst calculate
the moment of inertia tensor deﬁned by [71,72]
1 
Iαβ =
(ri,α − rj,α )(ri,β − rj,β ),2N 2 i j(34)

where ri is the position of bead i, and α, β = x, y, z. The
sum is taken over bead positions in a given conﬁguration. The three eigenvalues of Iαβ are ordered according

349

Fig. 3. Top, side and bottom views of the equilibrated conﬁgurations of adhering spherical shells when the sets of the scaled
elastic constants (Cs /, Cb /) are (a) (100, 50), (b) (100, 10),(c) (100, 5), (d) (100, 2), (e) (100, 1.2), (f) (100, 1.1).

to magnitude λ1 ≤ λ2 ≤ λ3 . The directions of the principal axes are given by the eigenvectors corresponding tothese eigenvalues.
As a quantitative measure of the asphericity of the deformed shell, we have calculated the following three quantities [68]:
λ1
Γ1 =
,
(35)
λ3
∆=

λ21 + λ22 + λ23 − (λ1 λ2 + λ2 λ3 + λ3 λ1 ),(λ1 + λ2 + λ3 )2andS=(λ1 − λ̄)(λ2 − λ̄)(λ3 − λ̄)
,
2λ̄3

(36)

(37)

where λ̄ = (λ1 +λ2 +λ3 )/3 in S. The range of each value is0 ≤ Γ1 ≤ 1, 0 ≤ ∆ ≤ 1, or − 81 ≤ S ≤ 1. The value of Γ1 is

350

The European Physical Journal E

Fig. 4. The anisotropic factor Γ1 deﬁned in equation (35)
as a function of the scaled bending constant Cb / for Cs / =100, 300, 500, 700, and 900.

Fig. 5. The anisotropic factor ∆ deﬁned in equation (36)
as a function of the scaled bending constant Cb / for Cs / =100, 300, 500, 700, and 900.

unity when the shell is completely isotropic. Conversely, ∆vanishes when the shell is isotropic, while it deviates fromzero for an anisotropic conﬁguration. Negative S meansthat the shell is oblate, while it is positive when the shellis prolate. In Figures 4–6, we have respectively plottedΓ1 , ∆, and |S| as a function of Cb / for various diﬀerent
values of Cs / ranging from 100 to 900.
We begin by discussing Figure 4. For Cs / = 100, Γ1
decreases monotonically as Cb / becomes smaller. In thiscase, the buckling occurs at around Cb / ≈ 10 when theslope of the curve changes drastically. A similar buckling
behavior is observed both for Cs / = 300 and 500. Unlessthe shell is buckled strongly as in the anisotropic buckling
region, the largest and the second largest eigenvalues arealmost equal; λ2 ≈ λ3 . This means that the deformation
is isotropic in the xy-direction (see Figs. 3a to d).

Fig. 6. The absolute value of the anisotropic factor |S| deﬁned
in equation (37) as a function of the scaled bending constant
Cb / for Cs / = 100, 300, 500, 700, and 900.

In Figure 5, the buckling of the shell is manifested
in the sharp increase of ∆ as Cb / is decreased. There
are jumps of ∆ for larger values of Cs / , which indicatesthe occurrence of a discontinuous buckling transition. For
Cs / = 100, on the other hand, ∆ changes continuously
and a continuous buckling takes place. It is worthwhile tomention that ∆ attains its minimum at Cb / larger thanits threshold value of the buckling. For Cs / = 300, theminimum and the discontinuous jump of ∆ occur at diﬀerent Cs / , whereas they coincide with each other for larger
Cs / . In the continuum limit, ∆ should vanish for a perfectspherical shell. However, due to our ﬁnite discretization of
the shell, a small asphericity exists even in the initial undeformed conﬁguration, i.e., ∆ ≈ 3.1 × 10−5 for N = 812.
Since this value is slightly larger than the minimum valuesof ∆, the appearance of minima in Figure 5 can partiallybe an artifact of the discretization. Although the origin of
the minima in ∆ is not completely clear, we note that it
has nothing to do with the global deformation of the shell.
The measured value of S is always negative, and in
Figure 6, we have plotted the absolute value of S as afunction of Cb / . Negative S reﬂects the oblate pancakelike shape of the deformed shells. We see here that the
behaviors of ∆ (Fig. 5) and |S| (Fig. 6) are quite analogous. From Figures 5 and 6, we conclude that the buckling
can occur both in continuous and discontinuous manners.
4.3 Minimized energy
Here we look at the minimized energy of the deformed
shell. Figure 7 shows all the energies equations (28–30)
and the total energy Etot as a function of Cb / whenCs / = 900. In this case, the shell exhibits a discontinuous buckling transition at around (Cb / )∗ ≈ 3.5 as indicated by the dashed line. When the value of Cb / crosses
this critical value from above, both the stretching energy
Es and the bending energy Eb increases abruptly. TheseS. Komura et al.: Buckling of spherical shells adhering onto a rigid substrate

Fig. 7. The minimized total energy Etot / as a function of
Cb / for Cs / = 900. The three energies Es , Eb , and W contributing to Etot are also shown. The discontinuous bucklingtransition occurs at (Cb /)∗ ≈ 3.5 indicated by the dashed line.

losses in the elastic energies are compensated by the gainin the van der Waals energy W which decreases discontinuously at the transition point. In other words, the shellbuckles at the expense of the elastic energy when the adhesive force is strong enough. Interestingly, however, thediscontinuity in Etot around the transition point is verysmall.
For diﬀerent values of Cs / , the behaviors of eachenergy are qualitatively similar when the discontinuous
buckling occurs. However, the discontinuities in Eb and W
at the transition are smaller when Cs / becomes smaller.
To see this more clearly, we have plotted in Figure 8 thebending energy Eb / as a function of Cb / for various
Cs / as before. When the discontinuous buckling occurs
for Cs / ≥ 300, all the data fall onto a single curve in thelarge Cb / region. Moreover, Eb / is almost proportional
to Cb / because the shell deforms only slightly (small deformation regime). By contrast, the continuous bucklingtakes place when Cs / = 100 for which the data deviate
from others.
Figure 9 is a similar plot of the adhesion energy W/
as a function of Cb / . When Cs / is larger than 300, W/
is almost independent of Cb / in the unbuckled region,
but starts to decreases as the buckling takes place. For
Cs / = 100, however, W/ decreases continuously as Cb /is reduced.
4.4 Indentation length
To investigate the nature of the buckling transition in
more detail, we have measured the indentation lengths
H1 and H2 as deﬁned in Figure 10. In Figure 11, we plot
H1 /R as a function of Cb / for various Cs / ranging from100 to 900. Here the radius of the undeformed spherical shell is R = 0.75σ when N = 812 (see Sect. 4.6 or
Tab. 1). In accordance with the aforementioned discussion, H1 changes discontinuously at the transition point

351

Fig. 8. The scaled bending energy Eb / as a function of thescaled bending constant Cb / for Cs / = 100, 300, 500, 700, and
900.

Fig. 9. The scaled adhering energy W/ as a function of thescaled bending constant Cb / for Cs / = 100, 300, 500, 700, and900.

for larger Cs / , revealing the ﬁrst-order nature of thebuckling transition. This discontinuous buckling transition takes place between the disk formation regime andthe isotropic buckling regime (Figs. 2b and c). Hence thecontact region changes from a disk to a ring at the transition point.
In the same way, we plot H2 /R as a function of Cb /in Figure 12. The length H2 deviates from zero only if
the shell buckles for which the behaviors of H1 and H2
are almost identical to each other. This result indicatesthat the geometry of the buckled region is represented bya mirror image of the original undeformed shell.
Figure 13 shows the variation of the total indentation
length H deﬁned by
H = H1 + H 2 .

(38)

The discontinuous jump in H/R becomes smaller as Cs /is decreased, and ﬁnally vanishes at around Cb / ≈ 6.1

352

The European Physical Journal E

Fig. 10. Notation of the indentation lengths H1 and H2 in theisotropic buckling regime. The total indentation length is given
by H = H1 + H2 . R is the radius of the undeformed sphericalshell. r is the radius of the circular bending strip whose widthis denoted by w.

Fig. 11. The scaled indentation length H1 /R as a function of
the scaled bending constant Cb / for Cs / = 100, 300, 500, 700,
and 900.

Fig. 12. The scaled indentation length H2 /R as a function of
the scaled bending constant Cb / for Cs / = 100, 300, 500, 700,
and 900.

Fig. 13. The total indentation length H/R as a function of the scaled bending constant Cb / for Cs / =100, 300, 500, 700, and 900. The ﬁlled circle being locatedroughly at (Cb /, H/R, Cs /) = (6.1, 0.18, 210) indicates thepoint at which the discontinuity vanishes.

4.5 Scaling relation(the ﬁlled circle). The corresponding critical indentation
length is H1 /R ≈ 0.18, and that of the stretching constantis Cs / ≈ 210. Below this value of Cs / , the buckling occurs continuously rater than discontinuously.
We immediately note that Figure 13 is very reminiscent of the isotherms of non-ideal gases in the pressurevolume plane. Analogous to the liquid-gas coexistence region of van der Waals ﬂuids, the region of discontinuoustransition has been shaded in Figure 13. In the present
model, the parameter Cs / plays a role similar to thetemperature of van der Waals ﬂuids. Another similar phenomenon is the volume transition of gels which is induced
either by changing the temperature or the ionic strength.

Here we analyze the geometry of the buckled shell from adiﬀerent aspect. In Figure 14, we have plotted the relation
between the scaled ring radius r/R and the total indentation length H/R (see Fig. 10 and Eq. (38)) for various
combinations of Cb / and Cs / when N = 812. Diﬀerent points represented by the same symbol correspond to
diﬀerent Cb / values having the same Cs / values. Interestingly, most of the data collapse onto a single line, andwe ﬁnd that a scaling relation r/R ∼ (H/R)1/2 holds inthis regime. This scaling relation results from a simple
geometrical consideration. As we will discuss in the next
section, the buckled region is almost a mirror image of theoriginal undeformed shell.

S. Komura et al.: Buckling of spherical shells adhering onto a rigid substrate

Fig. 14. Relation between the scaled ring radius r/R deﬁned
in Figure 10 and the scaled total indentation H/R for various
combinations of Cs / (ranging from 100 to 500) and Cb /. Mostof the data collapse onto a single line which gives the scaling
relation r/R ∼ (H/R)0.5 .

The data deviate from a straight line when H/R issmall because the deformation of the shell cannot be described by the geometry in Figure 10. Moreover, the above
scaling relation does not hold in the anisotropic buckling regime in which the polygonal ridges are formed (seeFigs. 3e or f).
4.6 Size dependence
So far, we have discussed only the results when N = 812.
Even the shell size N is varied, the qualitative properties of the deformation are unchanged from the case of
N = 812. For example, the total indentation length H
behaves similarly to Figure 13 although the location of
the critical point shifts systematically. Table 1 summarizes the values of the critical point for four diﬀerent sizes
N = 92, 272, 812 and 2432. We have simultaneously listed
the values of L0 /σ, R/σ, and R/L0 for each size.
Roughly speaking, the critical values of the elastic constants (Cs / )c and (Cb / )c are larger for larger shells. On
the other hand, the values of (H/R)c decreases for larger
shells. This tendency holds true as long as the potential
range satisﬁes R/σ ≤ 1. When R/σ is much larger thanunity, the buckling does not occur. It is reasonable to thinkthat larger shells with a small curvature can be easily deformed due to the reduced coupling eﬀect.
4.7 Hysteresis
If the initial shape of the shell is far from a sphere, wesometimes could not ﬁnd the global minimum within the
conjugate gradient method. For complete spherical shells,
on other hand, we could always obtain reasonable equilibrium conﬁgurations as depicted before. For the initial

353

Fig. 15. The total indentation length H/R as a function of
the scaled bending constant Cb / for two diﬀerent initial conﬁgurations. The case (a) is the result when a spherical shell isused as the initial conﬁguration. In the case of (b), the strongly
buckled ﬁnial conﬁguration obtained when Cs / = 500 andCb / = 1 is used as the initial conﬁguration. There is a small
hysteresis.

spherical conﬁguration, the stretching energy Es is relaxedin the absence of adhesion although the bending energy
Eb is inherent. We remind that the eﬀect of spontaneous
curvature is not included in our model.
For certain parameter choices, however, the numerical results seem to depend on the initial conﬁguration. In
Figure 15, we have plotted the equilibrium total indentation length H/R obtained from the two diﬀerent initial
conﬁgurations but having the same elastic parameters.
The case (a) is the result when a spherical shell is used as
the initial conﬁguration. In the case of (b), the strongly
buckled ﬁnial conﬁguration obtained when Cs / = 500
and Cb / = 1 is used as the initial conﬁguration. Althoughmost of the results obtained from these two cases coincidewith each other, there is a slight diﬀerence in bucklingtransition point, which results in a small hysteresis. The
observed hysteresis becomes more remarkable for smaller
Cs / and/or Cb / , but we did not investigate it systematically since it is impossible to scan all the allowed initial
conﬁgurations.

4.8 Multi-buckling
Finally, we show a peculiar type of adhesion-induced
deformation which cannot be classiﬁed into the four
regimes as described in Section 4.1. Figure 16 showsthe equilibrium conﬁguration when the parameters are(Cs / , Cb / ) = (30, 1) corresponding to a relatively strong
adhesion regime. Interestingly, the buckling transition occurs twice in this case, i.e., a new buckling takes place inside the original buckled region. We call this phenomenon
as the multi-buckling transition which is observed when

354

The European Physical Journal E

Fig. 16. Top, side and bottom views of the equilibrated conﬁgurations of an adhering spherical shell when the sets of thescaled elastic constants are (Cs /, Cb /) = (30, 1).

the adhesion is strong enough to induce the second buckling. For much larger shells, one would expect to see
the multi-buckling which exhibits more than triple bucklings. Although we are not aware of such an unique shapein the real life, it would be very interesting to ﬁnd themulti-buckled state experimentally. Biological cells adhering onto a rigid substrate may be one of the possible systems to observe the multi-buckling because the binding
energy between the cell membranes and substrate can be
fairly large. We comment here that a cascade of buckling
was observed by compressing thin plates [73].

Fig. 17. Notation of the indentation lengths H in the disk
formation regime. R is the radius of the undeformed sphericalshell. d is the radius of the circular disk.

1
(44)
ν= ,
3
√
3
Cb =
Y h3 .
(45)
16
Combining these three relations, h and Y are given by

5 Scaling theoryh = (8Cb L20 /Cs )1/2 ,Based on the continuum shell theory, as described in
Section 2, we now interpret the deformation of the shellwithin the scaling argument [18]. Attention will be paid
to the cases in Figures 2b and c, i.e., the disk formation
regime (case I) and the isotropic buckling regime (case II).
First, we discuss how the parameters in the discretizedmodel are related to those in the continuum theory suchas the Young’s modulus or the Poisson’s ratio. We remindagain that both Cs and Cb in the discretized model havethe dimension of energy. Then according to equation (25),
we can relate them as

1 Cs
Y = √
·
6 Cb1/2 L30

Yh
Cs
∼
,
L20
1 − ν2

(39)

Y h3
,
(40)
1 − ν2
except the prefactors. Note that the three-dimensional
Young’s modulus Y has the dimension of energy per volume, and the Poisson’s ratio is dimensionless. From theabove relations, the eﬀective thickness h and the Young’s
modulus are given byCb ∼h ∼ (Cb L20 /Cs )1/2 ,1/2

Y ∼ Cs3/2 /Cb L30 .

(41)
(42)

In order to determine the numerical factors, the geometryof the network should be speciﬁed. For a two-dimensional
triangular lattice, Seung and Nelson showed that the following relations holds [63,74]:
√
Cs
3
Y h,
(43)
=
L20
2

(46)

3/2

(47)

In the disk formation regime (case I), the eﬀect of adhesion is weak so that the shell deforms only slightly atthe bottom as we have seen in Figure 2b. This situation
is analogous to the case of a shell subjected to a small localized force [18]. Let d be the dimension of the deformed
region which is caused by the contact between the shelland the substrate as depicted in Figure 17. The out-ofplane displacement ζ (see Eq. (11)) in the deformed region
can be identiﬁed as the indentation length H, i.e., ζ ∼ H.
Following the continuum treatment, the strain tensor isof the order of ζ/R ∼ H/R. Hence the total stretchingenergy is
Es ∼ Y hH 2 d2 /R2 ,(48)
where we have multiplied the area of the deformed region
d2 . The fact that ζ varies considerably over a distance dgives the curvature ζ/d2 ∼ H/d2 . Then the total bending
energy behaves as
Eb ∼ Y h3 H 2 /d2 .

(49)

Note that the stretching energy increases and bending energy decreases with increasing d.
The size d is provided by the condition that these twoenergies balance:d ∼ (hR)1/2 ∼ (Cb L20 /Cs )1/4 R1/2 ,(50)

where we have used equation (41). Hence the area of thecontact region S (I) for the case I scales as
S (I) ∼ d2 ∼ hR ∼ (Cb L20 /Cs )1/2 R.

(51)

S. Komura et al.: Buckling of spherical shells adhering onto a rigid substrate

The minimized total elastic energy Ee = Es + Eb scalesasEe(I) ∼ Y h2 H 2 /R
∼ (Cb Cs /L20 )1/2 H 2 /R.

(52)

Varying this with respect to H gives the force:f (I) ∼ (Cb Cs /L20 )1/2 H/R,

(53)

which is proportional to H. This result indicates the linear
Hooke’s law of the deformation.
In the isotropic buckling regime (case II), on the other
hand, the adhesion is strong enough for the shell to undergo the buckling as in Figure 2c. Then most of the elastic energy is concentrated over a narrow bending strip ofwidth w and radius r as deﬁned in Figure 10. The buckled
region is assumed to be a spherical cap which is a mirror
image of its original shape. We remind that the assumption of vanishing spontaneous curvature plays an important role here. Then the following relation holds according
to the simple geometrical reason [18]:r ∼ H 1/2 R1/2 ,(54)

where H is the total indentation length. This explains
the scaling relation which we found in our simulation (seeFig. 14).
Since the order of magnitude of the displacement ofa point within the bending strip is ζ ∼ wr/R, the strain
is given by ζ/R ∼ wr/R2 , and the curvature is ζ/w2 ∼r/Rw. Then the stretching energy and bending energy
scale asandEs ∼ Y h(wr/R2 )2 wr ∼ Y hw3 r3 /R4 ,(55)

Eb ∼ Y h3 (r/Rw)2 wr ∼ Y h3 r3 /R2 w,(56)

respectively. Here wr is the area of the bending strip. Minimizing these two energies with respect to w, we obtainw ∼ (hR)1/2 ∼ (Cb L20 /Cs )1/4 R1/2 .

(57)

Note that the scaling of d in equation (50) and that of w
are the same. From equations (54) and (57), the area ofthe bending strip that contacts with the substrate as in
Figure 2c is given by
S (II) ∼ wr ∼ (Cb L20 /Cs )1/4 H 1/2 R.

(58)

The minimized total elastic energy is given byEe(II) ∼ Y h5/2 r3 /R5/2∼ Y h5/2 H 3/2 /R
3/4∼ Cb (Cs /L20 )1/4 H 3/2 /R.

(59)

In this case, the required force f is
3/4

f (II) ∼ Cb (Cs /L20 )1/4 H 1/2 /R.

(60)

355

In contrast to equation (53), this relation is non-linear.
So far the discussion is valid as long as H is ﬁxed andgiven. We now consider how the indentation length H canbe related to the strength of adhesion. Let v be the van der
Waals energy per unit area. It was shown in reference [51]that v can be approximately given by v ∼ A/(12πδ 2 ),where A is the Hamaker constant and δ is an atomic cutoﬀ.
Then the total adhesion energy is estimated by
Ea ∼ vS,(61)

where S is the contact area. In the disk formation regime(case I), we use equation (51) for the contact area. If thedeformations are driven by van der Waals adhesion, theadhesion energy Ea is expected to balance with the elastic(I)energy Ee given by equation (52). By setting Ea ∼ Ee ,we arrive at the estimate for the indentation length H:
H (I) ∼ v 1/2 (L20 /Cs )1/2 R,(62)

for given v and R. In the isotropic buckling regime(case II), we use equation (58) for the contact area. By(II)setting Ea ∼ Ee , we get−1/2

H (II) ∼ vCb(L20 /Cs )1/2 R2 ,(63)

which is a diﬀerent scaling.
Comparing equations (51) and (58), we see that thecontact area of the bending strip becomes larger thanthat of the disk (S (I) < S (II) ) when the relation H >(Cb L20 /Cs )1/2 holds. Since the right hand side of this inequality scales similarly with the eﬀective thickness h (seeEqs. (41) or (46)), we see that the transition from the disk
formation regime to the isotropic buckling regime occurs
typically for a deformation H ≥ h. The increase in thecontact area between the shell and the substrate resultsin the gain in the van der Waals adhesion energy, andhence W decreases when the buckling takes place as seenin Figure 9.
Let us check if this scaling argument holds true in our
numerical simulation. As a rough estimate of the transition point, we pay attention to the critical point in
Figure 13 which separates the discontinuous and continuous buckling behaviors. We mentioned in the previous
section that the critical point appears when the combination of the elastic constants are (Cb / )c ≈ 6.1 and(Cs / )c ≈ 210 for N = 812 (see also Tab. 1). By using
equation (46), we can deduce the eﬀective thickness to be(h/R)c ≈ 0.064. On the other hand, the numerically obtained critical indentation length is (H/R)c ≈ 0.18. Bytaking the ratio between H and h, we obtain the relation
H ≈ 2.8h at the critical point. This result indeed conﬁrms
the fact that the buckling transition takes place when the
indentation length exceeds the shell thickness.
We have performed the similar analysis for other shells
which have diﬀerent sizes N . The results are summarizedagain in Table 1. The ratio H/h (the last column) isroughly 2–3 with the largest uncertainty for N = 92. Remarkably enough, this value is almost independent of N ,although the location of the critical point diﬀerers significantly between the diﬀerent sizes. It is very interesting356

The European Physical Journal E

to see that the universal property of the shell manifests
itself at the critical point of the buckling transition. As
we shall discuss in the next section, the general condition
for buckling transitions seems to hold in various systems
ranging from nanoscale to macroscale.

6 Discussion
In order to bridge between our results and real materials,
we ﬁrst give some typical numbers to the model parameters. In the case of a layered material made of carbon, thetwo-dimensional Young’s modulus and the bending rigidity are roughly 1.3 × 105 erg/cm2 and 1.6 × 10−12 erg,
respectively [75]. Assuming that the adhesion energy
is the order of thermal energy kB T , we can deduce themodel parameters as Cs / ≈ 480 and Cb / ≈ 40. According to Figure 13, the adhesion of a single-walled fullerene(with radius R/σ ≈ 0.75) should correspond to a point
well above the critical point. In this regime, the fullerene
may deform as in Figure 2b, which is consistent with theprevious prediction [51].
Other example of a spherical shell is a polyelectrolyte
multilayer capsule [42–45]. Such a material is produced
using layer-by-layer coating of dissolvable colloids andsubsequent dissolution of the core material. These capsules oﬀer the advantage that they can be preparedwith well-controlled radius and shell thickness. The threedimensional Young’s modulus of the capsule was measured
to be 500–750 MPa [42,43] or 1.5–2 GPa [44,45], but its
thickness h tends to be in the 10 nm range. Since this gives
fairly large bending rigidity of the order of Cb / ≈ 104 ,van der Waals adhesion only would not cause a considerable deformation of a microcapsule and belong to thesmall deformation regime as in Figure 2a.
However, it was shown in reference [45] that other attractive interaction such as electrostatic interaction leads
to a strong adhesion of microcapsules. In fact, anionicmicrocapsules on cationic glass resulted in a truncatedsphere topology with a circular adhesion disk. This situation obviously corresponds to the disk formation regime
in our simulation. The dependence of the adhesion diskon the shell thickness is found to be in agreement withthe previous theoretical prediction [18,51]. Furthermore,
microcapsules become unstable and buckle due to the osmotic pressure diﬀerence between inside and outside theshell [42,43]. The critical osmotic pressure depends on thecapsule radius and the shell thickness.
As brieﬂy mentioned in Introduction, there are several controlled mechanical experiments which deal withthe buckling of spherical shells. For example, Dubreuil
et al. compressed the above mentioned polyelectrolyte microcapsules using AFM [44]. The shape of the deformed
shell was monitored by reﬂection interference contrast microscopy. They measured the relation between force anddeformation, and revealed that the capsule ﬁrst deforms
only weakly. As the deformation becomes larger, an increase in the contact area is observed, which is followedby the buckling transition. In a more recent investigation on the same system, both the isotropic buckling andanisotropic buckling are distinguished [76]. Although ahigh hysteresis between the loading and unloading curve
was detected, the capsule stayed elastic. The observed sequence of deformation caused by the compression is very
similar to what we see in our simulation.
Pauchard and Rica studied the deformation of a pingpong ball which is forced to be in contact with a rigid
plate [37,38]. In their work, the boundary of the halfsphere was ﬁxed in order to avoid non-axisymmetric deformations. For low applied forces, the shell ﬂattens against
the horizontal plate. For higher compression forces, a discontinuous buckling transition occurs when the deformation is close to twice the thickness of the shell; H/h ≈ 2.4.
This value cannot be directly compared with the corresponding ratio at the critical point in our simulation(H/h ≈ 2.8) since the buckling is discontinuous for a pingpong ball. However the fact that the deformation becomeslarger than the shell thickness is the required conditioneven for the buckling of a ping-pong ball. Moreover, thesequence of the deformation is in good agreement with
our simulation results, although their experimental set-up
is not identical to our model of adhesion. Interestingly,
the polygonal structures associated with the anisotropic
buckling as in Figure 2d were also formed when a localized point force was applied to the shell [37,38].
In a smaller scale experiment, the microrheology ofself-assembled actin-coated vesicles was studied using optical tweezers and single-particle tracking [35]. The actin
ﬁlaments mimic cytoskeletal networks in cells, and they increase the bending modulus of the membrane up to around
100 kB T . A buckling instability was observed when a large
localized force of the order of 0.5 pN is applied perpendicular to this vesicle [36]. This deformation involves boththe stretching and bending contributions, and has beeninterpreted in terms of the shell theory as in the present
paper. The thickness of the actin-coated vesicle is roughlyh ≈ 100 nm, and it buckles when the deformation exceeds
H ≈ 200 nm. Hence the ratio between the two lengths is
H/h ≈ 2 at the threshold of the buckling. We therefore seethat the condition H > h determines the onset of various
types of buckling transitions in diﬀerent length scales.
As a result of the Delaunay triangulation of the spherical surface, there are always 12 grid points which have ﬁveneighbors as explained in Section 3. Recently, the facetingof spherical shells associated with 12 isolated points ofﬁvefold symmetry was argued by Lidmar et al. [70]. They
introduced the so called Föpple-von Kármán number of aspherical shell deﬁned by
τ=

Ŷ R2
,
κ

(64)

where Ŷ is the two-dimensional Young’s modulus [74].
From the relation Ŷ = Y h and equation (26), we notethat τ is proportional to the square of the ratio between
the radius R and thickness h of the shell; τ ∼ (R/h)2 . It isreported in reference [70] that, in the absence of adhesion,a signiﬁcant deviation from a perfect spherical shape takes
place when τ becomes of the order of 103 . This instability

S. Komura et al.: Buckling of spherical shells adhering onto a rigid substrate357

results in the faceting of the shell, which is manifested in
suﬃciently large viruses composed of protein capsomers.

7 Conclusion

As a rough comparison with this prediction, we estimate the Föpple-von Kármán number τ at the critical point for each size N . Using the numbers of (h/R)c
listed in Table 1, we can estimate τ to be less than 103 for
N = 92, 272, 812, but becomes roughly 103 for N = 2432.
Hence the faceting of the shell may be irrelevant for mostof the cases as long as the buckling is concerned. One difference between our model and that used in reference [70]is that not all the springs have the equal natural length
in the present case. More precisely, the springs which areconnected to the ﬁve-handed beads have slightly smaller
natural length than those connected to the six-handedbeads. Hence the spherical shape is more stable in our
numerical simulation and the defect-induced buckling is
suppressed. However, it is possible that the anisotropic
buckling is triggered by the singular disclinations as mentioned in Section 4. This can be important when the buckled region creates a polygonal structure.

We have investigated the deformation of the elastic shell
adhering onto the substrate both numerically and theoretically. The sum of the stretching, bending, and adhesion energies is minimized using the conjugate gradientmethod. The deformation of the shell is characterized by
the dimensionless parameters Cs / and Cb / . There are
four diﬀerent regimes of deformation: (i) small deformation regime, (ii) disk formation regime, (iii) isotropic buckling regime, and (iv) anisotropic buckling regime. As forthe buckling transition, there are both discontinuous andcontinuous cases for large and small Cs / , respectively.
These diﬀerent cases are separated by the critical point.
According to the scaling arguments, the buckling transition takes place when the indentation length exceeds
the eﬀective shell thickness, which is in good agreementwith our numerical results. Moreover, the ratio between
the indentation length and its thickness close to the critical point is roughly 2 even for diﬀerent shell sizes. This
general condition seems to hold in various experimental
systems ranging from nanoscale to macroscale.

In a recent experiment by Pauchard and Couder, thebuckling of shell-shaped membranes was observed using
droplets of suspension [40]. As evaporation goes on, aspherical droplet on a super-hydrophobic substrate ﬁrst
ﬂattens at the top. Then the buckling starts at the top of
the droplet, and the inverted region grows into an invagination. Although this behavior is dynamic in its nature,
the sequence of deformation is similar to what we havedescribed in the present paper. In the last stage, a transition to a toroidal shape was observed, which is interestingin the context of gastrulation of embryos. This phenomena is attributed to the inhomogeneity of the shell, i.e.,
the elastic constants of the ﬂattened part is smaller than
the rest of the shell. The inhomogeneity in the elastic constants can play an important role such as in the domain
formation in lipid bilayers [77]. Generalization to take into
account the inhomogeneous elastic constants is straightforward, and will be examined in the future.
For red blood cells, it is reported that a strong adhesion produces a ﬁnite membrane tension [78,79]. Such a
spreading-induced tension can cause the rupture of cells.
We have not included the eﬀect of tension in our simulation, but can make the following argument. Since thepresence of a positive tension tends to shrink the total
area of the shell, we expect that it will eﬀectively reduce
the natural length of each spring. Since this results in asmaller shell thickness (see Eq. (41)), it it possible that ashell with a tension can buckle easier than a tensionlessshell. In reality, buckling may induce rupture since theelastic energies will be localized in a small region.
As a ﬁnal remark, we note that the notion of spontaneous curvature of shells is diﬀerent from that of ﬂuid
membranes. For elastic shells, one needs to introduce defects to produce a preferred curvature, which depends on
how the shell is prepared. Hence the spontaneous curvature of shells is induced by a kinetic eﬀect. For ﬂuid
membranes, on the other hand, spontaneous curvature isindeed a material constant.

We thank A. Fery, R. Lipowksy, and S.A. Safran for useful discussions. This work is supported by the Ministry of Education,
Culture, Sports, Science and Technology, Japan (Grant-in-Aid
for Scientiﬁc Research No. 15540395).

References
1. T. Hwa, E. Kokufuta, T. Tanaka, Phys. Rev. A 44, R2867(1991)
2. X. Wen, C. Garland, T. Hwa, M. Kardar, E. Kokufuta, Y.
Li, M. Orkisz, T. Tanaka, Nature 355, 426 (1992)
3. M.S. Spector, E. Naranjo, S. Chiruvolu, J.A. Zasadzinski,
Phys. Rev. Lett. 73, 2867 (1994)
4. L. Bourdieu, J. Daillant, D. Chatenay, A. Braslau, D.
Colson, Phys. Rev. Lett. 72, 1502 (1994)
5. A.
Saint-Jalmes,
F.
Graner,F.
Gallet,B.
Houchmandzadeh, Europhys. Lett. 28, 565 (1994)
6. A. Saint-Jalmes, F. Gallet, Eur. Phys. J. B 2, 489 (1998)
7. C.F. Schmidt, K. Svoboda, N. Lei, I.B. Petsche, L.E.
Berman, C.R. Saﬁnya, G.S. Grest, Science 259, 952 (1993)
8. A.A. Boulbitch, Phys. Rev. E 57, 2123 (1998)
9. E.M. Kramer, T.A. Witten, Phys. Rev. Lett. 78, 1303
(1997)
10. B.A. DiDonna, T.A. Witten, S.C. Venkataramani, E.M.
Kramer, Phys. Rev. E 65, 016603 (2001)
11. T.A. Witten, H. Li, Europhys. Lett. 23, 51 (1993)
12. A.E. Lobkovsky, S. Gentges, H. Li, D. Morse, T.A. Witten,
Science 270, 1482 (1995)
13. A.E. Lobkovsky, Phys. Rev. E 53, 3750 (1996)
14. A.E. Lobkovsky, T.A. Witten, Phys. Rev. E 55, 1577
(1997)
15. B.A. DiDonna, T.A. Witten, Phys. Rev. Lett. 87, 206105(2001)
16. B.A. DiDonna, Phys. Rev. E 66, 016601 (2002)
17. K. Matan, R.B. Williams, T.A. Witten, S.R. Nagel, Phys.
Rev. Lett. 88, 076101 (2002)

358

The European Physical Journal E

18. L.D. Landau, E. M. Lifshitz, Theory of Elasticity(Pergamon, Oxford, 1986)
19. M. Ben Amar, Y. Pomeau, Proc. R. Soc. Lond. A 453,
729 (1997)
20. E. Cerda, L. Mahadevan, Phys. Rev. Lett. 80, 2358 (1998)
21. S. Chaı̈eb, F. Melo, J.-C. Géminard, Phys. Rev. Lett. 80,
2354 (1998)
22. S. Chaı̈eb, F. Melo, Phys. Rev. E 60, 6091 (1999)
23. E. Cerda, S. Chaı̈eb, F. Melo, L. Mahadevan, Nature 401,
46 (1999)
24. B. Audoly, Phys. Rev. Lett. 83, 4124 (1999)
25. B. Audoly, B. Roman, A. Pocheau, Eur. Phys. J. B 27, 7(2002)
26. B. Audoly, A. Boudaoud, Phys. Rev. Lett. 91, 086105(2003)
27. A. Boudaoud, P. Patrı́cio, Y. Couder, M. Ben Amar,
Nature 407, 718 (2000)
28. T. Mora, A. Boudaoud, Europhys. Lett. 59, 41 (2002)
29. E. Cerda, K. Ravi-Chandar, L. Mahadevan, Nature 419,
579 (2002)
30. E. Cerda, L. Mahadevan, Phys. Rev. Lett. 90, 074302(2003)
31. S. Komura, R. Lipowsky, J. Phys. France II 2, 1563 (1992)
32. Z. Zhang, H.T. Davis, D.M. Kroll, Phys. Rev. E 48, R651
(1993)
33. S. Komura, in Vesicles, edited by M. Rosof (Marcel Dekker,
1996), pp. 198–236
34. H. Yoon, J.M. Deutch, Phys. Rev. E 56, 3412 (1997)
35. E. Helfer, S. Harlepp, L. Bourdieu, J. Robert, F.C.
MacKintosh, D. Chatenay, Phys. Rev. Lett. 85, 457 (2000)
36. E. Helfer, S. Harlepp, L. Bourdieu, J. Robert, F.C.
MacKintosh, D. Chatenay, Phys. Rev. Lett. 87, 088103(2001)
37. L. Pauchard, Y. Pomeau, S. Rica, C. R. Acad. Sci. Paris
324, 411 (1997)
38. L. Pauchard, S. Rica, Phil. Mag. B 78, 225 (1998)
39. L. Pauchard, C. Allain, Europhys. Lett. 62, 897 (2003)
40. L. Pauchard, Y. Couder, Europhys. Lett. 66, 667 (2004)
41. N. Tsapis, E.R. Dufresne, S.S. Sinha, C.S. Riera, J.W.
Hutchinson, L. Mahadevan, D.A. Weitz, Phys. Rev. Lett.
94, 018302 (2005)
42. C. Gao, S. Leporatti, S. Moya, E. Donath, H. Möhwald,
Langmuir 17, 3491 (2001)
43. C. Gao, E. Donath, S. Moya, V. Dudnik, H. Möhwald, Eur.
Phys. J. E 5, 21 (2001)
44. F. Dubreuil, N. Elsner, A. Fery, Eur. Phys. J. E 12, 215(2003)
45. N. Elsner, F. Dubreuil, A. Fery, Phys. Rev. E 69, 031802(2004)
46. O.I. Vinogradova, J. Phys.: Condens. Matter 16, R1105(2004) and references therein
47. T. Hertel, R. Martel, P. Avouris, J. Phys. Chem. B 102,910 (1998)
48. T. Hertel, R.E. Walkup, P. Avouris, Phys. Rev. B 58,
13870 (1998)
49. M.-F. Yu, T. Kowalewski, R.S. Ruoﬀ, Phys. Rev. Lett. 86,
87 (2001)

50. S. Komura, K. Tamura, T. Kato, Eur. Phys. J. E 13, 73(2004)
51. U.S. Schwarz, S. Komura, S.A. Safran, Europhys. Lett. 50,
762 (2000)
52. K. Miura, S. Kamiya, N. Sasaki, Phys. Rev. Lett. 90,
055509 (2003)
53. K. Tamura, S. Komura, T. Kato, J. Phys.: Condens.
Matter 16, L421 (2004)
54. H. Hertz, J. Reine Angew. Math. 92, 156 (1881)
55. K.L. Johnson, K. Kendall, A.D. Roberts, Proc. R. Soc.
Lond. A 324, 302 (1971)
56. J.N. Israelachivili, Intermolecular and Surface Forces(Academic, London, 1991)
57. F.I. Niordson, Shell Theory (North Holland, New York,1985)
58. The factor 12 is conventional. See also Ref. [57]59. M. Tanemura, T. Ogawa, N. Ogita, J. Comp. Phys. 51,
191 (1983)
60. J. M. Augenbaum, C.S. Peskin, J. Comp. Phys. 59, 171(1985)
61. A. Baumgärtner, J.-S. Ho, Phys. Rev. A 41, 5747 (1989)
62. S. Komura, A. Baumgärtner, Phys. Rev. A 44, 3511 (1991)
63. H.S. Seung, D.R. Nelson, Phys. Rev. A 38, 1005 (1988)
64. Y. Kantor, D.R. Nelson, Phys. Rev. Lett. 58, 2774 (1987)
65. Y. Kantor, D.R. Nelson, Phys. Rev. A 36, 4024 (1987)
66. S.A. Safran, Statistical Thermodynamics of Surface,
Interfaces, and Membranes (Addison-Wesley, New York,1994)
67. W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T.
Vetterling, Numerical Recipes (Cambridge University,
New York, 1989)
68. G. Gompper, D.M. Kroll, Phys. Rev. A 46, 7466 (1992)
69. A.R. Bausch, M.J. Bowick, A. Cacciuto, A.D. Dinsmore,
M.F. Hsu, D.R. Nelson, M.G. Nikolaides, A. Travesset,
D.A. Weitz, Science 299, 1716 (2003)
70. J. Lidmar, L. Mirny, D.R. Nelson, Phys. Rev. E 68, 051910(2003)
71. Y. Kantor, M. Kardar, D.R. Nelson, Phys. Rev. Lett. 57,
791 (1986)
72. Y. Kantor, M. Kardar, D.R. Nelson, Phys. Rev. A 35, 3056
(1987)
73. B. Roman, A. Pocheau, Europhys. Lett. 46, 602 (1999)
74. Notice that the two-dimensional Young’s modulus Ŷ(which is used in Ref. [63]) and the three-dimensional
Young’s modulus Y are related by Ŷ = Y h. Poisson’s ratiois the same in the both dimensions
75. J. Tersoﬀ, Phys. Rev. B 46, 15546 (1992)
76. A. Fery (private communication)
77. T. Baumgart, S.T. Hess, W.W. Webb, Nature 425, 821(2003)
78. A. Hategan, R. Law, S. Kahn, D.E. Discher, Biophys. J.
85, 2746 (2003)
79. A. Hategan, K. Sengupta, S. Kahn, E. Sackmann, D.E.
Discher, Biophys. J. 87, 3547 (2004)