PHYSICAL REVIEW E 86, 061401 (2012)

Viscoelasticity of two-layer vesicles in solution
C.-Y. David Lu*
Department of Chemistry and Department of Physics, Center of Theoretical Physics, National Taiwan University, Taipei 106, Taiwan

Shigeyuki Komura
Department of Chemistry, Graduate School of Science and Engineering, Tokyo Metropolitan University, Tokyo 192-0397, Japan

Kazuhiko Seki
National Institute of Advanced Industrial Science and Technology (AIST), Tsukuba 305-8565, Japan(Received 15 August 2012; published 5 December 2012)
The dynamic shape relaxation of the two-layer vesicle is calculated. In addition to the undulation relaxation
where the two bilayers move in the same direction, the squeezing mode appears when the gap between the twobilayers is small. At a large gap, the inner vesicle relaxes much faster, whereas the slow mode is mainly due tothe outer-layer relaxation. We have calculated the viscoelasticity of the dilute two-layer-vesicle suspension. It
is found that for a small gap, the applied shear drives the undulation mode strongly while the slow squeezing
mode is not much excited. In this limit, the complex viscosity is dominated by the fast-mode contribution. On theother hand, the slow mode is strongly driven by shear for a larger gap. We have determined the crossover gap,
which depends on the interaction between the two bilayers. For a series of samples where the gap is changed
systematically, it is possible to observe the two amplitude switchings.
DOI: 10.1103/PhysRevE.86.061401

PACS number(s): 82.70.Uv, 83.60.Bc, 87.16.D−, 83.80.Qr

I. INTRODUCTION

Vesicle dynamics has long been investigated both experimentally and theoretically [1–6]. The decay rate of thethermally excited shape fluctuation provides information ofvesicle elasticity. The knowledge of the microscopic relaxation
can also be used to predict the macroscopic rheologicalproperty of the vesicle solution [6].
The vesicles are mostly nonequilibrium system. The external perturbation, be it thermal, electrical, sonication [7,8], orflow [9,10], transforms the lamellar structure into the vesicles.
In such situation, one often obtains mixture of unilamellar
vesicle and multilamellar vesicles (MLV) of various sizes.
However, in sharp contrast with the very detailed calculations
and scattering experiments on the unilamellar vesicles, there
are relatively few studies of the MLV dynamics in the literature.
In this paper, we consider the two-bilayer vesicle whichis the simplest MLV. We calculate the dynamic relaxation
rates, as well as the rheological response of the two-layer
vesicles. From the experimental point of view, it is hard to
prepare vesicles with exactly two bilayers. Nonetheless, sucha concrete calculation provides a clear picture of how theinteraction between the two membranes affect the relaxation
rates, as well as the squeezing (lubrication) flow which arises
exclusively in MLV. Since our method can be extend to vesicles
with more than two bilayers, a more specific calculation canbe performed similarly.
The squeezing dynamics of the sandwiched solvent between the two layers has been analyzed in two related
problems. This relaxation process of the flat lamellar system
was calculated by Brochard and de Gennes as the “slip*cydlu@ntu.edu.tw

1539-3755/2012/86(6)/061401(11)mode” [11] or later measured and called the “baroclinicmode” [12–14]. In the soap film system, the “squeezing mode”
dispersion has been calculated and measured [15–17]. In thiswork, we obtain a similar relaxation, in which the solvent
is squeezed between the two adjacent bilayers to relax thebilayer curvature energy and the mutual interaction energy
between the two bilayers. As this mode arises only in MLV,
we are interested in its dispersion relation and its couplingstrength with the applied shear. Due to the strong lubrication
resistance, the squeezing mode is often the slowest relaxation
mode. This makes the squeezing mode an important candidate
for the vesicle rheology. However, not every relaxation modeis equally excited by the applied shear. Therefore whether the
shear can drive the squeezing mode with a large amplitude is
an important problem to consider. In this work, we will try tobuild our understanding of the coupling strength through theconcrete calculation.
Based on the analysis below, we find that the squeezing
mode makes the dominant contribution to the complex viscosity for strongly interacting bilayers. Interestingly, when the
gap between the two bilayers is smaller than a characteristic
crossover gap, the fast undulation mode becomes the dominantmode for the complex viscosity. We find that the crossovergap gets very small when the bilayer interaction is strong.
On the opposite limit where the bilayer has extremely weakinteraction, the crossover gap becomes comparable to the innervesicle radius.
In Sec. II, we define our model. Section III discusses
the elastic force. In Sec. IV and the appendices, the flow
resistance between the bilayers is solved. In Sec. V, weshall analyze the relaxation rate. In Sec. VI, we consider theviscoelasticity of the dilute vesicle suspension. Section VII
includes further discussions of the results. Finally, inSec. VIII, we provide a nontechnical summary of ourresults.

061401-1©2012 American Physical Society

C.-Y. DAVID LU, SHIGEYUKI KOMURA, AND KAZUHIKO SEKI

FIG. 1. The vesicle consists of two bilayers with the radius r1 andr2 . The surrounding solvent viscosity is η everywhere.

PHYSICAL REVIEW E 86, 061401 (2012)of the strain. At large d, the curvatures and the areas forthe two bilayers differ markedly. Presumably, a microscopic
model may derive a more suitable weighted interaction energy,
which is proportional to [u2 − (r1 /r2 )β u1 ]2 with a weighting
exponent β. Below we present the calculation without thisextra weighting factor. Nonetheless, the calculation procedure
is exactly the same for the case β = 0.
At the dilute phase boundary of the lamellar phase, thesurface tension γn should vanish [18]. Inside the lamellar
phase, the tension depends on the applied osmotic pressure. To
simplify the discussion, we shall only consider the case where
the surface tension vanishes γn = 0. In fact, the calculation
with finite surface tension can be carried out in the same way.
The layer displacement un is related to the radial velocity vr
at r = rn by

II. THE MODEL

We consider a two-layer vesicle with two bilayers locatedat the spheres with mean radius r1 and r2 as shown in Fig. 1.
The deformation of the bilayer shapes are described by theradial layer displacements u1 and u2 relative to the two
reference spheres, respectively. Both of the bilayer membranes
are surrounded by a solvent of viscosity η. Let a denote thearea per molecule projected on the reference sphere and a0 theaveraged projected area per molecule. For each membrane,
we define the dimensionless surface density φn (n = 1,2)
given by the ratio a0 /a. Notice that φn becomes unity atequilibrium. In a fluctuating vesicle, φn is not uniform ingeneral. In this work we propose the free energy of a two-layer
vesicle as
  
 3

