Phase structure of a surface model on dynamically triangulated spheres
  with elastic skeletons
Phase structure of a surface model on dynamically triangulated spheres with elastic
skeletons
Hiroshi Koibuchi∗
Department of Mechanical and Systems Engineering,
Ibaraki National College of Technology, Nakane 866 Hitachinaka, Ibaraki 312-8508, Japan
(Dated: August 12, 2021)
We find three distinct phases; a tubular phase, a planar phase, and the spherical phase, in a
triangulated fluid surface model. It is also found that these phases are separated by discontinu-
ous transitions. The fluid surface model is investigated within the framework of the conventional
curvature model by using the canonical Monte Carlo simulations with dynamical triangulations.
The mechanical strength of the surface is given only by skeletons, and no two-dimensional bending
energy is assumed in the Hamiltonian. The skeletons are composed of elastic linear-chains and rigid
junctions and form a compartmentalized structure on the surface, and for this reason the vertices
of triangles can diffuse freely only inside the compartments. As a consequence, an inhomogeneous
structure is introduced in the model; the surface strength inside the compartments is different from
the surface strength on the compartments. However, the rotational symmetry is not influenced by
the elastic skeletons; there is no specific direction on the surface. In addition to the three phases
mentioned above, a collapsed phase is expected to exist in the low bending rigidity regime that was
not studied here. The inhomogeneous structure and the fluidity of vertices are considered to be the
origin of such variety of phases.
PACS numbers: 64.60.-i, 68.60.-p, 87.16.Dg
I. INTRODUCTION
A crumpling of surfaces has been investigated on the
basis of the singularity analysis, and progress has been
recently made on understanding the crumpling phenom-
ena; the universal structure on the crumpled thin sheets
was found in the formations of singularity of ridges and
cones [1, 2]. A similar transition to this phenomena was
also found experimentally between the smooth state and
the crumpled state in an artificial membrane, which is
partly polymerized [3].
Studies have also been focused on the transition in the
surface model of Helfrich, Polyakov and Kleinert (HPK)
[4, 5, 6] from the viewpoint of statistical mechanics [7,
8, 9, 10, 11, 12, 13]. The bending rigidity is known to
be stiffened by the thermal fluctuation of the surface,
and this was confirmed in the statistical mechanics of
membranes [14, 15, 16, 17, 18]. Numerical studies were
made to understand the transition in triangulated surface
models [19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. The
transition was reported as first-order in recent numerical
studies [31, 32].
On the other hand, the concern with inhomogeneous
surfaces has been growing over the past decade [33, 34]. A
homogeneous artificial membrane that is coated by elas-
tic skeletons is also considered to be an inhomogeneous
membrane. Some of the mechanical properties of such
membranes were revealed experimentally [35]. The hop
diffusion of membrane protein or lipids was observed, and
as a consequence the compartment of cytoskeletons was
∗Electronic address: koibuchi@mech.ibaraki-ct.ac.jp
confirmed to be in biological membranes [36]. It is also
well known that the microtubule, which is an element
of the cytoskeleton, gives a mechanical strength to the
surface of the biological membranes.
However, the surface collapsing phenomena and the
surface fluctuation phenomena are almost unknown in
such inhomogeneous models for membranes. Therefore,
it is worthwhile to study an inhomogeneous fluid sur-
face model within the framework of the conventional sur-
face model of HPK. We note that the inhomogeneity in
our model corresponds to the cytoskeletons in biologi-
cal membranes as stated above. The fluidity realized by
dynamical triangulations in the inhomogeneous model,
as well as the fluidity in the homogeneous surface mod-
els, corresponds to the lateral diffusion of lipids in mem-
branes.
In this paper we study a compartmentalized surface
model by Monte Carlo (MC) simulations. The Hamil-
tonian of the model includes no two-dimensional bend-
ing energy but a one-dimensional bending energy. The
model is defined on dynamically triangulated surfaces,
where the free diffusion of vertices is confined inside the
compartments. The mechanical strength of the surface is
given only by the compartment boundary, which is com-
posed of one-dimensional elastic chains and rigid junc-
tions. Because the collapsed phase is expected to appear
at sufficiently small bending rigidity b[kT ]→0(b 6=0), we
concentrate on the phase structure at relatively large b
in this paper. Consequently, information on the phase
boundary at b→0 remains unanswered.
We recently reported numerical results of three types
of surface models [37, 38], which are similar to the model
in this paper. Then, we should comment on the similar-
ity/difference between the model in this paper and the
http://arxiv.org/abs/0704.0493v1
mailto:koibuchi@mech.ibaraki-ct.ac.jp
models in [37, 38]. Firstly, the lattice structure of the
model in this paper is very similar to that of the first
model in [37] and that of the model in [38], and is iden-
tical to that of the second model in [37]. Secondly, the
lattice in this paper and that of the first model in [37]
are the dynamically triangulated one, while the lattice of
the second model in [37] and that in [38] are the fixed-
connectivity one. Thirdly, the Hamiltonian is different
from the one in the first model in [37]. The Hamiltonian
of the model of this paper includes only one-dimensional
bending energy, which is defined on the compartment
boundary, while the Hamiltonian of the first in [37] in-
cludes only a two-dimensional bending energy, which
is defined all over the surface, and no one-dimensional
bending energy is given to the compartment boundary.
Therefore, the model in this paper is different from the
three models in [37, 38].
Our results obtained in this paper show that the model
undergoes a first-order transition between the smooth
phase and the crumpled phase. Moreover, the smooth
phase can be divided into the spherical phase and the
planar phase, and the crumpled phase can also be divided
into the tubular phase and the collapsed phase, which is
expected to appear at sufficiently small b because no self-
avoiding property [39, 40, 41] is assumed in the model.
It must be emphasized that such variety of phases can
be seen neither in the conventional surface models nor in
the compartmentalized models such as those in [37, 38].
One remarkable result is the appearance of planar sur-
faces. The echinocytic shapes of erythrocytes were ex-
tensively studied, and they are currently known to be
described by many models such as the area difference
bilayer model [34]. The shape of membranes is also sen-
sitive to the flow fields [42]. Our model in this paper indi-
cates that one possible origin of such planar shape comes
from the inhomogeneity due to the cytoskeltal structure
and the fluidity of lateral diffusion of vertices.
II. MODEL
Figure 1(a) shows a triangulated surface of size
(N,NS , NJ , L) = (2322, 600, 42, 6), where N is the total
number of vertices including the junctions, NS is the total
number of vertices on the chains, NJ is the total number
of junctions, and L is the length of chains between the
two nearest-neighbor junctions. It should be noted again
that NJ is included in N ; junctions are counted in the
total number of vertices. The junctions are assumed as
rigid plates; twelve of them are pentagon and all the oth-
ers are hexagon. The junction size in Fig.1(a) is drawn
many times larger than that of the lattices for the sim-
ulations; and it will be discussed in the last part of this
section. Thick lines on the surface in Fig.1(a) denote the
chains, which are terminated at the junctions.
The construction of the lattices is as follows: Let us
start with the icosahedron. Every edge of the icosahe-
dron is divided into ℓ pieces of uniform length, and then
(a) (N,NS , NJ , L)=
(2322, 600, 42, 6)
(b) A rigid junction
with the chains
θ(ij)θ(ij)
FIG. 1: (Color online) (a) Starting configuration of surfaces
of size (N,NS, NJ , L)=(2322, 600, 42, 6), and (b) angles θ(ij)
in the bending energy S2 of Eq.(4). Thick lines in (a) denote
the compartment boundary composed of the linear chains and
the rigid junctions of the hexagonal and the pentagonal plates,
whose size is drawn many times larger than that of the lattices
for the simulations.
we have a triangulated surface of size N0 = 10ℓ
2+2 (=
the total number of vertices on the surface). The com-
partmentalized structures are constructed by dividing ℓ
further into m pieces (m=1, 2, · · · ). Thus, we have the
chains of uniform length L= (ℓ/m)−2 when m divides
ℓ. The reason for the subtraction −2 is because of the
junctions at the two end points of the chain. Because
the compartmentalized structure is a sublattice, the to-
tal number of junctions NJ is given by NJ = 10m
The total number of bonds in the sublattice is 3NJ−6,
and each bond contains L−1 vertices, then NS is given
by NS =(3NJ−6)(L−1), which can be written as NS =
30m(ℓ−3m). The hexagonal (pentagonal) rigid junctions
are composed of 7 (6) vertices, then NJ−12 hexagonal
rigid junctions and 12 pentagonal rigid junctions reduce
the total number of vertices N0 by (NJ−12)×6 and 12×5.
Therefore, we have N=N0−6NJ+12, which can also be
written as N=10ℓ2−60m2+2. The thermodynamic limit
of our model is defined by N→∞, NS→∞, and NJ →∞
under the condition that L is finite. We have the thermo-
dynamic limit at ℓ→∞ and m→∞. The lattice of size
(N,NS , NJ , L)=(2322, 600, 42, 6) in Fig.1(a) is given by
two independent integers (ℓ,m)=(16, 2).
The surfaces can be characterized by the length L. In
this paper, we assume three values for L such that
L = 6, L = 8, L = 11. (1)
The value of L has a one to one correspondence with the
total number of vertices n in a compartment; in fact, the
values of L in Eq.(1) correspond to n= 21, n= 36, and
n=66, respectively [37]. We note that the effective phys-
ical meaning of increasing (decreasing) L can be consid-
ered as the increasing (decreasing) temperature. In fact,
the surface fluctuation mainly comes from the thermal
fluctuation of vertices inside the compartments. Because
no bending energy is assumed inside the compartments,
the fluctuation of vertices becomes large not only in the
in-plane directions (free diffusion) but also in the direc-
tion perpendicular to the surface. Thus, we consider that
the fluctuations are expected to grow with increasing n,
i.e., increasing L.
We use the surfaces of size (N,NS , NJ) listed in Ta-
ble I. Three different sizes (N,NS , NJ) are assumed for
each L. The corresponding integers (ℓ,m) are as follows:
(16, 2), (24, 3), and (32, 4) for the L=6 surfaces, (10, 1),
(20, 2), and (30, 3) for the L = 8 surfaces, and (13, 1),
(26, 2), and (39, 3) for the L=11 surfaces.
TABLE I: The surface size assumed in the simulations. Three
sizes (N,NS , NJ ) are assumed for each L.
L (N,NS , NJ ) (N,NS , NJ ) (N,NS , NJ )
6 (2322,600,42) (5222,1350,92) (9282,2400,162)
8 (942,210,12) (3762,840,42) (8462,1890,92)
11 (1632,300,12) (6522,1200,42) (14672,2700,92)
The model is defined by the partition function
dXi exp [−S(X, T )] , (2)
S(X, T ) = S1 + bS2,
where S1 is the Gaussian bond potential, which is de-
fined all over the surface, and S2 is the one-dimensional
bending energy, which is defined on the compartment
boundary and will be given below. The parameter b is
the bending rigidity. The integration symbol
in Eq.(2)
denotes that the center of mass of the surface is fixed.
denotes the sum over all possible triangulations T ,
which are performed by the bond flip technique keeping
the compartments unflipped. The bond flip procedure
will be given in the following section.
The integration measure
i=1 dXi is given by the
product
dXi =
qαj(i)
(α = 3/2, 0), (3)
where N ′ (=N−NJ) is the total number of vertices ex-
cluding the junctions,
i=1 dXiq
i denotes the integra-
tion over the 3D translational degrees of freedom (DOF)
of the vertices i, and
i=1 dXi
j(i) q
j(i) denotes those
of the 3D translational DOF and the 3D rotational DOF
of the junctions i. The co-ordination number qi is the
total number of bonds meeting at the vertex i, and qj(i)
is the total number of bonds meeting at the corner j(i)
of the junction i.
The parameter α was chosen to be α=3/2 in [43, 44],
while α = 0 in many previous simulations on dynami-
cally triangulated surfaces in the literatures. It is easy to
understand that large positive α suppresses the configu-
rations with large coordination number. Therefore, it is
interesting to see the dependence of the phase structure
on α.
We chose both α = 3/2 and α = 0 for the weight
qαi [43, 44], and see whether the phase structure of the
model depends on α or not. If the parameter is chosen
to α=3/2, then the coordination number qi serves as a
weight of the integration dXi, while α=0 gives the uni-
form weight. The weight
i=1 q
i can also be written
i=1 q
i = exp(α
log qi), and therefore,
i=1 q
is considered to be the co-ordination dependent term
log qi in the Hamiltonian; −α
log qi changes
its value only on dynamically triangulated surfaces.
The Gaussian term S1 and the bending energy term
S2 are defined by
(Xi −Xj)
, S2 =
1− cos θ(ij)
, (4)
where
(ij) in S1 is the sum over bonds (ij) connecting
the vertices i and j, and
(ij) in S2 is also the sum over
bonds (ij). θ(ij) in S2 is the angle between the bonds
i and j, which include virtual bonds. The virtual bonds
denote the lines between the center and the corners of the
junction; the hexagonal (pentagonal) junction contains
six (five) virtual bonds.
Figure 1(b) is a junction and the chains linked to the
junction on a fluctuating surface. Triangles are elimi-
nated from the figure. One θ(ij) shown at a corner of
the junction is defined by using a virtual bond and a real
bond in a chain, and the other θ(ij) shown at a vertex is
defined by real bonds on the same chain.
The size of the junctions can be characterized by the
edge length R, which is fixed to
R = 0.1 (edge length of the junctions). (5)
The value R = 0.1 is relatively smaller than the mean
bond length 0.707, which corresponds to the relation
S1/N = 1.5 satisfied in the equilibrium configuration of
surfaces without the rigid junctions. As we will see later,
the relation S1/N =1.5 is slightly violated in the model
of this paper because of the rigid junctions.
Here we comment on the unit of physical quantities.
Let a be the length scale of the model, then the unit
of physical quantity that has the length unit can be ex-
pressed by a; the unit of S1 is [a
2]. The surface tension
coefficient λ in λS1+bS2 has the unit [kT/a
2] and as-
sumed to be λ=1[kT/a2], and the bending rigidity b has
the unit of [kT ] as described above.
Note that the bending rigidity b in the Hamiltonian is
a microscopic quantity from the view point of statistical
mechanical model, and therefore b is not always identical
to the macroscopic bending rigidity of real physical mem-
branes. However, the microscopic value b of real mem-
branes can effectively be varied with the temperature,
because b has the unit of kT . Therefore, it is possible to
consider that the phase structure described in terms of b
in the surface model corresponds to the phase structure
described in terms of T in real physical membranes. The
length scale a in the model is also a microscopic quantity
and, we consider that a is sufficiently smaller than the
membrane size.
III. MONTE CARLO TECHNIQUE
A sequence of random numbers called Mersenne
Twister [45] is used in the canonical MC simulations.
The Metropolis technique is applied to update X and T ,
where the variable X denotes the position of the ver-
tices and that of the junctions. The vertex position
X is shifted so that X ′ = X + δX , where δX is ran-
domly chosen in a small sphere. The new position X ′ is
accepted with the probability Min[1, exp(−∆S)], where
∆S = S(new)−S(old). The position X of a hexagonal
(or pentagonal) junction, which is not a point but a rigid
plate, is also integrated out by performing 3D random
translations and 3D random rotations.
Thus, the variable X is updated by a random N ′ (=
N−NJ) shifts of vertices, a random NJ translations of
junctions, and a randomNJ rotations of junctions. These
updates are denoted by (N ′, NJ , NJ) updates of X . The
N ′ shifts of X can be divided into NS shift of the vertices
on the linear chains and N ′−NS shifts of all the other
vertices, which are those inside the compartments.
The radius of the small sphere for δX is fixed at the
beginning of the MC simulations in order to maintain
about 50% acceptance rate. The vertices on the linear
chains carry the bending energy S2 in Eq.(4), while all the
other vertices inside the compartments does not. There-
fore, the acceptance rate is independently controlled in
the two-groups of vertices. The radius for the random
translation of the junctions and that for the random rota-
tion are also independently chosen so that the acceptance
rates are both about 50%.
The summation over T in Z of Eq.(2) is performed by
using the standard bond flip technique [22, 23]. The flip
is accepted with the probability Min[1, exp(−∆S)]. The
acceptance rate for the bond flip is not under control and
is about 75%, which is almost independent of b.
The bonds are labeled with sequential numbers. The
total number of bonds is denoted by N ′B, which excludes
the number of bonds on the linear chains because the
bonds on the linear chains remain unflipped.
The bond flip is performed as follows: Firstly, the odd-
numbered bonds are sequentially chosen to be flipped for
the N ′B/2 updates of T , and after that, the (N,NJ , NJ)
updates of X are performed. Secondly, the remaining
even-numbered bonds are chosen to be flipped for the
N ′B/2 updates of T , and after that, the (N,NJ , NJ) up-
dates ofX are performed. Thus, the (N,NJ , NJ) updates
of X and the N ′B/2 updates of T are consecutively per-
formed, and these make one MCS (Monte Carlo Sweep).
We introduce the lower bound 1× 10−8 to the area of
triangles. No lower bound is imposed on the bond length.
IV. RESULTS OF SIMULATION
A. α = 3/2
As mentioned in Section II, we assume the value of
α in Eq.(3) as α = 3/2 and α = 0. In this subsection,
we present the results obtained under α= 3/2 by using
snapshots and figures, and in the next subsection we will
show some of the results under α=0.
The thermalization MCS is 1× 107 in almost all cases.
However, more than 1 × 108 thermalization MCS were
done close to the transition point in such cases that the
surface is trapped in one phase at first and then changes
its phase to a more stable one under a given condition.
The total number of MCS for the production of samples
is 0.8 × 108 ∼ 1.3 × 108. At the transition point, about
2 × 108 MCS was performed after the thermalization in
some cases.
(a) b= 21.2 (b) b= 21.8 (c) b= 22
(d) The section (e) The section (f) The section
FIG. 2: (Color online) The snapshots of surfaces of size
(N,NS , NJ , L) = (8462, 1890, 92, 8) obtained at (a) b = 21.2
(tubular phase), (b) b = 21.8 (planar phase), and (c) b = 22
(spherical phase), and (d),(e),(f) are the surface sections of
(a),(b),(c), respectively. α=3/2.
We show snapshots of the (N,NS , NJ , L) =
(8462, 1890, 92, 8) surface in Figs.2(a)–2(c). They were
obtained at (a) b = 21.2, (b) b = 21.8, and (c) b = 22,
which respectively corresponds to the tubular phase, the
planar phase, and the spherical phase. The snapshot of
Fig.2(b) at b=21.8 was the final configuration produced
after 2× 108 MCS including 1× 108 thermalizaion MCS;
the planar surface was stable after the thermalization
MCS. The surface sections are shown in Figs.2(d)–2(f);
the sections in Figs.2(d) and 2(e) were obtained by slicing
the surfaces perpendicular to the vertical axis, and the
section in Figs.2(f) was obtained by slicing the surface
perpendicular to a horizontal axis. All of the snapshots
were drawn in the same scale. The axis of the tubular
surface Fig.2(a) as well as the axis perpendicular to the
planar surface Fig.2(b) is spontaneously chosen.
The planar phase is stable only on the L=8 surfaces,
while it seems unstable on the L=6 surfaces and on the
L=11 surfaces. Even if the planar phase once appears on
the surfaces of L=6 and L=11 of size at least N≤9282
and N ≤ 14672, respectively, it eventually collapses into
the tubular phase. Therefore, we find that no planar
phase can be seen on the L=6 and the L=11 surfaces;
the tubular phase and the spherical phase are connected
by a discontinuous transition on those surfaces. Thus,
we understand that the planar phase appears depending
on the size of the compartments. We should note that
the planar surface may bend and fluctuate in the limit of
N→∞, and the tubular surface may also bend and wind
in the same limit.
20 21 22 23
S2/NS'
N=8462
α=1.5
tubular
planar
spherical
9 9.5 10
S2/NS'
N=9282
α=1.5
tubular spherical
46 48 50 52
spherical
S2/NS'
N=6522
α=1.5
tubular
FIG. 3: The one-dimensional bending energy S2/N
S against
b obtained on the surfaces of (a) L = 6, (b) L = 8, and (c)
L=11. N ′S(=NS+6NJ −12) is the total number of vertices
where S2 is defined.
Figures 3(a),3(b), and 3(c) show the bending energy
S of Eq.(4) against b, which were obtained on the
surfaces of L = 6, L = 8, and L = 11, respectively.
N ′S(=NS+6NJ−12) is the total number of vertices where
S2 is defined. 6NJ −12 is the total number of corners
of the junctions, which include 12-pentagons. The solid
lines on the data were drawn to guide the eyes. Dashed
lines drawn vertically denote the phase boundary be-
tween the tubular and the spherical phases, the bound-
ary between the tubular and the planar phases, and the
boundary between the planar and the spherical phases.
The discontinuous change of S2/N
S between the tubu-
lar phase and the spherical (or the planar) phase is very
clear in the figures and considered to be a sign of the
first-order transition.
In order to see the difference between S2/N
S in those
three phases, we plot in Figs.4(a),4(b), and 4(c) the varia-
tion of S2/N
S against MCS obtained at b=21.2, b=21.4,
and b = 21.8 on the (N,NS , NJ , L) = (8462, 1890, 92, 8)
surface. The thermalization MCS were not discarded;
they were included only in those variations. S2/N
b=21.2 in Fig.4(a) shows a jump from the spherical phase
to the planar phase and a jump from the planar phase to
the tubular phase; the corresponding MCS at the jumps
were indicated with the dashed vertical lines. We also
find in Fig.4(b) a jump from the spherical phase to the
planar phase. A jump is also seen in S2/N
S at b=21.8
0 0.5 1
planar
N=8462
L=8b=21.4
spherical
[x108]0 0.5 1
N=8462
b=21.2
tubularplanarspherical
[x108] 0 1 2
N=8462
b=21.8
spherical
planar
0.05 0.06
b=21.2
(d) 0.05 0.06
b=21.4
(e) 0.05 0.06
b=21.8
FIG. 4: The variation of S2/N
S against MCS, which were
obtained on the (N,NS , NJ , L)=(8462, 1890, 92, 8) surface at
(a) b=21.2, (b) b=21.4, and (c) b=21.8. The dashed lines
denote the MCS where the jumps occurred. The correspond-
ing normalized histogram h(S2) obtained at (d) b=21.2, (e)
b = 21.4, and (f) b = 21.8. The parameter α was fixed to
α=3/2.
in Fig.4(c) from the spherical phase to the planar phase.
The value of b=21.2 corresponds to the tubular phase,
whereas b = 21.4 and b= 21.8 correspond to the planar
phase, because the final states are considered to be stable
states. The surfaces at b=21.2 and b=21.8 can be seen
in the snapshots in Figs.3(a) and 3(b).
The distribution of S2/N
S are shown as the normal-
ized histograms h(S2) in Figs.4(d)–4(f), which respec-
tively correspond to the variations in Figs.4(a)–4(c). We
see that h(S2) in Fig.4(d) has three peaks; two of them
are almost overlapping and the other one is distinctly
separated from the previous two. Those three peaks in
h(S2) correspond to the spherical phase, planar phase,
and the tubular phase. Two almost overlapping peaks
can also be seen in h(S2) in Figs.4(e) and 4(f), and they
are corresponding to the spherical phase and the planar
phase. We remark that the surfaces hardly change not
only from the tubular phase to the smooth (= spherical
or planar) phase but also from the planar phase to the
spherical phase on the L = 8 and L = 11 surfaces. For
this reason, we find in Figs.4(a)–4(c) no jump-back from
a higher S2 state (such as the tubular state) to a lower
S2 state (such as the planar state).
The two-dimensional bending energy is defined by
(1− ni · nj) , (6)
where ni is the unit normal vector of the triangle i, and
ni ·nj is defined on the common bond (ij) of the triangles
i and j. S3 is not included in the Hamiltonian and is
defined even on the edges of the rigid junctions. Figures
5(a)–5(c) show S3/NB against b obtained on the surfaces
of L=6, L=8, and L=11, where NB is the total number
of bonds including the edges of the junctions. The jump
of S3/NB in Fig.5(b) is clearly seen between the tubular
20 21 22 23
S3/NB
N=8462
α=1.5
tubular
planar
spherical
9 9.5 10
S3/NB
N=9282
α=1.5
tubular spherical
46 48 50 52
spherical
S3/NB
N=6522
α=1.5
tubular
FIG. 5: The two-dimensional bending energy S3/NB against
b obtained on the surfaces of (a) L = 6, (b) L = 8, and (c)
L=11. NB is the total number of bonds where S3 is defined.
phase and the planar phase. On the contrary, S3/NB
in the planar phase in Fig.5(b), as well as S2/N
S in the
planar phase in Fig.3(b), is not so clearly distinguishable
from that in the spherical phase.
20 21 22 23
1.512
1.514
1.516
N=8462
α=1.5
tubular
planar
spherical
9 9.5 10
1.514
1.516
1.518
N=9282
α=1.5
sphericaltubular
46 48 50 52
1.508
1.512
N=6522
α=1.5
tubular spherical
FIG. 6: The Gaussian bond potential S1/N against b obtained
on the surfaces of (a) L=6, (b) L=8, and (c) L=11.
It is expected that the Gaussian bond potential S1/N
is influenced by the phase transitions. The potential
S1/N should be S1/N ≃ 3/2, which is satisfied in the
model without the rigid junctions because of the scale
invariant property of the partition function in that case.
However, the junction size R in Eq.(5) is finite in the
model of this paper, and therefore S1/N can slightly de-
viate from 3/2.
Figures 6(a)–6(c) show S1/N against b obtained on the
surfaces of (a) L=6, (b) L=8, and (c) L=11. Discontin-
uous changes in S1/N shown in the figures are consistent
with the discontinuous transitions of the model, although
the changes are very small compared to the value of S1/N
itself. We find also the expected deviation of S1/N from
3/2 in the figures.
Figures 7(a)–7(c) show the mean square sizeX2, which
is defined by
Xi − X̄
, X̄ =
Xi, (7)
where X̄ is the center of mass of the surface. We see
that the phase transition is not reflected in X2 on the
L=6 surfaces in Fig.7(a), and the transition is also not
reflected in X2 on the L= 8 surfaces in Fig.7(b) at the
transition point between the planar phase and the spher-
ical phase. To the contrary, X2 discontinuously changes
20 21 22 23
N=942
α=1.5
N=8462
N=3762
tubular spherical
planar
9 9.5 10
α=1.5
N=2322
N=9282
N=5222
tubular spherical
46 48 50 52
N=6522
α=1.5
N=1632
N=14672
tubular spherical
FIG. 7: The mean square size X2 against b obtained on the
surfaces of (a) L=6, (b) L=8, and (c) L=11.
in Fig.7(b) at the transition point between the tubular
phase and the planar phase and also at the transition
point in Fig.7(c). All of these behaviors ofX2 at the tran-
sition points are consistent with those of S2/N
S , S3/NB,
and S1/N .
B. α = 0
In this subsection, we present some of the results ob-
tained under α=0.
(a) b= 20.9 (b) b= 21.4 (c) b= 21.8
(d) The section (e) The section (f) The section
FIG. 8: (Color online) The snapshots of surfaces of size
(N,NS , NJ , L) = (8462, 1890, 92, 8) obtained at (a) b = 20.9
(tubular phase), (b) b=21.4 (planar phase), and (c) b=21.8
(spherical phase), and (d),(e),(f) are the surface sections of
(a),(b),(c), respectively. α=0.
Snapshots of surfaces of α=0 are shown in Figs.8(a),
8(b), 8(c), which respectively correspond to the tubu-
lar phase (b = 20.9), the planar phase (b = 21.4), and
the spherical phase (b = 21.8). The surface size is
(N,NS , NJ , L)= (8462, 1890, 92, 8), which is identical to
that in Fig.2. The snapshot in Fig.8(b) at b = 21.4 is
the final configuration produced after 1.9× 108 MCS in-
cluding 1 × 107 thermalizaion MCS; the planar surface
was stable throughout the simulation. Thus, we find
that three distinct phases are seen also in the surfaces
of L = 8, and that the planar phase is unstable on the
surfaces of L= 6 and L= 11 under the condition α= 0.
Therefore, we consider that the phase structure of the
model is independent of whether α=3/2 or α=0.
20 21 22 23
S2/NS'
N=8462
tubular
planar
spherical
9 9.5 10
S2/NS'
N=9282
tubular spherical
46 48 50 52
spherical
S2/NS'
N=14672
tubular
FIG. 9: The one-dimensional bending energy S2/N
S against b
obtained on the surface of (a) L=6, (b) L=8, and (c) L=11.
N ′S(=NS+6NJ−12) is the total number of vertices where S2
is defined.
The one-dimensional bending energy S2/N
S obtained
under α=0 is shown in Figs.9(a)–9(c). A discontinuous
change can be seen in S2/N
S not only in Fig.9(b) at the
phase boundary between the tubular phase and the pla-
nar phase but also in Fig.9(c) at the phase boundary be-
tween the tubular phase and the spherical phase. A jump
of S2/N
S in Fig.9(b) at the transition point between the
planar phase and the spherical phase is very small, and
hence is hardly seen just the same as in Fig.3(b) under
α = 3/2 in the previous subsection. Thus, we find no
difference between S2/N
S of α=0 and that of α=3/2.
20 21 22 23
N=942
N=8462
N=3762
tubular spherical
planar
9 9.5 10
N=2322
N=9282
N=5222
tubular spherical
46 48 50 52
N=14672
tubular spherical
FIG. 10: The mean square size X2 against b obtained on the
surfaces of (a) L=6, (b) L=8, and (c) L=11.
The mean square size X2 are shown in Figs.10(a)–
10(c). A jump is also seen in X2 on the L = 8 and
L=11 surfaces in Figs.10(b) and 10(c), and it is hardly
seen on the L=6 surfaces of size up to (N,NS , NJ , L)=
(9282, 2400, 162, 6). These results are identical to those
observed in Figs.7(a)–7(c) under α=3/2.
Finally, we comment on the planar phase appeared
only on the L = 8 surface. The thermal fluctuation of
vertices inside the compartments disorders the surface
against the bending energy of the compartment bound-
ary. Therefore, the strength to disorder the surface in-
creases (decreases) with increasing (decreasing) L if N
remains fixed, as stated in Section II. On the other
hand, the mechanical strength of the surface increases
(decreases) with decreasing (increasing) L, because the
total number of junctions increases (decreases) with de-
creasing (increasing) L. Therefore, the strength to or-
der the surface increases (decreases) with decreasing (in-
creasing) L. Then, we expect that the surface is ordered
(disordered) at sufficiently small (large) L at given in-
termediate value of b. Moreover, it seems possible that
two competitive forces to order/disorder the surface are
balanced with each other at intermediate values of L and
consequently, some new phase appears depending on b
at those L. Note also that the possibility of the appear-
ance of planar phase is not completely eliminated on the
surfaces of L=6 and L=11 of sufficiently large size.
V. SUMMARY AND CONCLUSION
We have shown that a dynamically triangulated spher-
ical surface has three distinct phases; the tubular phase,
the planar phase, and the spherical phase, and that they
are separated by discontinuous transitions. The first-
order nature was very clear from the discontinuity in the
bending energies S2 and S3 not only at the transition
point between the tubular phase and the planar phase
but also at the transition point between the tubular phase
and the spherical phase. We know that the model has
the collapsed phase at sufficiently small b, since the self-
avoiding property is not assumed at least. Therefore, we
expect that the model has four different phases including
the collapsed phase, although the order of the transition
between the collapsed phase and the tubular phase is un-
known.
The mechanical strength of the surface is given only
by elastic linear-chains with rigid junctions. The
triangulated surfaces are characterized by the size
(N,NS , NJ , L), where N is the total number of vertices
including the junctions, NS is the total number of vertices
on the chains, NJ is the total number of junctions, and L
is the length of chains between the two nearest-neighbor
junctions on the starting configurations. These four pa-
rameters are not totally independent, because these are
given by two independent integers (ℓ,m), where m di-
vides ℓ. In fact, N =10ℓ2−60m2+2, NS =30m(ℓ−3m),
NJ =10m
2+2, and L=(ℓ/m)−2.
We assumed three different values for L such that L=
6, L=8, and L=11 in the simulations. The edge length
R of the rigid junction was fixed to be R = 0.1. The
parameter α, which represents a weight for the three-
dimensional integrations of the partition function, was
assumed as α=3/2 and α=0.
It is remarkable that the model has the planar phase,
which is stable only on the surfaces with a specific struc-
ture. In fact, the planar phase can be seen on the surfaces
of L=8, and it is unstable on the L=6 and L=11 sur-
faces. The planar phase appears in a narrow region on
the b-axis between the tubular phase and the spherical
phase, and it is distinguishable from the spherical phase
because a small but finite discontinuity can be seen in
the bending energies S2/N
S and S3/NB. The gap of the
bending energy S2 at the transition point is very small,
i.e., S2 in the planar phase is almost identical to that
in the spherical phase; however, the double peak struc-
ture was clearly seen in the histogram of S2, which is
included in the Hamiltonian. From this, we confirmed
that the transition between the planar phase and the
spherical phase is of first order. Our model in this pa-
per indicates that one possible origin of planar shape of
spherical membranes comes from the inhomogeneity due
to the cytoskeltal structure and the fluidity of lateral dif-
fusion of vertices.
We have confirmed that the results obtained at α=3/2
in Eq.(3) remain unchanged when α = 0. The phase
structure of the fluid surface model in this paper is in-
dependent of the choice of α at least for α = 3/2 and
α= 0. Large scale simulations should be performed. It
remains to be studied how large (ℓ,m) are sufficient for
the thermodynamic limit of the model.
Acknowledgments
This work is supported in part by a Grant-in-Aid for
Scientific Research from Japan Society for the Promotion
of Science.
[1] E. Cerda and L. Mahadevan, Phys. Rev. Lett. 80, (1998)
2358.
[2] R. da Silveira, S. Chaieb and L. Mahadevan, Science,
287, (2000) 1468.
[3] Sahraoui Chaieb, Vinay K. Natrajan, and Ahmed Abd
El-rahman, Phys. Rev. Lett. 96, 078101(1 - 4) (2006).
[4] W. Helfrich, Z. Naturforsch, 28c (1973) 693.
[5] A.M. Polyakov, Nucl. Phys. B 268 (1986) 406.
[6] H. Kleinert, Phys. Lett. 174B (1986) 335.
[7] D. Nelson, in Statistical Mechanics of Membranes and
Surfaces, Second Edition, edited by D. Nelson, T.Piran,
and S.Weinberg, (World Scientific, 2004), p.1.
[8] F. David, in Two dimensional quantum gravity and ran-
dom surfaces, Vol.8, edited by D. Nelson, T. Piran, and
S. Weinberg, (World Scientific, Singapore, 1989), p.81.
[9] D. Nelson, in Statistical Mechanics of Membranes and
Surfaces, Second Edition, edited by D. Nelson, T.Piran,
and S.Weinberg, (World Scientific, 2004), p.149.
[10] K. Wiese, in: C.Domb, J.Lebowitz (Eds.), Phase Transi-
tions and Critical Phenomena, Vol. 19, Academic Press,
London, 2000, p.253.
[11] M. Bowick and A. Travesset, Phys. Rep. 344 (2001) 255.
[12] G. Gompper and M. Schick, Self-assembling amphiphilic
systems, In Phase Transitions and Critical Phenomena
16, C. Domb and J.L. Lebowitz, Eds. (Academic Press,
1994) p.1.
[13] J.F. Wheater, J. Phys. A Math. Gen. 27 (1994) 3323.
[14] L. Peliti and S. Leibler, Phys. Rev. Lett. 54 (15) (1985)
1690.
[15] F. David and E. Guitter, Europhys. Lett, 5 (8) (1988)
[16] M. Paczuski, M. Kardar, and D. R. Nelson, Phys. Rev.
Lett. 60 (1988) 2638.
[17] M.E.S. Borelli, H. Kleinert, and Adriaan M.J. Schakel,
Phys. Lett. A 267 (2000) 201.
[18] M.E.S. Borelli and H. Kleinert, Phys. Rev. B 63 (2001)
205414.
[19] Y. Kantor and D.R. Nelson, Phys. Rev. A 36 (1987) 4020.
[20] J.F. Wheater, Nucl. Phys. B 458 (1996) 671
[21] M. Bowick, S. Catterall, M. Falcioni, G. Thorleifsson,
and K. Anagnostopoulos, J. Phys. I France 6 (1996) 1321;
M. Bowick, S. Catterall, M. Falcioni, G. Thorleifsson,
and K. Anagnostopoulos, Nucl. Phys. Proc. Suppl. 47
(1996) 838;
M. Bowick, S. Catterall, M. Falcioni, G. Thorleifsson,
and K. Anagnostopoulos, Nucl. Phys. Proc. Suppl. 53
(1997) 746.
[22] A.Baumgartner and J.S.Ho, Phys. Rev. A, 41, (1990)
5747 .
[23] S.M. Catterall, Phys. Lett. 220B, 253 (1989).
[24] S.M. Catterall, J.B. Kogut, and R.L. Renken, Nucl. Phys.
Proc. Suppl. B 99A (1991) 1.
[25] J. Ambjorn, A. Irback, J. Jurkiewicz, and B. Petersson,
Nucl. Phys. B 393 (1993) 571.
[26] K. Anagnostopoulos, M. Bowick, P. Gottington, M. Fal-
cioni, L. Han, G. Harris, and E. Marinari, Phys. Lett.
317B (1993) 102.
[27] M. Bowick, P. Coddington, L. Han, G. Harris, and E.
Marinari, Nucl. Phys. Proc. Suppl. 30 (1993) 795;
M. Bowick, P. Coddington, L. Han, G. Harris, and E.
Marinari, Nucl. Phys. B 394 (1993) 791.
[28] H. Koibuchi, Phys. Lett. A 300 (2002) 582;
H. Koibuchi, N. Kusano, A. Nidaira, K. Suzuki, and
M.Yamada, Phys. Lett. A 319 (2003) 44;
H. Koibuchi, N. Kusano, A. Nidaira, and K. Suzuki,
Phys. Lett. A 332 (2004) 141.
[29] H. Koibuchi, Eur. Phys. J. B 45 (2005) 377; Eur. Phys.
J. B, 52 (2006) 265.
[30] H. Koibuchi, A. Nidaira, T. Morita, and K. Suzuki, Phys.
Rev. E 68 (2003) 011804;
H. Koibuchi, Z. Sasaki, and K. Shinohara, Phys. Rev. E
70, (2004) 066144.
[31] J-P. Kownacki and H. T. Diep, Phys. Rev. E 66 (2002)
066105.
[32] H. Koibuchi, N. Kusano, A. Nidaira, K. Suzuki, and M.
Yamada, Phys. Rev. E 69 (2004) 066139;
H. Koibuchi and T. Kuwahata, Phys. Rev. E 72 (2005)
026124;
I. Endo and H. Koibuchi, Nucl. Phys. B 732 [FS] (2006)
[33] Ling Miao, Udo Seifert, Michael Wortis, and Hans-
Gunther Dobereiner Phys. Rev. E 49, 5389 - 5407 (1994).
[34] Marija Jari, Udo Seifert, Wolfgang Wintz, and Michael
Wortis, Phys. Rev. E 52 (1995) 6623 - 6634.
[35] E.Helfer, S.Harlepp, L.Bourdieu, J.Robert,
F.C.MacKintosh, and D. Chatenay, Phys. Rev. Lett. 87
(2001) 088103.
[36] K. Murase, T. Fujiwara, Y. Umehara, K. Suzuki, R.
Iino, H. Yamashita, M. Saito, H. Murakoshi, K. Rito-
hie, and A. Kusumi, Ultrafine Membrane Compartments
for Molecular Diffusion as Revealed by Single Molecule
Techniques, Biol. J. 86 (2004) 4075 - 4093 .
[37] H.Koibuchi, submitted; the first model and the second
model are identical to that in cond-mat/0607225 and that
in cond-mat/0607508, respectively.
[38] H.Koibuchi, J. Stat. Phys. in press, cond-mat/0607225.
[39] G. Grest, J. Phys. I (France) 1 (1991) 1695.
[40] M. Bowick and A. Travesset, Eur. Phys. J. E 5 (2001)
[41] M. Bowick, A. Cacciuto, G. Thorleifsson, and A. Traves-
set, Phys. Rev. Lett. 87 (2001) 148103.
[42] H. Noguchi and G. Gompper, Phys. Rev. Lett. 93,
258102 (2004).
[43] F. David, Nucl. Phys. B 257 [FS14], 543 (1985).
[44] D.V. Boulatov, V.A. Kazakov, I.K. Kostov and A.A.
Migdal, Nucl. Phys. B 275 [FS17], 641 (1986).
[45] M. Matsumoto and T. Nishimura, ”Mersenne Twister:
A 623-dimensionally equidistributed uniform pseudoran-
dom number generator”, ACM Trans. on Modeling and
Computer Simulation Vol. 8, No. 1, January (1998) pp.3-
http://arxiv.org/abs/cond-mat/0607225
http://arxiv.org/abs/cond-mat/0607508
http://arxiv.org/abs/cond-mat/0607225