r2 − r13 B
u2 − u1 2
d
F =
3
2
r2 − r12  κEγn + Hn2 + (φn − 1)2 dAn , (1)
+
2
2
n=1
where B is the layer compression modulus, γn the surface
tension, κ the bending modulus, and E the area stretching(compression) modulus. Here both κ and E are taken to be thesame for the two membranes. Moreover, d is the differential
solid angle, and the surface area element is approximately
given by dAn ≈ [1 + (∇⊥ un )2 /2]rn2 d. The mean curvature
Hn is given by
Hn ≈ −

22un+ 2 + ∇⊥2 un ,
rn
rn

∂un
= vr (rn ),∂twhere we neglect the solvent permeation. The surface densityobeys
2vr (rn )∂φn=−φn − ∇⊥ · [v⊥ (rn )φn ],
∂t
rn

(4)

where ∇⊥ is the two-dimensional (2D) surface derivativeand v⊥ (rn ) is the tangential velocity at r = rn . Here wehave neglected the amphiphile exchange flux between theneighboring bilayer (the amphiphile permeation).
The flow field obeys the Stokes equation, which is presented
here as the formη∇ 2 v − ∇p = 0,(5)where η is the viscosity and p the pressure. Here we write
down the force balance conditions that are satisfied on the twolayers. The normal force balance is given by−δF+ σrr (rn+ ) − σrr (rn− ) = 0,
δun

(6)

where σrr = −p + 2η∂r vr (r being the radial distance) and thesuperscripts + and − indicate that the stresses are evaluatedat the exterior and interior of the bilayer, respectively. On theother hand, the tangential force balance is given by(2)

up to linear order in un .
Two comments should be made about the first term ofEq. (1). The compression modulus B should be a functionof r1 and r2 . In principle, a microscopic statistical model
for MLV should provide the functional form. Here we focus
on the layer dynamics; to build such a model is beyond
the scope of this work. As a rough estimate for discussion,
below in Sec. III B, we will use the B, which is calculated
from the flat layers, and replace the distance d between twolayers by r2 − r1 . This approximation is justified for d  r2 ,while it may deviate considerably when d is comparable
to r2 .
The interaction term proposed here is proportional to(u2 − u1 )2 , which arises naturally at small d from the square

(3)

−

δF
+ σ⊥ (rn+ ) − σ⊥ (rn− ) = 0,
δxn

(7)

where xn is the 2D tangent displacement of the layers, and the
2D stress is defined as σ⊥ = θ̂ σrθ + ϕ̂σrϕ (θ and ϕ being thepolar and azimuthal angles, respectively).
The viscoelastic response of a dilute two-layer-vesicle
suspension can be calculated by considering the stress responseto an external flow at large distances from the vesiclesv∞ (r,t) = ∇[r 2 Y20 (θ,ϕ)]eiωt ,(8)where is the strength of the elongational flow, Ylm (θ,ϕ) are
the spherical harmonics, and ω is the angular frequency. Eachsuspended vesicle contributes to the averaged stress. In the
dilute solution, the effective complex viscosity is calculated

061401-2

VISCOELASTICITY OF TWO-LAYER VESICLES IN SOLUTION

as [6,19]


II
p20
∗
η =η 1−
φv ,4η r23(9)where φv = (4π/3)cr23 is the volume fraction occupied by theIIisvesicles (c is the number density of the vesicles) and p20
the frequency-dependent complex coefficient of the sphericalharmonic expansion of the pressure [see Eqs. (A10) and (B6)
later]. The effective viscosity can also be expressed in termsof the complex modulus as G∗ = iωη∗ .
III. THE ELASTIC FORCES
A. Force expressions

The elastic forces are calculated by evaluating the derivatives −δF /δun and −δF /δxn . From Eq. (4), we see that thebilayer tangential displacement and the radial displacement
produce the first-order surface density perturbation asδφn = −

2un
− ∇⊥ · xn .
rn

(10)

PHYSICAL REVIEW E 86, 061401 (2012)

together with the counterion and coion entropy produce thefree energy per area [20–22],qψ
64Cs kB T
e−κD d ,
tanh2
(15)
Ve =
κD
4kB T
where Cs is the salt concentration, qψ the potential energy of
the counterion q at the surface, κD the inverse of the Debye
screening length, d = r2 − r1 the distance between the two
surfaces, and kB T the thermal energy.
The van der Waals attraction potential per unit area can becalculated by summing the dipoles to obtain

1
A
1
2
(16)
VvdW −+− 2 ,
12π (d − δ)2(d + δ)2d
where A is the Hamaker constant and δ is the membranethickness. A more complicated Lifschitz theory calculation
can provide a more accurate description [20–23]. The sumof these two interactions consists the standard DLVO theory.
When the electrostatic repulsion and van der Waals attraction
stabilize MLV, B can be evaluated from
B=d

The tangential elastic force is given by−δF= ∇⊥ (δγn ),
δxn

(11)

where the tension perturbation δγn is induced by the surface
density perturbation, i.e., δγn = −E(φn − 1). Notice that the
2D gradient has the usual component form,∇⊥ (δγn ) = θ̂

1 ∂(δγn )1 ∂(δγn )
+ ϕ̂
.
r ∂θ
r sin θ ∂ϕ(12)

Up to linear order in un and δγn , the normal force on bilayer
1 is given by−

B(r23 − r13 ) u2 − u12(γ1 + δγ1 )
δF
=
−
2
2
δu1
3(r2 − r1 )
r1
r

1
2
+ (γ1 − κ∇⊥2 ) ∇⊥2 + 2 u1 ,
r1whereas that for bilayer 2 is
B r23 − r13 u2 − u12(γ2 + δγ2 )
δF
=−
−
−
δu2
3(r2 − r1 )2 r22
r22+ (γ2 − κ∇⊥2 ) ∇⊥2 + 2 u2 .
r2

∂ 2 (Ve + VvdW )
.
∂d 2

(17)

For a flexible bilayer where the undulation entropy depends
strongly on the membrane separation d, Helfrich estimated thefree energy per area and B by a self-consistent argument toobtain [21,24](kB T )2
,
(18)
κd 3
where c0 is an order unity numerical constant. In the following
examples below, we use c0 = 36/π 2 . When the van der Waals
interaction becomes important, a more elaborate calculationis required to combine the undulation entropy and the van der
Waals attraction [25].
B = c0

IV. THE NORMAL FORCE BALANCE

(13)

(14)

In the large stretching modulus limit (E → ∞), the tension
perturbations δγn become Lagrange multipliers, so they ensurethat the right-hand side of Eq. (10) vanishes. By using Eqs. (7)and (11), the values of δγn are determined from the viscousstress σ⊥ (see Appendix C).
B. The bilayer interactions

There are several interactions which contribute to the layer
compression modulus B. In this subsection, we indicate their
physical origins and give the simplest formulae to describe
them. For charged bilayers, the electrostatic interaction,

Using the spherical harmonic expansion, the Stokes equation can be solved in terms of the radial functions, which aresimple polynomial of the radius. The detailed calculation ispresented in Appendix A. In general, the solution depends
on the boundary velocities, which are the velocities of themembranes and the applied external shear. In the usual casewhere the bilayers have a large stretching modulus E, we cansimplify the discussion by considering the large stretchingmodulus limit E → ∞. In this limit, the surface density
approaches unity, and the combinations E(φn − 1) become
the Lagrange multipliers to ensure that surface densities areconstant. Hence, the constant E can be eliminated. The detailed
calculation is presented in Appendix B. At the boundaries, thislimit also turns the tangential velocity into the function of thenormal velocity. This means that the normal stresses σrr (rn+ )and σrr (rn− ), as well as the tension perturbation δγn , are all
proportional to the bilayer normal velocities and the external
shear. The detailed calculation is presented in Appendix C.
The bilayer displacements and velocities are expanded
as the sum of the spherical harmonics Ylm (θ,ϕ) times thetime-varying amplitudes, i.e., un (θ,ϕ) = un Ylm (θ,ϕ) and

061401-3

C.-Y. DAVID LU, SHIGEYUKI KOMURA, AND KAZUHIKO SEKI

vr (rn ,θ,ϕ) = vn Ylm (θ,ϕ). Without cluttering the notation withthe angular quantum numbers l and m, hereafter we use un
and vn to express the time-varying amplitudes for an arbitrary
set of (l,m). After performing some calculations to express
σrr and δγn by vn [see Eqs. (C1), (C5), (B4), (B6), and (B9)],
we find that the normal force balance Eq. (6) at r1 and r2 canbe written asE·

r12 u1
r22 u2+D·

r12 v1
r22 v2= 20η eiωt δl2 δm0 ê2 ,(19)where ê2 = (0,1). Here we prefer to use the variables r12 u1
and r22 u2 which make the matrices E and D symmetric. With
this variable choice, the free energy per solid angle is thenthe quadratic form of E. The components of the matrix E aregiven by 2 2
γ1 r2 ρ + κ L̂2 (L̂2 − 2)B(1 − ρ 3 )+,
E11 = 3
3r2 (1 − ρ)2 ρ 4
r26 ρ 6B(1 − ρ 3 ), E21 = E12 ,3r23 (1 − ρ)2 ρ 2
 2
γ2 r2 + κ L̂2 (L̂2 − 2)B(1 − ρ 3 )+,
E22 = 3
3r2 (1 − ρ)2r26

E12 = −

where have defined ρ = r1 /r2 ( 1) and


∂
1 ∂2
1 ∂
L̂2 = −sin θ−,
sin θ ∂θ∂θ
sin2 θ ∂ϕ 2

(20)

(21)

+ (l 2 − 1)(2l − 1)(2l + 3)ρ 2l−1
+ (l − 1)2 (4l 2 + 8l + 3)ρ 2l−3 + (8l 2 + 8l − 4)ρ −2 ],ηr 4l+1 (4l + 2)[(2l 3 + 5l 2 + l − 2)(ρ 3l+2 − ρ l−1 )D12 = 2 2G0 (l + l)(25)

where dRj are the (nonorthogonal) right eigenvectors. To workwith the orthogonal eigenvectors, here we define the symmetricdecay rate matrix,
M = D−1/2 · E · D−1/2 ,(26)

which has the (normalized) eigenvectors dj ,M · dj = j dj .

(27)

(28)

which will be used in the later discussion. Hereafter, the fasterand the slower rates are denoted by 1 and 2 , respectively.
The full expression of the decay rates are straightforward
but complicated. Hence, we shall obtain some simplified
expressions to gain some understanding on the nature of
the relaxation. We, first, discuss the fast mode 1 . When
ρ = 1, the fast undulation mode has a well-known dispersion
relation [3,4,6],2κ(l − 1)l 2 (l + 1)2 (l + 2),
ηr23 (2l + 1)(2l 2 + 2l − 1)

(29)

E11
.
(30)
D11
This corresponds to a simple picture that the fast mode consists
mainly the inner bilayer relaxation, whereas the outer bilayer
does not move much.
Next we consider the slow mode 2 . When ρ ≈ 1, the slow
mode can be approximately obtained by the series expansionof 1 − ρ. Based on the relation 1 + 2 = tr(D−1 · E), in
which 2 only starts from the second-order term, we get1 tr(D−1 · E) for the zeroth and the first-order terms.
Because the product 1 2 is given by det(D−1 · E), we use
the leading two terms of det(D−1 · E)/ 1 to obtain1 ≈+ 2l(l + 2)(2l − 1)(2l + 3)ρ 2l+2
(22)

where
G0 = r24l+4 [4ρ + 4ρ 4l+3 − (2l + 1)2 ρ 2l+4(23)

Notice that the components of D are the product of η/r23 andthe functions of ρ and l.
When the two bilayers are well separated, i.e., r1  r2 , they
are not hydrodynamically coupled. In this limit, D12 and D21
become small, and we recover the isolated vesicle damping
given by [4–6]η(2l + 1)(2l 2 + 2l − 1).
rn3 l(l + 1)

D−1 · E · dRj = j dRj ,where the bending modulus is doubled as there are twobilayers. For ρ slightly less than unity, a series expansion
can be made if desired. For ρ  1, the decay rate can beapproximated by

D21 = D12 ,ηr 4l+1 (2l + 1)[−(l + 2)2 (4l 2 − 1)ρ 2l+4
D22 = 2 2G0 (l + l)

Dnn =The relaxation rates, denoted as j (j = 1,2), are theeigenvalues of the matrix D−1 · E. The normal force balance
condition Eq. (6) on the two bilayers gives the eigenvector
equation0 =+ (2l 3 + l 2 − 3l)(ρ l+1 − ρ 3l )],− (6 − 8l − 8l 2 )ρ 2l+2 − (2l + 1)2 ρ 2l ].

V. RELAXATION SPECTRUM

M = 1 d1 d1 + 2 d2 d2 ,

ηr24l+1 (2l + 1)[−(l + 1)2 (4l 2 − 1)ρ 2l+1G0 (l 2 + l)

+ l 2 (4l 2 + 8l + 3)ρ 2l + (8l 2 + 8l − 4)ρ],For the more general bilayer interaction which is proportional to B[u2 − (r1 /r2 )β u1 ]2 , the similar calculation will giverise to an extra factor ρ 2β for the B term in E11 , and an extra
factor ρ β for the B terms in E12 and E21 . The form of E22 is
not affected.

Here the eigenvalues are the same as in Eq. (25), and dj is
proportional to D1/2 dRj . We can decompose the decay ratematrix as

whereas the components of the matrix D are
D11 =PHYSICAL REVIEW E 86, 061401 (2012)(24)061401-42 ≈Bl(l + 1)(1 − ρ)2 (2 − ρ)12ηκ(l − 1)l 2 (l + 1)2 (l + 2)
(1 − ρ)3 .
+24ηr23(31)

VISCOELASTICITY OF TWO-LAYER VESICLES IN SOLUTIONPHYSICAL REVIEW E 86, 061401 (2012)

The leading term (B/12η)(1 − ρ)2 l(l + 1) is similar to theslip mode of the planar smectic, which has the decay rate
(B/12η)d 2 q⊥2 , where q⊥ is the wave vector projected onto thebilayer plane [11]. This expression is also analogous to thesqueezing-mode dispersion obtained for soap film [15,16]. In
the other extreme of ρ = 0, the slow mode has the dispersion
such thatthe numerical calculated decay rates. The upper one is the fastmode 1 , which has the limit 0 at ρ = 1. It coincides with theapproximation Eq. (30) (dotted line) for ρ < 0.6. The lower
solid line represents the slow mode 2 . The approximate slow
rate Eq. (34) follows qualitatively the exact value of 2 forthe full range of ρ. It also coincides with Eq. (32) in the limitof ρ = 0. Incidentally, the inner layer relaxation Eq. (30) alsofits the slow mode 2 at ρ ≈ 1.
Figure 2 is useful to illustrate the nature of the relaxation
for the fast and slow modes. In a real system, B is a function of
r1 and r2 (or ρ and r2 ), which is difficult to keep as a constant.
One will make a plot with a function for B which is suitable
for the specific MLV system. If one uses the flat layer formulas
described in Sec. III B as approximations, caution should be
kept for the accuracy at small ρ. Given an accurate function
for B, together with the proper weighting factor (to combineρ 2β and ρ β factors in E), the above decay rate approximations
should work better.

2 ≈B + 3(l − 1)l(l + 1)(l + 2))κ/r23 l(l + 1).
3η(2l + 1)(2l 2 + 2l − 1)

(32)

If we set the shear perturbation as l = 2, the numerator
contains a factor B + 72κ/r23 . Hence, the dimensionlessparameter Br23 /κ becomes important when its value is muchlarger than 72.
We now discuss another approximate expression for theslow mode valid for the whole range of ρ. Since the sum of
the two decay rates 1 + 2 = tr(D−1 · E) is dominated by
the fast mode for any ρ, the slow mode can be approximatedby the ratio
det(D−1 · E)tr(D−1 · E)
E11 E22 − E12 E21.
=
E11 D22 − E12 D21 + E22 D11 − E21 D12

VI. VISCOELASTICITY

2 ≈

(33)

When the inner bilayer 1 relaxes fast, one can apply the
adiabatic approximation to the fast relaxing inner membrane.
The approximated expression has even a simpler denominatorgiven by2 ≈E22 − E21 (E12 /E11 ).
D22 − D21 (E12 /E11 )(34)

In Fig. 2, we plot the decay rates as a function of ρ = r1 /r2for l = 2, keeping a constant B. The two solid lines represent
8

j

/

0

6

4

2

0
0

0.2

0.4

0.6

0.8

1

FIG. 2. The scaled relaxation rate j / 0 as a function of thedimensionless size ratio ρ = r1 /r2 between the two layers. Here 0is the vesicle relaxation rate with the rigidity 2κ given by Eq. (29). The
two solid lines represent the two vesicle relaxation rates 1 and 2
obtained numerically. The dotted line is an approximation for the fastmode given by Eq. (30). The dashed line represents the approximate
slow relaxation rate 2 given by Eq. (34). In this plot, we set
B = 600 J/m3 , κ = kB T = 4 × 10−21 J, η = 10−3 Pa s, r2 = 10−7 m,
l = 2, so Br23 /κ = 150.

We now consider two-layer vesicles under the external
oscillatory shear flow. Rearranging Eq. (9) together with
Eq. (B6), we obtainη∗ /η − 1
v2
5
= −
.
φv
2 r2 eiωt(35)

It is apparent that the dilute hard-sphere limit is recovered
when v2 = 0. By putting un = vn / iω and solving Eq. (19) for
v2 , we obtain20ηη∗ /η − 1
5
= − iω 3 ê2 · (E + iωD)−1 · ê2 ,
φv
2
r2

(36)

where one should set l = 2 for the matrices E and D. Hereê2 appears twice, because the shear directly affects the bilayer
2 through Eq. (19), and the velocity of the bilayer 2 carries
the stress contribution of the vesicle according to Eqs. (9)and (B6).
Considering the high-frequency limit, we now define
the dimensionless shear coupling strength S and the shear
deformation unit vector ŝ by20η 1/2 −1/2
1/2
S ŝ =
D· ê2 .
(37)
r23
Notice that S and ŝ depend only on the ratio ρ as long as we
fix to l = 2. As shown in Fig. 3, S is weakly dependent onρ because its value varies only between 2.1 and 2.19. Hencethe right-hand side of Eq. (36) varies between 0.31 and 0.4,
whereas the hard-sphere case gives 2.5. This means that, at
high frequencies, the coupling of a vesicle to the flow is weakercompared to that of a hard sphere.
The high-frequency viscosity asymptotically approachesη∞ = η 1 + 52 − S φv = η0 − ηSφv ,(38)
where η0 = η[1 + (5/2)φv ]. Using the eigenmode decomposition I = d1 d1 + d2 d2 , we can separate the viscosity
contributions from each mode as

061401-5

ηj = Sη(ŝ · dj )2 ,(39)

C.-Y. DAVID LU, SHIGEYUKI KOMURA, AND KAZUHIKO SEKI

PHYSICAL REVIEW E 86, 061401 (2012)

2.2

0.75

2.18

0.5

S

/

2.16

0.25

2.14
2.12

0
2.1
0

0.2

0.4

0.6

0.8

0

1

FIG. 3. The shear coupling strength S defined by Eq. (37) as afunction of the dimensionless size ratio ρ between the two layers.

so we have
ηS = η1 + η2 .

(40)

The full expression for the complex viscosity then is given byη∗ (ω) = η0 − φv

2

j =1

ηj

iω
.
iω + j(41)

The relative strength of the two modes are shown in Fig. 4.
Alternatively (but equivalently), the complex modulus can be
expressed asG∗ (ω) = iωη∞ + φv

2

j =1

Gj

ω2 + iωj,ω2 + 2j(42)

where the two modes have the positive amplitudes Gj = j ηj .
Even for different bilayer interaction strength B, we always
find that the slow mode has the larger viscosity amplitudeat ρ ≈ 0, and smaller one at ρ ≈ 1. The latter is reasonablebecause the shear perturb the two bilayers similarly for ρ ≈ 1,1j/( 1+ 2)

0.8
0.6
0.4
0.2
0
0

0.2

0.4

0.6

0.8

0.2

0.4

0.6

0.8

1

FIG. 5. The polar angle α (divided by π ) of various 2D vectors as
a function of the dimensionless size ratio ρ between the two layers.
The solid line represents ŝ, while the dotted and the dashed lines
represent d1 and d2 , respectively. We chose the value Br23 /κ = 150.

where ŝ is along the (1,1) direction. In terms of the polar
angle α between ŝ and the first axes on the (r12 u1 ,r22 u2 )plane, the (1,1) direction corresponds to α = π/4. Notice that
this direction also corresponds to the undulation eigenmodedirection. Since the sum of the two viscosity amplitudes isroughly a constant, as shown in Eq. (40), the slow squeezing
mode takes a small viscosity amplitude η2 for ρ ≈ 1. In theother limit of ρ ≈ 0, the shear mainly perturbs the outer layer.
In this limit, we have ŝ = (0,1) or α = π/2, and the slow modeis due to the outer layer relaxation. As a result, η2 becomes thedominant contribution for ρ ≈ 0. In Fig. 5, we plot the angle α
as a function of ρ. The shear vector ŝ always coincides with thefast mode at ρ = 1 and the slow mode at ρ = 0. Therefore, themode switching always takes place between 0  ρ  1. It isalso worth mentioning that Figs. 4 and 5, like Fig. 2, are plottedby assuming a constant B. These plots are used to analyze thebehavior of the mode amplitudes. The more physical plots willbe the ones using a suitable function for B.
Here we define a crossover size ratio ρ ∗ at which the twomodes have the same viscosity amplitude, i.e., η1 = η2 . In
Fig. 6, we plot the calculated ρ ∗ as a function of Br23 /κ. When
Br23 /κ  1, the crossover happens at ρ ∗ ≈ 0.52. This meansthat for the noninteracting case B = 0, the crossover happens
when d = r2 − r1 is roughly the same as r1 . When Br23 /κ
1,
on the other hand, ρ ∗ approaches unity. As mentioned in
Eq. (32), the parameter Br23 /κ will have a noticeable effectonly when it is greater than 72. At the lower right cornerof Fig. 6, where the layer interaction is strong and the layer
separation is not too small, the slow mode has the dominantviscosity contribution. At the upper left corner where the layerinteraction is relatively weak, the fast mode is more excited ascompared to the slow mode.

1

FIG. 4. The relative viscosity strengths defined by Eq. (39) as afunction of the dimensionless size ratio ρ between the two layers. The
solid and the dashed lines corresponds to η1 and η2 , respectively. We
chose the value Br23 /κ = 150. The crossing of the two modes occurs
at ρ ∗ ≈ 0.733.

VII. DISCUSSION

We have determined the crossover ration ρ ∗ in Fig. 6. The
actual bilayer interaction strength B depends on the separationd = r2 (1 − ρ). For the qualitative discussion at nonsmall ρ,
we use the flat layer results in Sec. III B. In Fig. 7 we

061401-6

VISCOELASTICITY OF TWO-LAYER VESICLES IN SOLUTION

1
0.9
0.8

*

fast mode
0.7
0.6

slow mode

0.5
0.4 -2
10

0

2

10

4

10

10

3

Br2 /
FIG. 6. The crossover size ratio ρ ∗ as a function of the dimensionless bilayer interaction Br23 /κ. For ρ > ρ ∗ , the fast mode is thedominant contribution to the viscosity, while the slower relaxation isdominant for ρ < ρ ∗ .

present both information in the (ρ,Br23 /κ) plot to compare
a series of systems with the same outer bilayer size r2 butwith different size ratios ρ. As one varies ρ, the dimensionlessbilayer interaction strength Br23 /κ changes according to the
DLVO theory Eq. (17) (two convex blue lines) or the Helfrich
repulsion Eq. (18) (two concave red lines). If such lines happen
to cross the line ρ ∗ (thick black line), one expects that thetwo viscosity amplitudes change their relative magnitudes. We
have indicated such scenarios by the dashed lines.
As shown in the convex blue dashed line in Fig. 7, anelectrostatically stabilized system may get large B by having

10

2

3

Br2 /

10

4

10

0

-2

10 0.4

0.5

0.6

0.7

0.8

0.9

1

FIG. 7. (Color online) The dimensionless bilayer interaction
Br23 /κ as a function of the size ratio ρ. The crossover ratio ρ ∗ in Fig. 6
is plotted by the thick solid black line. The two convex blue lines arefrom the DLVO theory Eq. (17) with Cs = 0.01 M, κD = 3 × 10−8 m,
A = 10−21 J, δ = 3 × 10−9 m, and kB T = 4 × 10−21 J. The surface
potentials are 12 and 24 mV for the convex blue thin solid anddashed lines, respectively. The concave red lines are the Helfrich
repulsion Eq. (18) with c0 = 36/π 2 , r2 = 10−7 m. The concave redthin solid and dashed lines are for κ = 0.75kB T = 3 × 10−21 J andκ = 7.5kB T = 3 × 10−20 J, respectively.

PHYSICAL REVIEW E 86, 061401 (2012)a high surface charge density. At small gap, the van der Waals
attraction may also lower B, causing an interesting multiple
crossing. This means that the viscosity amplitudes switch
their magnitudes more than once. For a sterically stabilized
system (concave red lines), we find that the dimensionlessinteraction parameter depends strongly on the bending rigidity
as Br23 /κ ∼ κ −2 [see Eq. (18)]. For a soft surfactant bilayer ofκ = 0.75kB T , we find that the slow mode always dominates
the viscosity. Whereas for a lipid bilayer whose bending
rigidity κ = 7.5kB T is one order of magnitude larger, theconcave red dashed line indicates that the squeezing mode
is not much excited by shear when ρ > 0.526.
Caution should be kept when analyzing the MLV when ρis close to unity. When the gap d = r2 − r1 is smaller than the
distortion amplitude |u1 | + |u2 |, our result should be modified
by the extra steric constraint, r1 + u1  r2 + u2 . Given a small
but finite gap d, our linear calculation should be valid for small
distortion amplitudes, |u1 | + |u2 | < d. This means that thelinear response regime is restricted by the gap. As ρ becomescloser to unity (or d → 0), the linear regime shrinks toward
zero. Nonetheless, at a given ρ (but ρ = 1), our result shouldstill be valid within the linear regime.
In this paper, we have only presented the calculation for
two-layer vesicles. For MLV with more than two layers, thecalculation can be performed in the same procedure but witha greater algebraic complexity [26]. For MLV with N layers,
we expect that there are N relaxation modes. When the gap
is small (rN−1 /rN ≈ 1), we expect that the majority of the N
relaxations to bear some resemblance to the squeezing modeor the spherical version of the “slip mode” [11]. For larger gap,on the other hand, the relaxations of the N layers may decouple
from each other. As for the viscoelasticity, we speculate thatthe viscosity amplitudes of the squeezing modes are small for
weakly interacting MLV. Neutral or weakly charged lipid MLV
should be an interesting system to investigate in this direction.
In polymer rheology calculation, one may consider a step
strain for t  0, where the time relaxation of the stress gives
the relaxation modulus G(t). This quantity can be further
converted to the complex modulus G∗ (ω) by the Fouriertransform. Right after the step strain, one often assumes that thepolymer deforms in an affine way [27]. Do we implicitly use
the affine deformation approximation for the MLV rheology?
A step strain for vesicle solution will induce a short but fast
bilayer movement, where the D terms dominate the left-handside of Eq. (19). Therefore, the initial layer displacements
will be proportional to D−1 · ê2 ∝ (−D12 /D11 ,1), which iscompatible with the incompressible constraints [see Eq. (B1)]
at both layers. In the small gap limit, D−1 · ê2 behaves like(ρ 2 ,1). Compared with the affine deformation un ∝ rn , orin terms of our chosen variables (r12 u1 ,r22 u2 ) ∝ (ρ 3 ,1), it isclear that our theory does not use the affine deformation
approximation.
VIII. SUMMARY

In this work, we have calculated the slow relaxation rates
and the viscoelasticity of a dilute two-layer-vesicle solution.
In the model, the two vesicles are assumed to have thespherical equilibrium shape, which makes the analytic flow
field tractable. We focus our discussion in the case where there

061401-7

C.-Y. DAVID LU, SHIGEYUKI KOMURA, AND KAZUHIKO SEKI

is no osmotic pressure so the tension for both layers vanishes.
We also restrict our discussion to the case where the fluid in
the three compartment has the same viscosity. To relax these
two restrictions does not require a new calculation. Should aneed arise, the present calculation procedure can be carried
to those cases without any technical difficulty. Nonetheless,the commonly encountered situations often obey these two
restrictions.
We have found the following results. (i) When the gap
between two vesicles is small, the two relaxation modes
consist of a combination between the squeezing and undulation
characters. The former involves the opposite layer distortions,where the layer interaction drive the lubrication flow between
the two layers as expressed in Eq. (31). The latter, calculated
in Eq. (29), describes a relaxation driven by the layer bending
elasticity, whereas the lubrication flow is not required. We findthat for MLV with a small gap (ρ ≈ 1), the slowest mode is
the squeezing mode for which the Poiseuille flow mobility
vanishes at the limit ρ → 1. The undulation mode appears to
be faster at a small gap.
(ii) At the other limit where the inner bilayer radius is smallρ ≈ 0, it is easier to understand the relaxation in terms of theoriginal variables u1 and u2 . Here the slow mode becomes therelaxation of the outer bilayer, Eq. (32), and the faster mode is
the relaxation of the inner bilayer, Eq. (30).
(iii) When one decreases ρ, the relaxation spectrum changes
from (i) to (ii). The two dispersions change their characterswith ρ to avoid the relaxation rate crossing. The slow modecan be approximated by Eq. (34) for the full range of ρ, wherean adiabatic approximation for the fast u1 variable is made.
When the gap is mall, the slow mode is the squeezing mode.
As the gap increases, the slow mode eventually becomes theouter layer relaxation mode.
(iv) At the low-frequency limit, the real part of the complex
viscosity approaches that of the corresponding hard-spherereference system. On the other hand, the real part of thecomplex viscosity shows that the high-frequency viscosity
increment is between 12 and 16% of the hard-sphere viscosity
increment. The difference between the two limits comes fromthe contributions of the two relaxation modes discussed above.
(v) For a fixed interaction strength Br23 /κ, we find that
there is a crossover size ratio ρ ∗ which describes the crossover
gap size. The slow mode is more difficult to excite when the
gap is smaller than ρ ∗ . At a small gap, ρ > ρ ∗ , the slow
(squeezing) mode has a small viscosity amplitude, while thefast (undulation) mode has a large viscosity amplitude. For
large gap, ρ < ρ ∗ , on the other hand, the slow mode has thedominant viscosity amplitude. This shows that, for the MLVwith small gaps, the slow mode is not the most important modefor viscoelasticity. In other words, the shear is not efficient toexcite the slow mode when the gap is too small.
(vi) The crossover size ratio ρ ∗ depends on the interaction
between the two bilayers. If Br23 /κ increases, ρ ∗ will increase
from 0.52 toward unity. The interaction parameter B can betuned, e.g., by adding block copolymers [28]. For stronglyinteracting MLV which has large Br23 /κ, the crossover gap isvanishingly small (ρ ∗ → 1) so the slow mode is important for
viscoelasticity. For noninteracting MLV, the crossover gap is
0.48r2 , about half of the outer bilayer radius. In this weakly
or noninteracting case there should be a higher possibility of

PHYSICAL REVIEW E 86, 061401 (2012)

finding MLV whose viscoelasticity is dominated by the fasterinterbilayer relaxation. For neutral bilayers whose interaction
arises from the Helfrich repulsion, the more rigid lipid systemwith a smaller gap might be a good candidate to observe thefast-mode-dominated viscoelasticity, as discussed in Fig. 7.
For charged bilayers, a weakly charged system can also be used
for such study. As suggested in a recent theoretical paper [26],the shear excited squeezing modes are important for MLV size
selection. It will be interesting to study a lamellar system where
the squeezing mode cannot be excited efficiently by shear.
ACKNOWLEDGMENTS

We thank M. Cates and S. Fujii for many useful discussions.
We also acknowledge the support from the National Science
Council of Taiwan and the Center of Theoretical Physics of National Taiwan University. S.K. also acknowledges support by
a Grant-in-Aid for Scientific Research (Grant No. 24540439)
from the MEXT of Japan, and the JSPS Core-to-Core Program“International research network for non-equilibrium dynamics
of soft matter.”
APPENDIX A: SOLUTION OF STOKES EQUATION

For the incompressible solenoidal flow, the velocity can be
expressed asv = ∇ × (∇ψ × r) + ∇ × (ζ r),(A1)

where the scalar functions ψ and ζ are the defining function
for the poloidal and toroidal flow fields, respectively. Notethat our definition of the defining function differs by a factor
r compared with the ones in the book by Chandrasekhar [29].
For an incompressible fluid, the divergence of the pressure
gradient vanishes, i.e.,∇ 2 p = 0.

(A2)

Therefore, the pressure gradient is also a solenoidal fieldand can be described by the above decomposition. Since
the toroidal part can be written as (∇ζ ) × r, the condition ∂θ (∂ϕ p) = ∂ϕ (∂θ p) requires that ∂θ (− sin θ ∂θ ζ ) =∂ϕ (∂ϕ ζ / sin θ ) or L̂2 ζ = 0, where


∂
1 ∂2
1 ∂
2
L̂ = −
sin θ
−
.
(A3)
sin θ ∂θ∂θ
sin2 θ ∂ϕ 2
Hence, the pressure gradient is only poloidal and written as∇p = ∇ × (∇ × r)
= r̂1 ∂ ∂(r)1∂ ∂(r)
L̂2 
+ θ̂
+ ϕ̂
,
r
r ∂θ ∂r
r sin θ ∂ϕ ∂r(A4)

where  is the defining function of ∇p. Comparing the θ̂ andϕ̂ directions, we can set
∂(r)
(A5)
∂r
for these two directions. The Laplace equation for the pressure
Eq. (A2) then implies that


2
1 ∂
L̂2 ∂(r)2 ∂ (r)
−
r
= 0.
(A6)
r 2 ∂r∂r 2
r 2 ∂r061401-8p=

VISCOELASTICITY OF TWO-LAYER VESICLES IN SOLUTIONPHYSICAL REVIEW E 86, 061401 (2012)

Therefore, the radial component of the gradient pressure canbe either expressed as L̂2 /r or ∂r2 (r). The latter expression
is consistent with the identification Eq. (A5).
For the velocity field, we will drop the toroidal part by
setting ζ = 0. This is justified because the tangential force field
is a surface gradient and drives only the poloidal flow. To obtain
the defining function ψ, we take the curl of Eq. (5) to get ∇ ×∇ 2 v = 0. Since ∇ 2 v = −∇ × ∇ × v holds for incompressible
flow, we haveIIII
= 0 and plm= 0 to drop the(0  r  r1 ), we set ψlm
Ifunctions which are singular at r = 0. We denote ψlm
and
I
I
I
plm as ψA and pA , respectively. The 2D incompressibility
condition Eq. (B2) relates the two remaining coefficients as∇ × ∇ × ∇ × ∇ × [∇ × (ψr)] = 0,

(A8)

To find the coupling between the pressure p and the velocity
poloidal function ψ in Eq. (5), we rewrite −∇ 2 v = ∇ × ∇ ×v = ∇ × [∇(−∇ 2 ψ) × r]. Equation (5) then is satisfied whenpAI = −ψCII = −∂pr= ηL̂2 ∇ 2 ψ.
(A9)∂r
Using the standard method of separation of variables, thegeneral solution of Eqs. (A2), (A8), and (A9) areI l
II −l−1plm
Ylm ,r + plmr(A10)p = p0 +
lm


I l
II −l−1
ψlm r + ψlm
Ylm
ψ=
r
lm

+


lm

η(2l + 3)(l − 1)r1−l−1
v1 .
l

r22
1
pCII + r25 δl2 δm0 .
2η(2l − 1)(l + 2)
12

(B4)

(B5)

Using the radial velocity amplitude v2 , we can express thepressure coefficient as
pCII =

η(2l − 1)(l + 2)r2l
v2 − 10ηr23 δl2 δm0 .
l+1

(B6)

For the region between the two bilayers (r1  r < r2 ),the expressions are more complex. The incompressibility
Eq. (B2), evaluated at r1 and r2 , provides two conditions
between the four coefficients,
ψBI = F11 pBI + F12 pBII ,I l+2
II −l+1 
plm
plm
r
r
+
Ylm ,
η(2l + 2)(2l + 3) 2ηl(2l − 1)

(B3)

For the region exterior to the bilayer 2 (r2  r < ∞), theII
coefficients ψlm
are zero, except ψ20, which needs to be chosen
to give the far flow Eq. (8). Comparing the radial componentI= /3. We also set
vr from Eq. (8) and L̂2 ψ/r, we set ψ20
III
plm = 0 so p does not diverge at r → ∞, and denote ψlm
and
II
II
II
plm as ψC and pC , respectively. Then the 2D incompressibility
condition Eq. (B2) relates the two remaining coefficients as− + η∇ 2 ψ = 0.
Here we prefer to use p instead of . Operating the above
equation by L̂2 /r and substituting the combination L̂2 /r by∂r p, we obtain

r12
pI .
2η(2l + 3)(l − 1) A

In terms of the radial velocity amplitude v1 , we can further
express the pressure coefficient as(A7)

where, in the square bracket, an equivalent form of ∇ψ × r isused. As detailed in Ref. [29], each curl will switch poloidal
and toroidal parts. The double curl will preserve the type andmodify the defining function by −∇ 2 . The curl of Eq. (5)becomes ∇ × (r∇ 4 ψ) = 0 or
∇ 4 ψ = 0.

ψAI =

ψBII = F21 pBI + F22 pBII ,

(B7)

where

(A11)

F11 = −

where Ylm (θ,ϕ) are the spherical harmonics.

r22 − r12,2η(2l − 1)(l − 1) r22l+1 − r12l+1

r12l+1 r22l+1 r22 − r12
F21 = −,2η(2l + 3)(l + 2) r22l+1 − r12l+1
F12 =APPENDIX B: THE SURFACE INCOMPRESSIBLE LIMITThe surface incompressibility condition also constrains thevelocity field as
2vr (rn )+ ∇⊥ · v⊥ (rn ) = 0(B1)
r
at both r = r1 and r2 . Because the fluid is incompressible, i.e.,∇ · v = 0, Eq. (B1) can also be expressed as the equivalentform,∂vr (rn )
= 0,
(B2)
∂r
at the bilayers. We prefer this condition because it simplifies
the calculation.
We now consider the relaxation rates of a small perturbation
described by the spherical harmonics Ylm (θ,ϕ). To simplify
the notation, we drop the subscript “lm” of the coefficients
I
II
I
II
, ψlm
, plm
, and plm
. For the region within the bilayer 1
ψlm

r22l+3 − r12l+3,
2η(2l + 3)(l − 1) r22l+1 − r12l+1

F22 = −

(B8)

r12 r22l+1 − r22 r12l+1.

2η(2l − 1)(l + 2) r22l+1 − r12l+1

We prefer to use the variables v1 and v2 instead of pBI and pBII .
The pressure coefficients then can be expressed as
pBI = G11 v1 + G12 v2 ,pBII = G21 v1 + G22 v2 ,
where
G11 =

061401-9η (−8l 2 − 4l + 12)r13l+2 r2
lG0

(B9)


− (4l + 6)r1l+1 r22l l(2l + 1)r12 − (l + 2)(2l − 1)r22 ,C.-Y. DAVID LU, SHIGEYUKI KOMURA, AND KAZUHIKO SEKI

G12 =

η 
(−8l 2 − 4l + 12)r1 r23l+2
lG0− (6l + 4)r12l r2l+1 l(2l + 1)r22 − (l + 2)(2l − 1)r12 ,G21 =(4l − 2)η −(2l 2 + 3l + 1)r13l+4 r22l
(l + 1)G0+ (2l 2 + l − 3)r13l+2 r22l+2 + (2l + 4)r1l+1 r24l+3 ,G22 =(4l − 2)η −(2l 2 + 3l + 1)r23l+4 r12l
(l + 1)G0+ (2l 2 + l − 3)r23l+2 r12l+2 + (2l + 4)r2l+1 r14l+3 ,
(B10)
with
G0 = r24l+4 [4ρ + 4ρ 4l+3 − (2l + 1)2 ρ 2l+4
− (6 − 8l − 8l 2 )ρ 2l+2 − (2l + 1)2 ρ 2l ].

(B11)When the two bilayer are well separated (r1  r2 ), both
G11 and G22 become small. In this limit, G12 becomes thecoefficient of Eq. (B4) with r1 replaced by r2 , whereas G21
becomes the coefficient of Eq. (B6) (without the term) withr2 replaced by r1 .
APPENDIX C: THE STRESS AND THE TENSIONPERTURBATIONThe radial component of the normal stress appears in
Eq. (6). From its component form σrr = −p + 2η∂r vr and the2D incompressibility condition Eq. (B2), it is just the negative
pressure. Therefore, the stress differences at the two bilayers
are
 

σrr (r1+ ) − σrr (r1− ) =pAI − pBI r1l − pBII r1−l−1 Ylm ,lm

σrr (r2+ ) − σrr (r2− ) =pBI r2l + pBII − pCII r2−l−1 Ylm .

lm

(C1)

[1] J. S. Huang, S. T. Milner, B. Farago, and D. Richter, Phys. Rev.
Lett. 59, 2600 (1987).
[2] B. Farago, D. Richter, J. S. Huang, S. A. Safran, and S. T. Milner,
Phys. Rev. Lett. 65, 3348 (1990).
[3] M. B. Schneider, J. T. Jenkins, and W. W. Webb, J. Phys. (France)
45, 1457 (1984).
[4] S. T. Milner and S. A. Safran, Phys. Rev. A 36, 4371
(1987).
[5] A. Onuki and K. Kawasaki, Europhys. Lett. 18, 729 (1992).
[6] K. Seki and S. Komura, Physica A 219, 253 (1995).
[7] N. Robinson, J. Pharm. Pharmacol. 12, 129 (1960).
[8] N. Robinson, J. Pharm. Pharmacol. 12, 193 (1960).
[9] T. D. Le, U. Olsson, K. Mortensen, J. Zipfel, and W. Richtering,Langmuir 17, 999 (2001).
[10] O. Diat and D. Roux, J. Phys. (France) II 3, 9 (1993).
[11] F. Brochard and P.-G. de Gennes, Pramana Suppl. 1, 1 (1975).
[12] F. Nallet, D. Roux, and J. Prost, J. Phys. (France) 50, 3147
(1989).

PHYSICAL REVIEW E 86, 061401 (2012)

The large bilayer stretching modulus E suppresses thesurface density perturbation. For the slow relaxation, we willtake the limit E → ∞ to eliminate the parameter E. Within
this limit, both φ1 and φ2 approach unity, so the tension
perturbations δγ1 and δγ2 become the Lagrange multipliers.
The values of the Lagrange multipliers then are determinedby Eqs. (7) and (11) instead of their original definitions, as
discussed below.
In Eq. (7), we are interested in the difference between thetangential stress across the bilayer. We express the tangential
stress as a function of ψ,
 2
∂ ψ
ψ
.
(C2)
σ⊥ = θ̂ σrθ + ϕ̂σrϕ = η∇⊥ r 2 + (L̂2 − 2)∂rr
We now replace ∂r2 ψ using Eq. (A9), and limit our discussion
to the l = 0 mode. The velocity field ∇ × (∇ψ × r) has
components similar to those of Eq. (A4). We can use
the r component vr = L̂2 ψ/r to eliminate ψ as (r/L̂2 )vr
so
 2


r ∂p
4
r ∂vr
vr − 2η
+η 2−
,
(C3)
σ⊥ = ∇⊥
2
2
L̂ ∂r
L̂
L̂2 ∂rwhere L̂2 → l(l + 1) for spherical harmonics with nonzero l.
Because vr is continuous across the bilayer and ∂r vr vanishes
on the bilayer, only the first term can differ across the bilayer.
This is the only important term for the tension perturbation in
Eqs. (7) and (11),  
rn2 ∂prn2 ∂p+.
(C4)δγn = −
L̂2 ∂r rn+ L̂2 ∂r rn−
We then obtain  pI − pI
pII
A
B 2l+1
r1 + B r1−l Ylm ,
δγ1 =
l+1
l
lm


 pI
pCII − pBII −l2l−1
Brr2 Ylm .
+
δγ2 =
l+1 2
l
lm

(C5)

[13] S. Ramaswamy, J. Prost, W. Cai, and T. C. Lubensky, Europhys.
Lett. 23, 271 (1993).
[14] G. Sigaud, C. W. Garland, H. T. Nguyen, D. Roux, and S. T.
Milner, J. Phys. (France) II 3, 1343 (1993).
[15] H. M. Fijnaut and J. G. H. Joosten, J. Chem. Phys. 69, 1022
(1978).
[16] C. Y. Young and N. A. Clark, J. Chem. Phys. 74, 4171 (1981).
[17] E. I. Kats, V. V. Lebedev, and S. V. Malinin, JETP 86, 1149
(1998).
[18] An equilibrium lamellar phase does not have tension. For alamellar phase that could swell, a metastable MLV will take upsome solvent so tension is developed. For a lamellar phase at
the dilute phase boundary, the layers can coexist with the excess
solvent without further swelling. Therefore, the metastable MLV
there is likely to have a negligible tension.
[19] G. K. Batchelor, J. Fluid Mech. 41, 545 (1970).
[20] J. N. Israelachvili, Intermolecular and Surface Forces, 2nd ed.
(Academic Press, New York, 1992).

061401-10

VISCOELASTICITY OF TWO-LAYER VESICLES IN SOLUTIONPHYSICAL REVIEW E 86, 061401 (2012)[21] S. A. Safran, Statistical Thermodynamics of Surfaces, Interfaces,and Membranes (Addison-Wesley, Reading, MA, 1994).
[22] W. A. Russel, D. A. Saville, and W. R. Schowalter, Colloidal Dispersions (Cambridge University Press, Cambridge, UK, 1989).
[23] V. A. Parseigian, Van der Waals Forces: A Handbook for
Biologists, Chemists, Engineers, and Physicists (Cambridge
University Press, Cambridge, UK, 2005).
[24] W. Helfrich, Z. Naturforsch. 33a, 305 (1978).

[25] R. Lipowsky and S. Leibler, Phys. Rev. Lett. 56, 2541 (1986).
[26] C.-Y. D. Lu, Phys. Rev. Lett. 109, 128304 (2012).
[27] M. Doi and S. F. Edwards, The Theory of Polymer Dynamics(Clarendon Press, Oxford, 1986).
[28] S. Fujii, S. Koschoreck, P. Lindner, and W. Richtering, Langmuir25, 5476 (2009).
[29] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability(Clarendon, Oxford, 1961).

061401-11