Eur. Phys. J. E (2016) 39: 52
DOI 10.1140/epje/i2016-16052-3

THE EUROPEAN
PHYSICAL JOURNAL E

Regular Article

Relaxation dynamics of two-component ﬂuid bilayer membranes
Ryuichi Okamoto1,2,a , Yuichi Kanemori1 , Shigeyuki Komura1,2 , and Jean-Baptiste Fournier2
1

Department of Chemistry, Graduate School of Science and Engineering, Tokyo Metropolitan University, Tokyo 192-0397,
Japan
2
Université Paris Diderot, Sorbonne Paris Cité, Laboratoire Matière et Systèmes Complexes (MSC), UMR 7057 CNRS, F-75205
Paris, France
Received 4 November 2015 and Received in ﬁnal form 8 March 2016c EDP Sciences / Società Italiana di Fisica / Springer-Verlag 2016
Published online: 6 May 2016 – 
Abstract. We theoretically investigate the relaxation dynamics of a nearly ﬂat binary lipid bilayer membrane by taking into account the membrane tension, hydrodynamics of the surrounding ﬂuid, intermonolayer friction and mutual diﬀusion. Mutual diﬀusion is the collective irreversible process that leads
to homogenization of the density diﬀerence between the two lipid species. We ﬁnd that two relaxation
modes associated with the mutual diﬀusion appear in addition to the three previously discussed relaxation
modes reﬂecting the bending and compression of the membrane. Because of the symmetry, only one of thetwo diﬀusive modes is coupled to the bending mode. The two diﬀusive modes are much slower than thebending and compression modes in the entire realistic wave number range. This means that the long time
relaxation behavior is dominated by the mutual diﬀusion in binary membranes. The two diﬀusive modes
become even slower in the vicinity of the unstable region towards phase separation, while the other modes
are almost unchanged. In short time scales, on the other hand, the lipid composition heterogeneity inducesin-plane compression and bending of the bilayer.

1 Introduction
Much attention has been paid to artiﬁcial lipid bilayer
membranes as model systems of biological cell membranes [1]. They exhibit a wide variety of complex phenomena in both statics and dynamics, since lipid densities,
membrane deformation and surrounding ﬂuids are coupled
to each other [2]. Dynamical properties of lipid membranes
near the equilibrium are characterized by the wave number dependent relaxation rates. In the early theoreticallystudies, the relaxation rate of a single-component membrane was discussed by regarding the membrane as anelastic sheet having out-of-plane deformation, and further
surrounded by a three-dimensional (3D) ﬂuid. Neglecting the bilayer structure, several authors predicted thatthe relaxation of the bending mode is dominated by thebending rigidity and the viscosity of the surrounding bulk
ﬂuid [3, 4].
Later, Seifert and Langer considered the inter-monolayer friction and the two-dimensional (2D) hydrodynamics of each monolayer, and obtained another relaxation
mode associated with the density diﬀerence between thetwo monolayers [5]. They found that the relaxation of thedensity ﬂuctuation is dominated by the inter-monolayer
friction and is relevant to the slow dynamics charactera

e-mail: okamotor@tmu.ac.jp

ized by large wave numbers, whereas the relaxation of thebending mode is relevant for small wave numbers if themembrane surface tension is not acting. A somewhat similar theory was also developed in ref. [6]. The predicted
mode crossing behavior has been supported by several experiments [7,8] and by molecular dynamics simulations [9].
More recently, some experiments reported a chemicallyinduced tubule growing from a giant unilamellar vesicle(GUV) [10, 11]. In these studies, they showed that theinterplay between the faster bending relaxation and theslower density relaxation on the scale of tens of micrometers plays an essential role.
In recent years, both the statics and dynamics of multicomponent lipid membranes have been extensively studied
because 2D phase separation takes place in a certain rangeof temperature and composition [12–14]. Studies on the
dynamics of multi-component membranes can be classiﬁed
into two categories; i) dynamics of lateral phase separation
below the phase separation temperature, and ii) dynamics
of concentration ﬂuctuations above the critical temperature. For the details of the domain growth dynamics inthe lower temperatures, which is not the subject of thepresent work, readers are referred to ref. [14]. Experimentally, Honerkamp-Smith et al. have investigated the dynamics of concentration ﬂuctuations in ternary GUVs andshowed that the dynamic critical exponent crosses overfrom a 2D value to a 3D one as the critical temperature isPage 2 of 21

approached from above [15]. The dynamics of concentration ﬂuctuations in membranes was ﬁrst modeled by Sekiet al. [16] and later extended by Inaura and Fujitani [17].
Ramachandran et al. used the general mobility tensor tonumerically calculate the eﬀective diﬀusion coeﬃcient ofconcentration ﬂuctuations [18]. In these theoretical works,
however, they did not take into account the out-of-plane
membrane deformation nor the membrane bilayer structure.
In biomembranes of living cells, the two monolayers have in general diﬀerent compositions, with a uniqueasymmetry between the inner and outer leaﬂets. Furthermore, the two leaﬂets are not independent, but ratherinteract strongly with each other due to various physical and chemical mechanisms [19]. Some experiments
have shown strong positional correlation and domain
registration between domains across the two membrane
leaﬂets [20, 21], while some papers reported the antiregistration of domains in diﬀerent leaﬂets [22–24]. Inspired by the experiments, several simulations have been
performed [25, 26] and some phenomenological models
have been proposed [27, 28] to describe the phase separation in such coupled leaﬂets. Hirose et al. considereda coupled bilayer composed of two modulated monolayers and discussed the static and dynamic properties ofconcentration ﬂuctuations above the transition temperature [29, 30].
The purpose of this paper is to investigate the relaxation dynamics of a binary lipid bilayer membrane in the
one-phase region (rather than phase separated membranes
in a two-phase state). In such membranes, the interplay ofvarious important eﬀects, such as inter-monolayer friction
and composition-deformation coupling, leads to a complex
behavior. In particular, as in usual 3D multi-component
ﬂuids, a chemical potential gradient, and thus mutual diffusion, are induced by the inhomogeneity of the density
diﬀerence between the two lipid species. Such mutual diffusion leads to homogenization of the density diﬀerencein each monolayer, as an irreversible process. In this paper, we take into account the membrane surface tension,
2D hydrodynamics of each monolayer, mutual diﬀusion
in each monolayer, 3D hydrodynamics of the surrounding ﬂuid, and inter-monolayer friction, whereas the ﬂipﬂop motion of the lipid molecules between the two leaﬂets
is not included. In our model, the sources of the energy
dissipation are the viscosities of the monolayers and thesurrounding ﬂuid, the inter-monolayer friction, and themutual diﬀusion in each monolayer. We ﬁnd that the tworelaxation modes associated with the mutual diﬀusion appear in addition to the three previously discussed relaxation modes [5]. These two diﬀusive modes turn out tobe much slower than the other hydrodynamic modes, andbecome even slower in the vicinity of the unstable region
towards the phase separation.
This paper is organized as follows. In sect. 2, we
present the free energy functional of a binary lipid bilayer
membrane by assuming that the membrane deformationand the density deviations from the respective average values are small. The whole set of dynamic equations are

Eur. Phys. J. E (2016) 39: 52

Fig. 1. Schematic representation of a two-component ﬂuid
bilayer membrane consisting of lipid A and lipid B. The hydrophobic chains are arranged in a back-to-back conﬁguration
to form a bilayer. The surface at which the upper and lower
monolayers are in contact with each other is deﬁned as the
mid-surface. The membrane shape is expressed by the heightz = h(x, y) of the mid-surface measured from the z = 0 plane.

introduced in sect. 3, and the surrounding ﬂow ﬁeld is
integrated out to obtain the relaxation equations for themembrane variables. In sect. 4, thermodynamic stabilityof the one-phase state and the wave number dependencies
of various relaxation modes are discussed in detail for both
small and moderate surface tension cases. We also present
our numerical study of the domain relaxation dynamics.
Finally, sect. 5 is devoted for summary and discussion.

2 Free energy
A binary lipid bilayer membrane consists of lipid A andlipid B as schematically presented in ﬁg. 1. In the presence of the surrounding ﬂuid, the hydrophobic tails of lipid
molecules face each other to form a bilayer structure, while
the hydrophilic heads are in contact with the outer ﬂuid.
The surface at which the hydrocarbon tails are in contact
with each other is deﬁned as the mid-surface. Using theheight h(x, y) of the mid-surface from the z = 0 plane in
the 3D Euclidean space, we express the shape of a nearly
ﬂat membrane using the Monge gauge, i.e., z = h(x, y).
Let us write ψJ± the areal mass densities of lipid J(J = A, B) in the upper (+) and lower (−) monolayers,respectively. The total free energy of the bilayer membrane
is generally given by the form√(1)
F = d2 x g ftot (H, ψJ± , ∇⊥ ψJ± ),where ∇⊥ is the gradient operator along the membranesurface, g the determinant of the metric tensor, and d2 x

Eur. Phys. J. E (2016) 39: 52

Page 3 of 21

denotes the integration with respect to x and y. Within
˜ where ∇˜ = (∂x , ∂y )the lowest order in h, we have ∇⊥  ∇is the 2D gradient operator in the projected plane. The
areal free energy density ftot depends on the mean curvature H, the densities and their spacial derivatives. In
general, the free energy also depends on the temperature,
but we shall not write the temperature dependence of any
quantities explicitly. For small membrane deformations˜  1), H and g are approximated as H  (∇˜ 2 h)/2(|∇h|˜ 2 , respectively, within the lowest orderand g  1 + (∇h)
in h.
We assume in this paper that the upper and the lower
monolayers
the sameof lipid molecules,
 number have√√namely, d2 x g ψJ+ = d2 x g ψJ− . We introduce thereference mass densities of the lipid molecules ψJ0 as the
spacial average of the densities for a ﬂat membrane (orprojected mass densities). Then the conservation law for
the lipid molecules is written as√
d2 x g ψJ± = ψJ0 d2 x.
(2)
We further deﬁne the normalized density deviations as
ρ±
J =

ψJ±
− 1.
ψJ0

(3)

With the aid of eq. (3), the conservation law eq. (2) canbe rewritten as


1
√
˜ 2,
d2 x (∇h)
d2 x g ρ±

−
(4)
J
2
up to the second order in h. Notice that the integral on
the left-hand side does not vanish exactly because ψJ0 is
the projected average density.

where ftot and its derivatives ∂ftot /∂ψJ are evaluatedat the reference state. In deriving eqs. (6) and (7), we√
˜ 2 /2, ∂ftot /∂ρ± =have made use of g  1 + (∇h)
J±(∂ftot /∂ψJ )ψJ0 and eq. (4). The right-hand side of eq. (7)can be identiﬁed as the (negative) in-plane pressure for a
ﬂat membrane1 .
The coupling part Fcoup consists of all the possible bilinear couplings between H and ρ±
J . For later convenience,
we introduce the normalized total mass density deviation
 ±

±
ψJ
±
J
J ψJ0 ρJ
−1= 
,
(8)
ρ =
J ψJ0
J ψJ0
and the normalized mass density diﬀerence±
φ± = ρ±A − ρB .

We express Fcoup in terms of bilinear couplings between
˜ 2 h)/2, ρ± and φ± rather than those between H
H  (∇±and ρJ . With this choice of variables, the dynamic equations will be simpliﬁed as we will show in the next section.
Since we have ﬁve independent variables, there should bein principle fourteen coupling parameters in Fcoup 2 . However, we can reduce the number of coupling parameters
by using the invariance of the system under the interchange of the upper and the lower monolayers. For instance, the coupling parameter for (ρ+ )2 should be the˜ 2 h)
same for (ρ− )2 . Also the coupling parameter for ρ+ (∇should have the same magnitude but with an opposite sign
˜ 2 h). Using these symmetric properties, weof that for ρ− (∇are left with eight coupling parameters. Furthermore, it isconvenient to absorb two of them, d and (dimensionless)ν, in the following redeﬁnitions of the variables:
˜ 2 h),
α± ≡ ρ± ± d(∇

2.1 Bilinear free energy
Hereafter we assume that the membrane is weakly deformed and the density deviations are small enough so√that h and ρ±g
J can be treated as small variables. Thenand ftot in eq. (1) can be expanded about the reference
state (h = 0, ψJ± = ψJ0 ) with respect to the small vari˜ ±ables h, ρ±
J , and ∇ρJ . The total free energy is given by
the sum of three contributions
F = Fdef + Fcoup + Fgrad ,(9)

˜ 2 h).
β ± ≡ φ± ± νd(∇

(10)

The two lengths d and νd can be interpreted as the distances between the membrane mid-surface and the two
eﬀective neutral surfaces [5]. Introducing the parameters
k and Λi (i = 1, · · · , 5), we can write Fcoup in the form



k
2
d x
(α )2 + Λ1 (β  )2 + Λ2 α β 
Fcoup =2=+,−
+Λ3 α+ α− + Λ4 β + β − + Λ5 (α+ β − + α− β + ) .

(5)

where Fdef is the deformation part, Fcoup the coupling
part, and Fgrad the gradient part. Each part will be explained in order.
First the deformation part Fdef is given byσ˜ 2 + κ (∇˜ 2 h)2 ,(6)
Fdef = d2 x (∇h)
22Here k has the dimension of areal compression modulus,and Λi are the dimensionless parameters of order unity.
Within the lowest order in the membrane deformations
˜ inand density deviations, we can approximate as ∇⊥  ∇where σ is the membrane surface tension and κ the bending rigidity. The surface tension σ is expressed in terms offtot in eq. (1) as  ∂ftot
ψJ0 ,(7)
σ = ftot −
∂ψJ
=+,−

Let f (ρ1 , · · · , ρk ) denote the free energy density of a kcomponent ﬂuid mixture as a function of the mass densities ρi(i = 1, · · · , k) of each component.
Then the pressure of thisP
mixture is given by p = i ρi (∂f /∂ρi ) − f .
2
Notice that we have already taken into account the term˜(∇2 h)2 in Fdef .

J=A,B

(11)

1

Page 4 of 21

Eur. Phys. J. E (2016) 39: 52

ftot . Then the gradient part Fgrad is given by the sum of
˜ ± and ∇φ
˜ ± . For simplicity, we
the scalar products of ∇ρneglect here the couplings between the diﬀerent leaﬂets
˜ − ). Using again the above symmetric
˜ + ) · (∇φsuch as (∇ρproperties, we have c˜  )2 + λ1 (∇φ˜  )2d2 x
Fgrad =
(∇ρ
2
=+,−

˜  ) · (∇φ˜ ) ,+λ2 (∇ρ(12)
where c has the dimension of energy and is comparable tothermal energy, λ1 and λ2 are the dimensionless parameters of order unity.
Some comments are in order. i) We have thirteen parameters in our free energy; σ, κ, k, d, ν, Λi (i = 1, · · · , 5),c, λ1 and λ2 . In fact they all depend on the temperature
T and the reference densities ψJ0 . In this paper, however,
we regard them as independent parameters although theycannot be varied independently in experiments. In the following sections, we investigate the behaviors of the relaxation rates as these parameters are varied, especiallywhen the instability boundary of the one-phase state isapproached.
ii) In the above total free energy F , terms which are
purely linear in H do not exist. They can be always eliminated by using the invariance of the system under theinterchange of the two leaﬂets, which ﬂips the sign of H.
Notice that the terms which are linear in ρ±
J have alreadybeen taken into account in the deﬁnition of the surface
tension σ in eq. (7).
iii) In principle, the free energy can include terms linear in Gaussian curvature K which is proportional to h2 .
However, without any topological change of the membrane, the integral of K depends only on the geodesic
curvature along the boundary of the membrane. As longas the topology and the geodesic curvature at the edge ofthe membrane are ﬁxed, the integral merely adds a constant to the free energy. For this reason, we do not include
any Gaussian curvature term in eq. (6).

The total free energy is alternatively expressed in termof the Fourier modes asF =d2 q 1 †a Aa + b† Bb ,(2π)2 2(18)

where † denotes the conjugate transpose. In the above, A
and B are symmetric matrices of 3 × 3 and 2 × 2, respectively. Owing to the rotational symmetry, their components depend only on the magnitude of the wave vector,
q = |q̃|, and are given by
A11 = σd2 q 2 + (κ + kd2 Ω0 )d2 q 4 ,

(19)

2

A12 = A21 = −kd Ω1 q ,(20)

A13 = A31 = −kd Ω2 q ,

(21)

2
2

2

A22 = k(2 − Λ3 ) + 2cq ,2(22)

A23 = A32 = k(Λ2 − Λ5 ) + cλ2 q ,(23)

A33 = k(2Λ1 − Λ4 ) + 2cλ1 q ,(24)

B11 = k(2 + Λ3 ) + 2cq 2 ,B12 = B21 = k(Λ2 + Λ5 ) + cλ2 q 2 ,B22 = k(2Λ1 + Λ4 ) + 2cλ1 q 2 .

(25)
(26)
(27)

2

2

and

Here we have introduced the following dimensionless combinations:Ω0 = 2 + 2ν 2 Λ1 + 2νΛ2 − Λ3 − ν 2 Λ4 − 2νΛ5 ,
Ω1 = 2 + νΛ2 − Λ3 − νΛ5 ,
Ω2 = 2νΛ1 + Λ2 − νΛ4 − Λ5 .

(28)
(29)
(30)

It is important to note that a and b are decoupled ineq. (18). This is due to the symmetry of the system under
the interchange of the two monolayers, i.e., a changes its
sign under this interchange while b does not.

2.2 Fourier representation
The in-plane Fourier transform of any function g(x̃) in themonolayer is deﬁned by(13)g(q̃) = d2 x g(x̃)e−iq̃·x̃ ,where x̃ = (x, y) and q̃ = (qx , qy ). It is convenient tointroduce the following new variables:ρ = (ρ+ − ρ− )/2,
φ = (φ+ − φ− )/2,ρ̄ = (ρ+ + ρ− )/2,φ̄ = (φ+ + φ− )/2,ĥ = h/d,
and deﬁne the column vectorsa = (ĥ, ρ, φ)T ,

(14)
(15)
(16)

b = (ρ̄, φ̄)T ,where “T” denotes the transpose.

(17)

3 Dynamic equations
In this section, we present the dynamic equations for a
two-component bilayer membrane surrounded by a viscous
ﬂuid. We shall take into account i) the ﬂows in the surrounding ﬂuid and in the membrane, ii) the frictional force
between the two monolayers, and iii) the mutual diﬀusion
in each monolayer. The surrounding ﬂuid is assumed to
be incompressible, while the membrane itself is compressible [5]. Our dynamic equations are based on the standard
irreversible thermodynamics [31, 32], and ensure that thedissipation in the whole system is non-negative deﬁnite(see appendix A). While our derivation presented in thissection is self-contained, they can be formulated in a more
systematic manner by using the so called Onsager’s variational principle (see appendix B) [33–35].

Eur. Phys. J. E (2016) 39: 52

Page 5 of 21

3.1 Hydrodynamic equations
We use v to denote the velocity ﬁeld of the surrounding
ﬂuid which is assumed to be incompressible and to havea low Reynolds number. Then v for z > 0 and z < 0 obeys
the Stokes equationη∇2 v − ∇p = 0,(31)where η is the shear viscosity, ∇ = (∂x , ∂y , ∂z ) the nabla
operator in 3D space, and p the pressure of the ﬂuid thatis determined by the incompressibility condition∇ · v = 0.

(32)

Let ṽJ± denote the ﬂow velocity of the lipid J in theupper (+) and the lower (−) monolayers. Here the ﬂow
velocity is deﬁned as the lipid mass ﬂux divided by themass density ψJ± . We consider the dynamic equations only
within the linear order in ṽJ± , h and ρ±
J . The average lipid
velocities ṽ ± in the upper and lower monolayers are deﬁned as± ±ψ ± ṽ ± + ψB
ṽB
,
(33)
ṽ ± = A A±±
ψA + ψ B
which can be approximated within the linear order asṽ ± =±±
ψA0 ṽA+ ψB0 ṽB.
ψA0 + ψB0

(34)

(35)

where use has been made of eq. (34) in the second equality.
It should be noted here that the diﬀusive ﬂux of lipid B isgiven by −jd± . Then the continuity equations for the lipids± ±˜A and B, ∂ψJ± /∂t = −∇·(ψJ ṽJ ), can be approximated as∂ρ±A
˜ · ṽ ± − 1 ∇˜ · j±,
= −∇
d
∂t
ψA0
±
∂ρB
˜ · ṽ ± + 1 ∇˜ · j±.
= −∇
d
∂t
ψB0

ij

j

i

On the other hand, the reversible force density due to thein-plane pressure is given by˜ δF ,
˜ δF = −∇f± = −
ψJ0 ∇
(42)
δρ±
δψJ±
J=A,B

up to the linear order [36]. The force balance equations in
the tangential direction of the monolayers are given by
±
±
± Tiz
∓ b(ṽi+ − ṽi− ) = 0,
fi± + ∂j τij

The diﬀusive ﬂux of lipid A is given by±±
jd± = ψA0 (ṽA− ṽ ± ) = −ψB0 (ṽB
− ṽ ± ),where Lφ > 0 is the Onsager coeﬃcient [31], and the sec±−
ond equality follows from the relation, μ± = δF/δψA
±
−1
−1
±
δF/δψB = (ψA0 + ψB0 )δF/δφ . In the deﬁnition of Lφ ,we have intentionally put the factor [(1/ψA0 )+(1/ψB0 )]−1
in order to make eq. (39) simpler. Equation (40) indicates
that, as in usual 3D multi-component ﬂuids, mutual diffusion occurs essentially due to the inhomogeneity of thedensity diﬀerence φ± between the lipid A and B in eachmonolayer. Furthermore, even if φ± are homogeneous, mutual diﬀusion can still be induced by the inhomogeneityof h and ρ± that are coupled to φ± via the free energy.
Next we discuss the force balance conditions. We regard each monolayer as a compressible 2D ﬂuid characterized by the shear viscosity μ and the bulk viscosity ζ.
±in the monolayers are
The 2D viscous stress tensors τijgiven by˜ · ṽ ± .
τ ± = μ(∂i ṽ ± + ∂j ṽ ± ) + (ζ − μ)δij ∇
(41)

(36)

for i = x, y. Here Tij± are the stress tensors of the surrounding ﬂuid Tij = −pδij + η(∂i vj + ∂j vi ) evaluated at
z → ±0. The last term in eq. (43) represents the frictional
forces between the upper and the lower monolayers, andb is the friction coeﬃcient [5, 6].
In the normal z-direction, the restoring force of themembrane is balanced with the normal force due to thesurrounding ﬂuid. Hence we have
+
−
− Tzz
=
Tzz

(37)

We further note that eqs. (36) and (37) can be expressedin simpler forms by using ρ± and φ± as∂ρ±
˜ · ṽ ± ,
= −∇
(38)
∂t
∂φ±
˜ · j±,
= −∇
(39)
φ
∂t
where the diﬀusive ﬂux associated with φ is now deﬁned−1−1 ±as jφ± = (ψA0
+ ψB0
)jd .
As in the standard irreversible thermodynamics, jφ± isassumed to be proportional to the gradient of the eﬀective±
chemical potential μ± = (μ±A /mA ) − (μB /mB ), where mJ±
and μJ are the molecular mass and the chemical potential per molecule for lipid J, respectively [32]. The chem±
ical potentials are given by μ±
J = mJ (δF/δψJ ). Then thediﬀusive ﬂux in eq. (39) becomes 1
1 −1 ˜ ±˜ δF ,
jφ± = −Lφ+(40)
∇μ = −Lφ ∇
ψA0
ψB0
δφ±

(43)

δF
˜ 2 h + (κ + kd2 Ω0 )∇
˜ 2∇
˜ 2h
= −σ ∇δh˜ 2 ρ + Ω2 ∇˜ 2 φ),+kd(Ω1 ∇(44)

where the last expression follows from eqs. (6), (11)and (28)–(30).
We further assume that the non-slip boundary condition holds at the upper and the lower monolayers. Let
v ± = (vx± , vy± , vz± ) denote the velocity of the surrounding
ﬂuid v evaluated at z → ±0. The tangential componentsof v ± should coincide with the average velocities of themonolayers(45)
vi± = ṽi± ,for i = x, y. On the other hand, the normal components vz±
should coincide with the time derivative of the membraneheight h
∂h
.
(46)
vz± =
∂t
Up to now, we have presented a set of dynamic equations given by eqs. (31), (32), (38), (39), (40), (43), (44),(45) and (46) to be solved. As mentioned before, they canbe also derived systematically by using the Onsager’s variational principle explained in appendix B.

Page 6 of 21

Eur. Phys. J. E (2016) 39: 52

3.2 Relaxation equations for membrane variables

Table 1. List of static and dynamic parameters taken fromthe literature [5, 8–11, 37–45] and used in sect. 4.

From the derived dynamic equations, we can integrateout the velocity ﬁelds v and ṽ ± to obtain the relaxation
equations for the spatially Fourier transformed dynamicalvariables, ρ± (q̃, t), φ± (q̃, t) and ĥ(q̃, t) = h(q̃, t)/d (seeeq. (13)). The details are described in appendix C and theresulting equations are∂a= −Γa (q)a,∂t∂b= −Γb (q)b,
∂t

(47)

⎞
c2 B12 q 2
c2 B11 q 2⎠,
Γb (q) = ⎝ 1
1
Lφ B12 q 2 Lφ B22 q 2
2
2

(50)

c1 = [4b + 4ηq + 2(μ + ζ)q 2 ]−1 ,

(51)

⎛

with
2 −1

c2 = [4ηq + 2(μ + ζ)q ]κ2(erg/cm )
−8

10

−4

or 10

η

k
2

(erg)

(erg/cm )

−12

10

70

b

d

c

(cm)

(erg)

−7

10−14

10

μ+ζ

Lφ

(erg · s/cm ) (erg · s/cm ) (erg · s/cm ) (cm /(erg · s))
3

10−2

4

2 × 107

2

10−7

4

1.4 × 10−9(48)

where a and b are deﬁned in eq. (17). In the above, thematrices Γa and Γb are given by
⎛
⎞
A11
A12
A13
⎜ 4ηd2 q
4ηd2 q4ηd2 q ⎟⎜⎟⎜22
Γa (q) = ⎜ c1 A12 q
c1 A22 q
c1 A23 q 2 ⎟⎟ , (49)⎝⎠
1
1
1
Lφ A13 q 2 Lφ A23 q 2 Lφ A33 q 2
2
2
2
and

σ

.

(52)

The eigenvalues of Γa and Γb correspond to the relaxation
rates of the binary bilayer membranes. The equations forthe ﬁve dynamical variables are split into the decoupled
two equations (47) and (48), where eq. (47) changes its
sign under the interchange of the two monolayers, while
eq. (48) does not. This is the consequence of the symmetryof the hydrodynamic equations as well as that of the freeenergy (the latter is discussed after eq. (30)). In the next
section, we will examine how these relaxation rates behaveas the coupling parameters Λi are varied.

4 Results
4.1 Parameter values
In table 1, we list the set of parameter values chosen
in our numerical calculations. Following previous experiments [37–39], the bending modulus κ is set equal to10−12 erg. As discussed after eq. (10), the lengths d anddν are comparable with the monolayer thickness. Then we
may set d = 10−7 cm, with ν of order unity. The combination cλi (i = 1, 2) in Fgrad is related to the line tension
ξ in a phase-separated membrane as ξ ∼ (kB T c)1/2 λi /d.
Since ξ has been measured to be several pN [40], we may
set c = 10−14 erg, with λi of order unity. The surface tension σ can take extremely wide range of values depending

on experimental conditions. For vesicles in a solution, it
can be controlled by changing the osmotic pressure diﬀerence between the inside and outside of the vesicles. In thefollowing, we will examine two cases, namely, the small
tension case with σ = 10−8 erg/cm2 and the moderatetension case with σ = 10−4 erg/cm2 . The coeﬃcients A22
and B11 in eqs. (22) and (25) can be interpreted as themoduli associated with the total densities in the upper
and the lower monolayers, respectively (to be more precise, the moduli of their linear combinations ρ and ρ̄ ineq. (14)). Then k(2 − Λ3 ) and k(2 + Λ3 ) are comparablewith the areal compression moduli. Following previous experiments [41, 42], we set k = 70 erg/cm2 with |Λ3 |  1.
The remaining parameters that have yet to be determined in the free energy are Λ1 , Λ2 , Λ4 and Λ5 . These
parameters depend on the temperature and the averagecomposition. We ﬁnd in the following, however, that the
behavior of the decay rates is not sensitive to these parameters, unless the reduced temperatures τa and τb deﬁned
below in eqs. (71) and (84) are very close to zero. Whenthey are close to zero (but positive), the associated diffusive modes become extremely slow. This point will be
discussed later in more detail.
We next discuss the kinetic parameters. The membrane viscosities μ and ζ appear only as a sum μ + ζ in
Γa and Γb . Since we could not ﬁnd any reliable value ofζ in the literature, we set μ + ζ = 10−7 erg · s/cm2 . (Themembrane bulk viscosity ζ was neglected in ref. [5].) The
Onsager coeﬃcient Lφ for the mutual diﬀusion is roughly
estimated as follows. Assuming that the mutual diﬀusion
constant is on the same order as the self-diﬀusion constant
D of a lipid molecule, we have D ∼ (Γa )33 /q 2 ∼ Lφ k(see eqs. (24), (48) and (49)). Using the value D ∼10−7 cm2 /s [43], we obtain Lφ = 1.4 × 10−9 cm4 /(erg · s).
Several authors have reported diﬀerent values of the friction coeﬃcient b [8–11, 44, 45]. Since they are in the rangeof 107 –3 × 108 erg s/cm4 , we set b = 2 × 107 erg · s/cm4 inthis paper.

4.2 Stability conditions
The wave number dependent susceptibilities χa (q) andχb (q) are deﬁned as the reciprocals of the eigenvalues ofthe matrices A and B in eq. (18), respectively. Then thethermodynamic stability of the one-phase state (withoutEur. Phys. J. E (2016) 39: 52

0.3

Page 7 of 21

necessary for stability

0.2
0.1
0

unstable

-0.1
-0.4

-0.2

0

0.2

0.4

Fig. 2. Stability diagram of the one-phase state at q = 0 in the(Λ1 , Λ4 )-plane as expressed by eq. (55). For the thermodynamic
stability of the one-phase state, Λ1 and Λ4 need to be in thegray region. The cross corresponds to the parameter values
Λ1 = Λ4 = 1/16 which we use in ﬁg. 3 and in the numericalanalysis in ﬁgs. 5–10.

any phase separation) is ensured when these susceptibilities are all positive. Since A and B are 3 × 3 and 2 × 2
matrices, respectively, there are three χa and two χb values which can be explicitly obtained in principle. However,since their full expressions are tedious, we discuss here
the conditions for the thermodynamic stability at q = 0
and ∞. More detailed discussions are given in appendix D
where we also show that the instability characterized byintermediate wave numbers does not occur as long as thestability conditions at q = 0 and ∞ are satisﬁed.
As q → ∞, we ﬁnd that the susceptibilities χa and χb
are both positive if and only if
0 < λ1 −λ22≡ Δλ .
4(53)

Hereafter we assume that the above condition is always
satisﬁed. The stability at q = 0, on the other hand, isensured by the positivity of χa (0) and χb (0), which is realized if and only if|Λ3 | < 2,
0 < Λ1 ,|Λ4 | < 2Λ1 ,

(54)
(55)

Fig. 3. (a) Stability diagram of the one-phase state at q = 0
in the (Λ2 , Λ3 , Λ5 )-space when Λ1 = Λ4 = 1/16 as marked inﬁg. 2. The one-phase state is stable if Λ2 , Λ3 and Λ5 are inside
the region enclosed by the surface. (b) The enclosing surface
consists of two surfaces; the blue one at which 1/χa (0) = 0
and the red one at which 1/χb (0) = 0. The cross sections of
the stable region on the (Λ2 , Λ5 )-plane when (c) Λ3 = −3/2
and (d) Λ3 = 15/8. The circles marked with (A)–(C) in (d)correspond to the parameter values used in ﬁgs. 5–13.

and|Λ2 − Λ5 | < [(2 − Λ3 )(2Λ1 − Λ4 )]1/2 ,(56)|Λ2 + Λ5 | < [(2 + Λ3 )(2Λ1 + Λ4 )]1/2 .

(57)

The conditions eqs. (54) and (55) are equivalent to A22 ,A33 , B11 , B22 > 0 at q = 0. In ﬁg. 2, we plot the condition eq. (55) on the (Λ1 , Λ4 )-plane. For the stability of the
one-phase state, Λ1 and Λ4 need to be within the gray region. Otherwise the system is unstable towards the phaseseparation. Given that Λ1 and Λ4 are ﬁxed at the values
satisfying eq. (55), the stability conditions for Λ2 , Λ3 andΛ5 are given by eqs. (54), (56) and (57).
In ﬁg. 3(a), we show the stable region in the (Λ2 , Λ3 ,Λ5 )-space when Λ1 = Λ4 = 1/16 as marked by a cross inﬁg. 2. The stable region is enclosed by a surface which consists of blue and red parts. On the blue surface (ﬁg. 3(b)left), one of the three χa -values diverges at q = 0 and its

corresponding mode becomes unstable, whereas on the red
surface (ﬁg. 3(b) right), one of the two χb -values diverges
at q = 0. One can conﬁrm from eqs. (56) and (57) that thecross section of the stable region on the (Λ2 , Λ5 )-plane atconstant Λ3 is given by an oblique rectangle whose center is at Λ2 = Λ5 = 0. In ﬁgs. 3(c) and (d), we present
the cross sections at Λ3 = −3/2 and 15/8, respectively.
As shown in (c), the apex coordinates of the rectangle Λ±
are given by
1 (2 + Λ3 )(2Λ1 + Λ4 ) ± (2 − Λ3 )(2Λ1 − Λ4 ) .
2(58)These values are Λ± ≈ 0.153 ± 0.234 and Λ± ≈ 0.426 ±0.0442 in ﬁg. 3(c) and (d), respectively.
As Λ3 is increased towards 2, the blue sides of thecross section become longer while the red sides become
Λ± =Page 8 of 21

Eur. Phys. J. E (2016) 39: 52

eigenvalues of Γa and Γb in eqs. (47) and (48), respectively. In the following calculations, we set the parametervalues to Λ1 = Λ4 = 1/16 and Λ3 = 15/8 as in ﬁg. 3(d),while (Λ2 , Λ5 ) are varied. For simplicity, we further set
λ1 = λ2 = ν = 1 which also satisfy eq. (53).
(a) undulated and anti-registered(b) flat and registered4.3.1 Eigenmodes of Γa : small tension caseunstable

unstableFig. 4. Schematic illustrations of the two possible instabilities. Note that for the sake of clarity (i) the diﬀerence in thelipid heights is not drawn, and (ii) a strong phase separationis represented. In this paper we investigate, however, only thehomogeneous phase in the vicinity of the phase separation. (a)Anti-registered instability. As one of the χa values diverges, thedensity diﬀerence between the two monolayers and the bending mode become unstable. (b) Registered instability. As oneof the χb values diverges, the sum of the densities in the upperand lower monolayers becomes unstable.

shorter. In the limit of Λ32, the stable region eventually turns out to be a line
√ segment whose√ endpoints
,
Λ
)
=
(±
2Λ
+
Λ
,
±
2Λ1 + Λ4 ) =are√given by
(Λ
2
5
1
4
√
−2, on the other(± 3/4, ± 3/4). In the limit of Λ3
hand, the stable region shrinks
with the√ to a line segment√
endpoints at (Λ2 , Λ5 ) = (± 2Λ1 − Λ4 , ∓ 2Λ1 − Λ4 ) =(±1/4, ∓1/4).
Even if we choose other Λ1 and Λ4 values in the stableregion of ﬁg. 2, the qualitative features of the stable region in the (Λ2 , Λ3 , Λ5 )-space remains the same. However,as the combination (Λ1 , Λ4 ) approaches the boundary of
the gray region, the stable region in the (Λ2 , Λ3 , Λ5 )-spacebecomes narrower, and eventually disappears just at theboundary of the stable region. In ﬁg. 4, we illustrate thecorresponding instabilities to take place. When one of theχa -values diverges, a certain linear combination of ρ, φand ĥ becomes unstable as in ﬁg. 4(a), while a linear combination of ρ̄ and φ̄ becomes unstable as in ﬁg. 4(b) when
one of the χb -values diverges. Hereafter we shall call the
instabilities of type (a) and (b) the “anti-registered instability” and the “registered instability”, respectively. Notethat these two types of instabilities are purely the consequences of the symmetry of the system (see also thesentences after eq. (30)).
So far, we have discussed the stability conditions of the
one-phase state at q = 0 and ∞. In appendix D, we showthat the instability does not take place for intermediatewave numbers as long as the membrane is stable at q = 0
and ∞. Hence the stability conditions are generally given
by eqs. (53)–(57).
4.3 Relaxation rates
As a main result of this paper, we next examine the relaxation rates (or decay rates) of the various hydrodynamic
modes in the one-phase state. They are obtained from theSetting the parameters as (Λ2 , Λ5 ) = (0.193, 0.233), weplot in ﬁg. 5(a) the three eigenvalues of Γa denoted by γai(i = 1, 2, 3 and γa1 > γa2 > γa3 ), and in (b) the three
diagonal elements of Γa denoted by (Γa )ii (i = 1, 2, 3 and(Γa )11 > (Γa )22 > (Γa )33 ) for a small surface tension,
σ = 10−8 erg/cm2 . Similar plots are given in ﬁg. 6 when(Λ2 , Λ5 ) = (−0.023, 0.063) with the same surface tension
value. These choices of the parameters are marked with(A) and (B) in ﬁg. 3(d). The system is far from and close
to the unstable region in ﬁgs. 5 and 6, respectively. In thelatter case, at least one of the eigenvalues of A becomesvery small, and the anti-registered instability shown inﬁg. 4(a) is about to take place.
In both ﬁgs. 5 and 6, the fastest decay rate γa1 is foundto be(Γa )22 (q  qmc ),
(59)
γa1 
(Γa )11 (q  qmc ).
Here the mode crossing wave number is given byqmc =ηk(2 − Λ3 )
,
κb

(60)

at which (Γa )11 = (Γa )22 holds. The positivity of qmc >0 follows from the stability condition eq. (54). Using
the present parameter values, we obtain qmc = 4.38 ×
103 cm−1 .
Let us introduce the “quasi-equilibrium” state of ρ forgiven ĥ and φ asρ(2)e (ĥ, φ, q) = −A12 ĥ + A23 φ
.
A22

(61)
(2)

We use the term “quasi-equilibrium” because ρe (ĥ, φ, q)
minimizes the free energy F under the condition that theother variables are ﬁxed at (ĥ, φ). It can be obtained byequating the second row of eq. (47) to zero, and solvingfor ρ. Then we can rewrite the second row of eq. (47) as(2)∂ρ/∂t = −(Γa )22 (ρ − ρe ). Hence γa1  (Γa )22 for q 
qmc indicates that ρ relaxes towards the quasi-equilibrium(2)
state ρe with the decay rate γa1 , while the other variablesĥ and φ are almost unchanged (frozen) during this process.
In other words, ρ relaxes much faster than ĥ and φ. In
this regime, we can approximate A22  k(2 − Λ3 ) andc1  1/(4b) in eqs. (22) and (51), respectively, and the
decay rate scales as γa1  (Γa )22  k(2−Λ3 )q 2 /(4b) ∼ q 2 .
Similarly, the decay rate γa1 for q  qmc corresponds
to the relaxation of ĥ to its quasi-equilibrium stateĥ(2)e (ρ, φ, q) = −A12 ρ + A13 φ
,
A11

(62)

Eur. Phys. J. E (2016) 39: 52

10 -6

Page 9 of 21

(a)

10 -6

10 -8
10 -10

10 -6

10 -8

10

103

105

10 -10

107

(b)

10 -6

10 -8
10 -10

(a)

10

103

105

107

103

105

107

(b)

10 -8

10

103

105

10 -10107

Fig. 5. Eigenmodes of Γa for small tension case; σ = 10−8
erg/cm2 . The parameter values are (Λ2 , Λ5 ) = (0.193, 0.233) as
marked (A) in ﬁg. 3(d), and the membrane is not close to theanti-registered instability boundary. (a) Plots of the relaxation
rates γai (i = 1, 2, 3) as a function of the wave number q. (b)Plots of the diagonal elements (Γa )ii (i = 1, 2, 3) of the matrix
Γa as a function of the wave number q (dashed color lines).
The eﬀective decay rate γφ∗ is plotted with a black dashed line.
For comparison, the relaxation rates γai in (a) are also plottedwith grey solid lines. For convenience, all the plotted quantities
are divided by q 2 .

while both ρ and φ are frozen during the relaxation of ĥ.
For q  q ∗ with

σ
∗
,
(63)
q =
κ
one can approximate eq. (19) as A11  (κ + kd2 Ω0 )d2 q 4 .
For the parameter values used in ﬁgs. 5 and 6, we haveqmc  q ∗ because q ∗ = 100 cm−1 (see also the sentences
below eq. (76)). Then the decay rate scales as γa1 (Γa )11  (κ + kd2 Ω0 )q 3 /(4η) ∼ q 3 for q  qmc .
The second fastest decay rate γa2 behaves as⎧(Γa )12 (Γa )21⎪
 (Γa )11⎪⎨(Γa )11 −(Γa )22γa2 ⎪⎪⎩(Γa )22 − (Γa )12 (Γa )21  (Γa )22(Γa )11(q  qmc ),(q  qmc ).
(64)

10

Fig. 6. Eigenmodes of Γa for small tension case; σ = 10−8
erg/cm2 . The parameter values are (Λ2 , Λ5 ) = (−0.023, 0.063)
as marked (B) in ﬁg. 3(d), and the membrane is close to theanti-registered instability. (a) Plots of the relaxation rates γai
as a function of the wave number q. (b) Plots of the diagonal elements (Γa )ii of the matrix Γa as a function of the wave
number q (dashed color lines). The eﬀective decay rate γφ∗ isplotted with a black dashed line. For comparison, the relaxation rates γai in (a) are also plotted with grey solid lines. For
convenience, all the plotted quantities are divided by q 2 .

Let us introduce the quasi-equilibrium states of ĥ and ρ
for given φ asĥ(1)e (φ, q) =A12 A23 − A13 A22φ,
A11 A22 − A212

(65)

ρ(1)
e (φ, q) =A12 A13 − A11 A23
φ,
A11 A22 − A212

(66)

which minimize the total free energy F under the condition that φ is ﬁxed. They are obtained by equating theﬁrst and the second rows of eq. (47) to zero, and solving
simultaneously for ĥ and ρ. Assuming that the relaxation(2)of ρ is much faster than that of ĥ, we substitute ρ  ρegiven by eq. (61) into the ﬁrst row of eq. (47) to obtain(Γa )12 (Γa )21  ∂ ĥ − (Γa )11 −ĥ − ĥ(1)e∂t(Γa )22
 −(Γa )11 ĥ − ĥ(1)
.
e

(67)

Hence the decay rate γa2 for q  qmc in eq. (64) corresponds to the relaxation of ĥ to the quasi-equilibrium

Page 10 of 21

Eur. Phys. J. E (2016) 39: 52

(1)

(2)

state ĥe , while φ is frozen and ρ instantly decays to ρe .
In this regime, we have γa2  (Γa )11  σq/(4η) ∼ q forq  q ∗ , and γa2  (Γa )11  (κ + kd2 Ω0 )q 3 /(4η) ∼ q 3 forq ∗  q  qmc . For q  qmc , on the other hand, γa2 is(1)
associated with the relaxation of ρ towards ρe , while φ(2)
is frozen and ĥ instantly decays to ĥe . In this regime, wehave γa2  (Γa )22  k(2 − Λ3 )q 2 /(4b) ∼ q 2 .
From eqs. (59) and (64), we see that the mode crossing
occurs around q  qmc ; the fastest mode is associated with
ρ for q < qmc while it is dominated by ĥ for q > qmc . Sucha mode crossing behavior between the density and thecurvature was predicted by Seifert and Langer for singlecomponent lipid bilayer membranes without any surface
tension [5]. In table 2(a), we present a list of the approximate expressions for γa1 and γa2 when the membrane tension is small (the threshold tension σt in the table caption
is deﬁned in eq. (76) below).
We now discuss the slowest decay rate γa3 . Assuming(1)ĥ and ρ vary much faster than φ, we substitute ĥ  ĥe(1)and ρ  ρe into the third row of eq. (47) to obtain∂φ −γφ∗ φ.
∂t

Table 2. Approximate expressions for the decay rates. (a) The
two fastest decay rates γa1 and γa2 associated with Γa for thesmall tension case σ < σt . (b) The two fastest decay rates γa1
and γa2 associated with Γa for the moderate (larger) tension
case σ > σt . The threshold tension σt is deﬁned in eq. (76).
(c) The slowest decay rate γa3 associated with Γa for both thesmall and the moderate tension cases. (d) The slowest decayrate γb2 associated with Γb for both the small and the moderatetension cases. The characteristic wave numbers q ∗ , qmc , qa andqb are deﬁned in eqs. (63), (60), (72) and (85) respectively.
(a)

(b)

q  q∗q∗  qγa1

σq4η

(κ + kd2 Ω0 )q 34ηγa3k(2 − Λ2 )q 2
4bq  qa1
Lφ kτa q 22

qa  qLφ cΔλ q 4

(69)
(d)

In the small and large wave number limits, its asymptotic
behaviors areLφ kτa q 2 /2 ∼ q 2 (q → 0),
∗
γφ →
(70)
Lφ cΔλ q 4 ∼ q 4 (q → ∞),where the reduced temperature τa is deﬁned by3(71)

and Δλ was deﬁned before in eq. (53). When the stability
conditions in eqs. (54)–(56) are satisﬁed, one can showthat τa is positive. As the unstable region is approached,τa becomes smaller and eventually vanishes just at theboundary. Then the anti-registered instability in ﬁg. 4(a)takes place at the boundary as well as in the unstableregion.
The crossover wave number qa between the two limits
in eq. (70) is given by

kτa
.
(72)
qa =
2cΔλ
3

k(2 − Λ2 )q 2
4b(κ + kd2 Ω0 )q 34η(c)(Λ2 − Λ5 )2,
2 − Λ3(κ + kd2 Ω0 )q 3
4η

σq
4η

With the aid of eqs. (65) and (66), the eﬀective decay rate
γφ∗ in the above equation can be obtained asqmc  q

2

γa2

γa2

τa = 2Λ1 − Λ4 −q ∗  q  qmc
k(2 − Λ2 )q
4b

γa1

(68)

Lφ (det A)q 2.
γφ∗ =2(A11 A22 − A212 )q  q∗

Because the parameters Λi depend on the temperature, one
may change the temperature in order to enter the unstable region at some temperature Tc . Since the combination τa (τb )changes sign at the transition, it is proportional to T − Tc .
Therefore, being dimensionless, it is a reduced temperature.
See also eqs. (D.11) and (D.12) in appendix D where the (normalized) susceptibilities for φ and φ̄ at q = 0 are given by τa−1and τb−1 , respectively.

γb2q  qb1
Lφ kτb q 22

qb  qLφ cΔλ q 4

In ﬁgs. 5(b) and 6(b), we have also plotted γφ∗ . We see thatγφ∗ provides a perfect ﬁt to the slowest mode γa3 . Thus
γa3 corresponds to the relaxation rate of φ, while ĥ and ρ(1)(1)
instantly change to their equilibrium values ĥe and ρe ,respectively. In table 2(c), the approximate expressions forthe slowest rate γa3 are summarized.
In ﬁg. 5(b), we see that the bare rate (Γa )33 almost
coincides with the eﬀective rate γφ∗  γ3a . This canbe understood as follows. For the parameters used inﬁg. 5, the reduced temperature is approximately given byτa  2Λ1 − Λ4 , while we have A33  k(2Λ1 − Λ4 ) when themembrane is far from the unstable region (see eq. (24)).
We thus have γφ∗  Lφ A33 q 2 /2 = (Γa )33 . The crossoverwave number given by eq. (72) is qa = 1.52 × 107 cm−1
which is microscopic and may not be measurable in experiments. On the other hand, the parameters used in ﬁg. 6yield τa = 3.33 × 10−3 and qa = 3.94 × 106 cm−1 . Hencethe crossover from γa3 ∼ q 2 to ∼ q 4 is measurable as inusual near critical ﬂuids [46]. Note that the q 4 -dependence
in large wave numbers is not due to the coupling with theother modes, but is just a consequence of diﬀusion when
there are squared-gradient terms in the free energy as ineq. (12).

Eur. Phys. J. E (2016) 39: 52

10 -6

Page 11 of 21

(a)

10 -6
10 -8

10 -8
10 -10

103

10

105

10 -10

107

(b)
10 -6

10 -6

10 -8

10 -8

10 -10

103

10

105

10 -10107

Fig. 7. Eigenmodes of Γa for moderate tension case; σ = 10−4
erg/cm2 . The parameter values are (Λ2 , Λ5 ) = (0.193, 0.233) as
marked (A) in ﬁg. 3(d), and the membrane is not close to theanti-registered instability boundary. (a) Plots of the relaxation
rates γai as a function of the wave number q. (b) Plots ofthe diagonal elements (Γa )ii of the matrix Γa as a function of
the wave number q (dashed color lines). The eﬀective decayrate γφ∗ is plotted with a black dashed line. For comparison,
the relaxation rates γai in (a) are also plotted with grey solidlines.

4.3.2 Eigenmodes of Γa : moderate tension caseIn ﬁgs. 7 and 8, we show (a) the relaxation rates γai and(b) the diagonal elements (Γa )ii of Γa for a moderate surface tension, σ = 10−4 erg/cm2 . In these plots, all theparameters except σ are the same as in ﬁgs. 5 and 6. In
the whole wave number range, the decay rates can be approximated asγa1  (Γa )11 ,γa2  (Γa )22 −(73)(Γa )12 (Γa )21
 (Γa )22 ,(Γa )11

γa3  γφ∗ ,

(74)
(75)

where γφ∗ was deﬁned in eq. (69). The fastest decay rate(2)

γa1 is associated with the relaxation of ĥ to ĥe , while ρand φ are frozen. The second decay rate γa2 corresponds to(1)
the relaxation of ρ to ρe , while φ is frozen and ĥ instantly(1)changes to ĥe . The slowest decay mode φ relaxes by theeﬀective decay rate γφ∗ , while ĥ and ρ instantly change to
(1)

ĥe

(1)

(a)

and ρe , respectively.

10

103

105

107

103

105

107

(b)

10

Fig. 8. Eigenmodes of Γa for moderate tension case; σ = 10−4
erg/cm2 . The parameter values are (Λ2 , Λ5 ) = (−0.023, 0.063)
as marked (B) in ﬁg. 3(d), and the membrane is close to theanti-registered instability. (a) Plots of the relaxation rates γai
as a function of the wave number q. (b) Plots of the diagonal
elements (Γa )ii of the matrix Γa as a function of the wave number q (dashed color lines). The eﬀective decay rate γφ∗ is plottedwith a black dashed line. For comparison, the relaxation rates
γai in (a) are also plotted with grey solid lines.

The slowest decay rate γ3a  γφ∗ in ﬁgs. 7 and 8 isalmost the same as in ﬁgs. 5 and 6 for which the membrane tension is very small (see table 2(c)). However, unlike in ﬁgs. 5 and 6, the mode crossing behavior between
the two fast (bending and density) modes does not occur
for the moderate tension case. Recently, the absence of themode crossing behavior due to the membrane tension has
been reported in the experiment [47], and theoretically discussed for single-component lipid bilayer membranes [48].
Since the minimum of (Γa )11 /q 2 is located around q ∼q ∗ and (Γa )22 /q 2 is almost constant, the condition qmc 
q ∗ gives the threshold surface tension2
1 kη(2 − Λ3 )
,
(76)
σt 
κ
b
below which the mode crossing occurs. For Λ3 = 15/8
and other parameter values, we can estimate σt ≈ 1.91 ×
10−5 erg/cm2 . Table 2(a) and (b) summarize the approximate expressions of the two fastest rates of Γa for small
tension case (σ < σt ) and for large tension case (σ > σt ),respectively. Notice that in the small tension case, we always have q ∗ < qmc .

Page 12 of 21

Eur. Phys. J. E (2016) 39: 52

102

102

(a)

10-2

10-2

10-6

10-6

10-10

10-10

10

103

105

102

107

10

10-2

10-2

10-6

10-6

10-10

10-10

103

105

107

4.3.3 Eigenmodes of Γb

1Tr Γb ± {Tr Γb + 2(Γb )−1
11 det Γb } ,
2

(77)
(78)

where eq. (78) follows from (Γb )11  (Γb )22 and (Γb )211 (Γb )12 (Γb )21 . Then we obtain approximatelyγb1  (Γb )11 ,Lφ q 2 det Bγb2 ≡ γφ̄∗ .
2B11

107

103

105

107

Fig. 10. Eigenmodes of Γb when the parameter values are(Λ2 , Λ5 ) = (0.445, 0.405) as marked (C) in ﬁg. 3(d), and themembrane is close to the registered instability boundary. (a)Plots of the relaxation rates γbi as a function of the wave number q. (b) Plots of the diagonal elements (Γb )ii of the matrix
Γb as a function of the wave number q (dashed color lines).
The eﬀective decay rate γφ̄∗ is plotted with a black dashed line.
For comparison, the relaxation rates γbi in (a) are also plottedwith grey solid lines.

Equating the right-hand side of eq. (48) to zero, we
obtain the quasi-equilibrium variables as

In ﬁgs. 9 and 10, we plot the eigenvalues and diagonal
elements of Γb in eq. (50). The parameters are chosen as(Λ2 , Λ5 ) = (0.193, 0.233) in ﬁg. 9 and (0.445, 0.405) inﬁg. 10. These two choices are marked with (A) and (C) in
ﬁg. 3(d). Since Γb is a 2 × 2 matrix, its eigenvalues γb can
be easily obtained as



105

(b)

10

Fig. 9. Eigenmodes of Γb when the parameter values are(Λ2 , Λ5 ) = (0.193, 0.233) as marked (A) in ﬁg. 3(d), and themembrane is not close to the anti-registered instability boundary. (a) Plots of the relaxation rates γbi (i = 1, 2) as a functionof the wave number q. (b) Plots of the diagonal elements (Γb )ii(i = 1, 2) of the matrix Γb as a function of the wave number
q (dashed color lines). The eﬀective decay rate γφ̄∗ is plottedwith a black dashed line. For comparison, the relaxation rates
γbi in (a) are also plotted with grey solid lines.

1(Γb )11 + (Γb )22
γb =
2


± {(Γb )11 − (Γb )22 }2 + 4(Γb )12 (Γb )21

103

102

(b)

10

(a)

(79)
(80)

B12
φ̄,
B11
B21
ρ̄.
φ̄e (ρ̄, q) = −B22ρ̄e (φ̄, q) = −

(81)
(82)

As in the previous subsections, the fastest decay rate γb1 (Γb )11 is associated with the relaxation of ρ̄ to ρ̄e while φ̄ isfrozen. However, as shown in ﬁgs. 9 and 10, the decay rate
γb1 is very large, and our theory, in which inertial eﬀect isneglected, may not properly describe the dynamics of sucha small time scale [5]. Hence we do not further discuss thewave number dependence of γb1 .
Nevertheless, we can discuss the slower relaxation ofφ̄ because ρ̄ relaxes rapidly to the quasi-equilibrium valueρ̄e (φ̄). With the aid of eq. (81), substitution of ρ̄  ρ̄e (φ̄)
into the second row of eq. (48) yields ∂ φ̄/∂t  −γφ̄∗ φ̄.
From eq. (80), we see that the slower decay rate γb2 
γφ̄∗ corresponds to the relaxation of φ̄, while ρ̄ instantly
changes to ρ̄e (φ̄, q). In the small and large wave number

Eur. Phys. J. E (2016) 39: 52

Page 13 of 21

limits, the asymptotic behaviors areLφ kτb q 2 /2 ∼ q 2 (q → 0),∗γφ̄ →(q → ∞),
Lφ cΔλ q 4 ∼ q 4

1.6
1.2

(83)

0.4

where the other reduced temperature τb is deﬁned by (seefootnote 3 on page 10)
τb = 2Λ1 + Λ4 −(Λ2 + Λ5 )2.
2 + Λ3

0.8

0
-0.4
10 5

(84)

When the stability conditions in eqs. (54), (55) and (57)are satisﬁed, τb is positive in the stable region. As the unstable region is approached, τb becomes smaller and eventually vanishes at the boundary where the registered instability in ﬁg. 4(b) starts to take place. The crossover wave
number qb between the two limits in eq. (83) is given by

kτb
.
(85)
qb =
2cΔλ
When the system is away from the unstable region (τb =0.141) as in ﬁg. 9, we have qb = 2.57 × 107 cm−1 whichis too large to be observed. However, when the system
is close to the unstable region (τb = 1.10 × 10−3 ) as inﬁg. 10, we have qb = 2.27 × 106 cm−1 which is measurable
in experiments. In table 2(d), the approximate expressions
for the slowest rate γb2 are summarized.

10 4 3
10
10 2

10

1

10 -1 0

1

2

3

4

5

6

4000
2000
0
-2000
-4000
10 5 410
10 3 2
1010

1

10 -1 0

1

5

6

3

4

2

6

4

5

2

3

1.2
0.8
0.4

4.4 Domain relaxation dynamics
In this subsection, we examine the relaxation dynamics ofa domain in which φ is larger than the outside. When thesystem is in the stable region as in eqs. (54) and (55), sucha domain should relax to a homogeneous state (ĥ, ρ, φ) =0. Let us assume that the initial state (ĥ0 , ρ0 , φ0 ) at t = 0
is described by one-dimensional proﬁlesĥ0 (x) = ρ0 (x) = 0,Δφ
Ld − |2x − L|
φc + tanh,
φ0 (x) =
2
2

(86)
(87)

while these proﬁles are homogeneous in y-direction. The
proﬁle φ0 (x) represents a patch centered at x = L/2, andits size and interfacial thickness are given by Ld and ,
respectively. The diﬀerence of φ between the inside andthe outside the initial domain is given by Δφ, whereas φcis determined so that the spacial average of φ0 vanishes.
We will not discuss the other variables (ρ̄, φ̄) because they
are not coupled to (ĥ, ρ, φ).
The three variables can be generally expressed as
Fourier series deﬁned by
"
!
∞

2πinx
g(qn , t),g(x, t) =
exp
(88)
L
n=−∞
where
qn =

2πn
.
L

(89)

0
-0.4
10 5 410
10 3 2
1010

1

10 -1 01

Fig. 11. Time evolutions of φ(x, t)/Δφ, ĥ(x)/Δφ and ρ(x, t)
Δφ from the initial state given by eqs. (86) and (87). The
parameters are the same as in ﬁg. 6, i.e., σ = 10−8 erg/cm−2and (Λ2 , Λ5 ) = (−0.023, 0.063).

Let a0 (qn ) = (0, 0, φ0 (qn )) denote the Fourier modes of
the initial state given by eqs. (86) and (87). Since the time
evolution of each Fourier mode is governed by eq. (47), we
can write a(qn , t) = e−Γa (q)t a0 (qn ) with q = |qn |. The
matrix Γa (q) can be diagonalized by using its eigenvalues
γi (q) and their respective eigenvectors wi (q) asΛa (q) = W −1 (q)Γa (q)W (q),(90)

where Λa = diag(γa1 , γa2 , γa3 ) is the diagonalized matrix
and W = (w1 , w2 , w3 ). Then a(x, t) can be generally written as
"
!
nc

2πinx
W e−Λa (q)t W −1 a0 (qn ),a(x, t) =
exp
L
n=−n
c

(91)

Page 14 of 21

0.3

Eur. Phys. J. E (2016) 39: 52

10 3

(a)

(a)

0.2

0.1
200
1

600

10 2
0

1000 1400

0.1

0.2 0.3 0.4 0.5
(b)

10-1

(b)

10 -1
10-2

10 -2
10 -3
0

10-3
0

10

20

30

40

20

30

40

50

Fig. 12. Semilogarithmic plots of the time evolutions of|ĥ(qn , t)|/Δφ, |ρ(qn , t)|/Δφ and |φ(qn , t)|/Δφ for (a) large tand (b) small t. Here we choose qn = 3.14 μm−1 . In (a), we
also plot |φ0 (qn )|e−γa3 t with a dashed line.

where we have introduced a cut-oﬀ wave number set by
the monolayer thickness d
π
2πnc
= .
L
d

10

(92)

In ﬁg. 11, we present the time evolution of φ(x, t),ĥ(x, t) and ρ(x, t) obtained from eq. (91) by setting L =6000,  = 10 and Ld = 1000 in nm. The other parameters
are the same as in ﬁg. 6 and the system is close to theanti-registered instability. Notice that a(x, t) divided by
Δφ is independent of Δφ since eq. (47) is linear in a(x, t).
For t ≤ 103 ms, |ĥ| and |ρ| increase while φ remains almost
the same. This means that, within a small time interval,
non-zero φ induces the bending ĥ and the density diﬀerence ρ which were initially both zero. For t ≥ 103 ms, all
the three variables become smaller and almost vanish fort ≥ 105 ms.
The above dynamics can be roughly understood bylooking at the time evolution of the Fourier modes at
qn  2π/Ld . In ﬁgs. 12(a) and (b), the time evolutions of |ĥ(qn )|, |ρ(qn )| and |φ(qn )| are presented at qn =3.14 μm−1 for which the decay rates are γa1 = 785 s−1 ,(1)
γa2 = 110 s−1 and γa3 = 0.157 s−1 . In ﬁg. 13, |ρ −ρe (φ0 )|(2)and |ĥ − ĥe (ρ0 , φ0 )| are plotted as a function of t for the

Fig. 13. Semilogarithmic plots of the time evolutions of (a)(2)|ĥ(qn , t) − ĥe (ρ0 , φ0 , qn )| (solid line) and e−γa1 t (dashed line),(1)and (b) |ρ(qn , t) − ρe (φ0 , qn )| (solid line) and e−γa2 t (dashedline).

same wave number as in ﬁg. 12. As for the long time be−1−1( γa1), ﬁg. 12(a) shows that all thehavior, t  γa2
three variables decay exponentially with a common de(1)cay rate γa3 . In this regime, we have ĥ(qn )  ĥe (φ, qn ),(1)
ρ(qn )  ρe (φ, qn ) and∗φ(qn )  φ0 (qn )e−γa3 t  φ0 (qn )e−γφ t ,(93)

as in the discussion after eq. (72). Substituting eq. (93)
into eqs. (65) and (66), we then have
A12 A23 − A13 A22
φ0 e−γa3 t ,A11 A22 − A212
A12 A13 − A11 A23
φ0 e−γa3 t ,ρ(qn ) 
A11 A22 − A212ĥ(qn ) 

(94)
(95)

which decay exponentially with the common rate γa3 .
For shorter times, on the other hand, |ĥ| and |ρ|rapidly vary while φ stays almost constant, as shown inﬁg. 12(b). In ﬁgs. 12 and 13, the chosen qn = 3.14 μm−1 ismuch larger than the mode crossing wave number qmc =0.438 μm−1 . Then the fastest decay rate γa1 corresponds(2)
to the relaxation of ĥ to ĥe while ρ and φ are frozen.
Notice that in ﬁg. 12(b) ĥ(qn , t) changes its sign around(2)t ≈ 20 [ms]. Hence ĥ − ĥe (ρ0 , φ0 ) decays exponentiallyEur. Phys. J. E (2016) 39: 52−1−1
with the rate γa1 for t  γa2( γa3). However, ﬁg. 13(a)(2)
shows a slight deviation between ĥ − ĥe (ρ0 , φ0 ) and−γa1 t. This is due to the fact that the ratio γa1 /γa2 = 7.15
eis not large enough to regard ρ as a completely frozen variable. The second mode γa2 in ﬁg. 13(b) is associated with(1)
the relaxation of ρ to ρe , while φ is frozen and ĥ rapidly(2)(1)
changes to ĥe . Hence we have ρ − ρe (φ0 ) ∼ e−γa2 t for−1−1
 t  γa1.
γa3

5 Summary and discussion
In this paper, we have theoretically investigated the relaxation dynamics of a binary lipid bilayer membrane bytaking into account i) the coupling between the height
and the density variables, ii) the hydrodynamics of thesurrounding ﬂuid, iii) the frictional force between the upper and lower leaﬂets, and iv) the mutual diﬀusion in eachmonolayer. In sect. 2, we have constructed the free energy
in terms of the membrane shape h, the total lipid density
ρ± , and the lipid density diﬀerence φ± up to quadratic order. The membrane surface tension σ, which was neglectedin the previous theory for single-component lipid bilayer
membranes [5], and taken into account recently [48], naturally appears in the expansion of the general free energy
in eq. (1).
In sect. 3, the dynamic equations have been formulatedon the basis of momentum and molecular number conservations. In appendix A, we have proved the non-negative
deﬁniteness of the dissipation in our formulation. We havealso presented an alternative derivation of the dynamic
equations by using the Onsager’s variational principle inappendix B. The derived equations for binary lipid bilayer
membranes are the generalization of those in the Seifert
and Langer model [5]. We have further obtained the relaxation equations for ﬁve variables by integrating out the velocity ﬁeld of the surrounding ﬂuid (see also appendix C).
The equations are separated into two independent sets ofequations; one for (ĥ, ρ, φ) and the other for (φ̄, ρ̄). The
former equations change their signs under the interchangeof the upper and lower leaﬂets, while the latter equations
are invariant.
In sect. 4, we have discussed the stability of the onephase state and found that there are two possible instabilities; the anti-registered instability of (ĥ, ρ, φ) and the registered instability of (φ̄, ρ̄). We have investigated in detailthe relaxation rates of the various hydrodynamic modes.
In the case of small surface tension σ < σt (see eq. (76)),ﬁgs. 5 and 6 show that the mode crossing between ρ andĥ takes place around the intermediate wave number qmc .
Such a mode crossing was originally predicted for tensionless single-component lipid bilayers [5]. When σ > σt ,however, the height variable ĥ is the fastest mode in thewhole wave number range, and the mode crossing does notoccur [48].
Unlike single-component membranes for which eitherĥ or ρ is the slowest mode, mutual diﬀusion in two-

Page 15 of 21

component membranes is the slowest mode both for small
and moderate surface tensions. While φ varies slowly, thefaster variables ĥ and ρ rapidly approach their respective
quasi-equilibrium states determined by φ. In all the examined cases, the eﬀective decay rate γa3  γφ∗ for φ (seeeq. (69)) is smaller than the bare decay rate (Γa )33 because of the faster slaved variables ĥ and ρ.
As the unstable region is approached, the slowdownof the eﬀective rate γa3 becomes even more signiﬁcant,and the crossover from γa3 ∼ q 2 to ∼ q 4 behaviors may bemeasurable in experiments. As for the faster dynamics, therelaxation of ĥ is controlled by the hydrodynamics of thesurrounding ﬂuid, and the corresponding decay rate is approximately given by A11 /(4ηd2 q) (see eqs. (49) and (67)).
The relaxation of ρ is dominated by the inter-monolayer
friction, and its decay rate is given by A22 q 2 /(4b) (see thesentences after eq. (67)). We have also examined the relaxation of a domain that is rich in φ when the membrane
is close to the unstable region. In the very early stage, thebending of the membrane is induced by a non-zero density variation of φ even the membrane is initially ﬂat. In
the late stage of the relaxation process, all the variables
decay with the common decay rate γa3  γφ∗ as mentioned
above.
The dynamics of (ρ̄, φ̄) is simpler than that of (ĥ, ρ, φ).
The fastest variable ρ̄ instantly approaches to its quasiequilibrium state ρ̄e . Then φ̄ relaxes with the eﬀective decay rate γb2  γφ̄∗ (see eq. (80)) which becomes even slower
as the unstable region is approached.
While the kinetic parameters and some of the staticparameters have been determined in sect. 4, the dimensionless parameters Λi (i = 1, 3, 4, 5) in the free energy
could not be estimated from the previous experimental
data. However, the behaviors of the relaxation rates, which
are summarized in table 2 and are described in ﬁgs. 5–10,
are not sensitive to these parameters, unless the reduced
temperatures τa and τb are very close to zero (τa or τb
are deﬁned in terms of Λi ’s in eqs. (71) and (84), respectively). In fact, besides the parameters determined fromthe experimental data, these reduced temperatures are theonly relevant parameters. In the case of Γa (resp. Γb ), thisis because the time scales of the diﬀerent modes characterized by the diagonal elements (Γa )ii (resp. (Γb )ii )are well separated, except in the vicinity of the characteristic wave number q  qmc where the values of two
fastest modes of Γa become close in the low tension case(qmc is independent of the parameters that could not be
estimated). The two reduced temperatures measure thedistances in the phase space from their respective critical points3 , and one can experimentally control them byvarying the average lipid composition and the temperature. As discussed above, when τa (resp. τb ) is close tozero, the anti-registered (resp. registered) diﬀusive mode
becomes very slow, and the associated rate is given by γφ∗(resp. γφ̄∗ ).
Finally, we give some remarks. i) We have constructed
our free energy as a power series expansion up to quadraticorder with respect to the deformation and the densities

Page 16 of 21

Eur. Phys. J. E (2016) 39: 52

about the reference state. Here the physical meaning andmicroscopic interpretation of some phenomenological coupling parameters such as Λi are not so obvious. It would be
ideal to construct a free energy from a microscopic model,
and perform a series expansion of the free energy withrespect to the densities and curvature. With such a procedure, a connection between our phenomenological parameters and the microscopic quantities can be made. Recently,
an attempt has been made for a ﬂat bilayer membrane by
Williamson and Olmsted who derived a mean ﬁeld freeenergy from a semi-microscopic lattice bilayer model. In
their model, the diﬀerence in length between the two different lipid species was taken into account [49].
ii) In real biological cells, inclusions in membranes suchas proteins play essential roles. It was recently discussed
that the proteins which span the bilayer give rise to a further constraint in the dynamics and an additional sourceof dissipation leading to anomalous diﬀusion [50]. Furthermore, the surrounding ﬂuid can be viscoelastic rather thanpurely viscous, and inclusions can be active in a sense thatthey consume energy and drive membranes out of equilibrium. Neglecting the bilayer structure, some authors
have investigated the membrane shape ﬂuctuations whenit contains active/non-active inclusions and is surroundedby viscoelastic media [51–53]. Generalization of our theory
to such situations is also interesting.
iii) As we further approach the unstable region or thecritical point, the dynamical non-linear coupling (modemode coupling) between the density variables and the velocity ﬁelds in the bilayer becomes important like in theordinary 3D critical ﬂuids [15–18]. It would be interesting
to investigate the eﬀects of the bilayer structure and friction on top of the mode coupling between the velocity andthe density ﬁelds.

Next we examine the contribution from the change of
ρ± . Using eqs. (38) and (43), we obtainδF ∂ρ±= − d2 x f ± · ṽ ±(Ḟ )ρ± ≡ d x ±
δρ ∂t


#↔±
$
= d2 x −D̃v± ± ṽ ± · T · ez − b(ṽ + − ṽ − ) ,


2

(A.2)
where ez is the unit vector in the z-direction, and D̃v± isthe viscous dissipation in the monolayers
D̃v± =

±(∂j ṽi± )τij
=

ij

μ

˜ · ṽ ±∂i ṽj± + ∂j ṽi± − δij ∇2
2
˜ · ṽ ± .
+ζ ∇ij

Appendix A. Dissipation function
In this appendix, we discuss the dissipation which is related to the change rate of the free energy, W = −dF/dt.
The contribution due to the mutual diﬀusion is given by
the change of φ±(Ḟ )φ± ≡d2 x

δF ∂φ±=−
δφ± ∂td2 x|jφ± |2
Lφ,where use has been made of eqs. (39) and (40).

(A.1)

(A.3)

From the boundary conditions eqs. (45) and (46), thevelocity in the monolayers ṽ ± can be expressed in terms of
the surrounding ﬂuid velocity v ± as ṽ ± = v ± −(∂h/∂t)ez .
Using this relation with eqs. (31) and (32), we obtain


↔±





∂h ±
T ,
∂t zz
±
(A.4)

where ± d3 x denotes the 3D integration in the ranges ofz > 0 and z < 0, and Dv is the viscous dissipation in thesurrounding ﬂuid
d2 x ṽ ± · T· ez = ∓

Dv =

η
ij

We thank D. Andelman, T. Hoshino, T. Kato, C.-Y.D. Lu,
P.D. Olmsted, P. Sens, M. Turner, K. Yasuda for useful discussions. R.O. and S.K. acknowledge support from the Grantin-Aid for Scientiﬁc Research on Innovative Areas “Fluctuationand Structure” (Grant No. 25103010) from the Ministry of Education, Culture, Sports, Science, and Technology of Japan, theGrant-in-Aid for Scientiﬁc Research (C) (Grant No. 15K05250)
from the Japan Society for the Promotion of Science (JSPS),
and the JSPS Core-to-Core Program “International Research
Network for Non-equilibrium Dynamics of Soft Matter ”.

22d3 x Dv −d2 x

2(∂i vj + ∂j vi ) .

(A.5)

Furthermore, we deﬁne the dissipation due to the friction
between the two monolayers
D̃f = b|ṽ + − ṽ − |2 .

(A.6)

Combining eqs. (A.1), (A.2) and (A.4), we ﬁnally obtain


 
δF ∂h
−
(Ḟ )ρ + (Ḟ )φ
W = − d2 xδh ∂t =+,−




 2
D̃v + L−1|j|= d2 x D̃f +
φ
φ
+

 
=+,−

=+,−

d3 x Dv .

(A.7)



Here we see that the dissipation occurs through i) the friction between the monolayers, ii) the shear and bulk viscosities of the monolayers, iii) the mutual diﬀusion in themonolayers, and iv) the shear viscosity of the surrounding ﬂuid. The positivity of η, μ, ζ, Lφ and b ensures the
positivity of the dissipation; W ≥ 0.

Eur. Phys. J. E (2016) 39: 52

Page 17 of 21

Appendix B. Onsager’s variational principle
In appendix A, we have derived the dissipation eq. (A.7),
starting from the dynamic equations given by eqs. (31),(32), (38), (39), (40), (43), (44), (45) and (46). Conversely,
these dynamic equations can be obtained by the variational principle provided that we know the dissipations,
namely, the right-hand side of eq. (A.7). This is called
Onsager’s variational principle [33–35]. It is applicable tomany dynamical problems in soft matter such as colloidal
dispersions, membranes and polymer solutions if inertial
eﬀects can be neglected [35, 48, 50, 54, 55].
More precisely, we can derive eqs. (31), (40), (43)and (44) by minimizing the Rayleighian
dF
W
+
2
dt


  δF
δF
W
δF
2
+ d x
ρ̇ +  φ̇
=
ḣ +

2
δh
δρ
δφ
=+,−

R=

,
(B.1)

with respect to v, ṽ ± , jφ± , ḣ = ∂h/∂t, ρ̇± = ∂ρ± /∂t andφ̇± = ∂φ± /∂t.
The incompressible condition of the surrounding ﬂuid
(eq. (32)), the continuity equations (eqs. (38) and (39)),
and the non-slip boundary conditions (eqs. (45) and (46))
are taken into account as the constraints under which theRayleighian is minimized. Hence we minimize the shifted
Rayleighian 
d3 x p∇ · v
R∗ = R −


=+,−

d2 x

+



 
=+,−


+

d2 x



˜ · ṽ  ) + Cφ (φ̇ + ∇˜ · jφ )Cρ (ρ̇ + ∇



⎡

⎣Cvz
(ḣ − vz ) +

=+,−



⎤

Cvi
(ṽi − vi )⎦ ,i=x,ywith respect to v, ṽ ± , jφ± , ḣ, ρ̇± , φ̇± and the Lagrange±
multipliers p, Cρ± , Cφ± , Cvi(i = x, y, z).
We ﬁrst consider an inﬁnitesimal variation of v as v →
v + δv. Then the ﬁrst variation δR∗ |v with respect to v
is given by 
d3 x −∇p + ∇2 v + ∇(∇ · v) · δv
δR∗ |v = −
=+,−

−

d2 x









[Cvi
+ Tiz] δvi ,Similarly, minimization of R∗ with respect to ṽ ± , jφ± , ḣ,ρ̇± and φ̇± yields±±
Cvi− ∂j τij− ∂i Cρ± ± b(ṽi+ − ṽi− ) = 0,

(B.5)

jφ± = Lφ ∇Cφ± ,

(B.6)

+
−
+ Cvz
=−
Cvz

(B.3)

=+,− i=x,y,z

where δvi± is the velocity variation δvi evaluated at z →±0. Hence the minimization of R∗ with respect to v (inthe bulk region) and p yields the Stokes equation for the

δF
,
δh

(B.7)

Cρ± = −

δF
,
δρ±

(B.8)

Cφ± = −

δF
,
δφ±

(B.9)

respectively. Substituting eqs. (B.4) and (B.8) into eq.
(B.5), we obtain the force balance equation of the upper and lower monolayers in the tangential direction,eq. (43). The force balance equation in the normal direction, eq. (44), is obtained by substituting eq. (B.4) intoeq. (B.7). Finally, substitution of eq. (B.9) into eq. (B.6)
yields the diﬀusive ﬂux, eq. (40).

Appendix C. Elimination of the velocity ﬁeld
In this appendix, we discuss the dynamics of a single
Fourier mode. Without loss of generality, we can take the(x, y)-coordinate so that the direction of the wave vectorq̃ coincides with the x-direction, i.e., q̃ = (qx , 0) withqx > 0. Then we have q = |q̃| = qx . Substitution ofp(x, z) = p(q, z)eiqx and v(x, z) = v(q, z)eiqx into eq. (31)and (32) gives
iqp + η(q 2 − ∂z2 )vx = ∂z p + η(q 2 − ∂z2 )vz = 0,

(C.1)

− ∂z2 )vy = 0,

(C.2)

∂z vz + iqvx = 0.

(C.3)

(q
(B.2)



surrounding ﬂuid, eq. (31). Here the Lagrange multiplierp can be identiﬁed as the pressure ﬁeld.
Furthermore, minimizing R∗ with respect to v at z =±0, we obtain±±= ∓Tiz
.
(B.4)
Cvi

2

We then solve these equations to havep = 2ηR± e∓qz ,
±

(C.4)
∓qz

vz = (R z + S)e
vx =

,

i ±
[R ∓ q(R± z + S)]e∓qz ,qvy = Q± e∓qz ,

(C.5)
(C.6)
(C.7)

where Q± , R± and S are integral constants, and the upper
and the lower signs indicate the ﬂuids for z > 0 and z < 0,respectively. In deriving eqs. (C.5) and (C.6), we have used
the boundary condition vz+ = vz− (eq. (46)).
Next we substitute ρ± (x) = ρ± (q)eiqx , φ± (x) =±φ (q)eiqx , h(x) = h(q)eiqx and ṽ ± (x) = ṽ ± (q)eiqx into

Page 18 of 21

Eur. Phys. J. E (2016) 39: 52

eqs. (43) and use eqs. (45) and (C.4)–(C.7). After somealgebra we obtain
R+ − R− − 2qS = −2iqc1 Fq̃ [fx+ − fx− ],(C.8)

R+ + R− = −2iqc2 Fq̃ [fx+ + fx− ],

(C.9)

±

Q = 0,(C.10)

where Fq̃ [· · · ] denotes the Fourier transform at wave number q̃ = (q, 0), and c1 and c2 are deﬁned in eqs. (51)and (52), respectively. Similarly, from eqs. (44), (C.4)
and (C.5), we obtain S = −Fq̃ [δF/δh]/(4ηq). Furthermore we use eqs. (46) and (C.5) to have ∂h(q)/∂t =vz (q, 0) = S. Then the time evolution of h(q) is given
by


1
δF
∂h
=−
Fq̃
.
(C.11)∂t4ηqδhSubstituting ρ± (x) = ρ± (q)eiqx and ṽ ± (x) = ṽ ± (q)eiqx into eq. (38) and using eqs. (45) and (C.6), we have∂ρ± (q)/∂t = −iqṽx± = R± ∓ qS. Then ρ(q) and ρ̄(q) deﬁned in eq. (14) obey the following equations:∂ρ= −iqc1 Fq̃ [fx+ − fx− ],∂t∂ ρ̄= −iqc2 Fq̃ [fx+ + fx− ],∂t(C.12)(C.13)

where use has been made of eqs. (C.8) and (C.9). The time
evolution of φ± (q) in eq. (39) now reads


δF
∂φ±
2
= −Lφ q Fq̃
.
(C.14)∂tδφ±in eq. (18). We deﬁne the susceptibilities χa (q) and χb (q)as the reciprocals of the eigenvalues of A(q) and B(q),respectively. For large wave numbers, they behave as(κ + kd2 Ω0 )d2 + σd2 q −2 + O(q −4 ),1
4=q ×
χa (q)
cM ± q −2 + O(q −4 ),
(D.1)
and
1
= q 2 [cM ± + O(q −2 )],(D.2)
χb (q)
where
M ± = 1 + λ1 ±

*
k) ±
δF
2α + Λ2 β ± + Λ3 α∓ + Λ5 β ∓ ,
=
±
δρ
2

c  ˜2 ±˜ 2 φ± ,−2∇ ρ + λ2 ∇
2

(D.3)

where χ±a are given by


+
1
1
2 + 4A2
A
=
+
A
±
(A
−
A
)
22
33
22
33
23 , (D.5)
2χ±aand Aij are evaluated at q = 0. In eq. (D.4), the ﬁrst line
vanishes as q → 0. This zero eigenvalue at q = 0 corresponds to the homogeneous translation of the membranein the z-direction, which costs no energy.
The eigenvalues of B are given by1
1= ± + O(q 2 ),χb (q)
χb

(D.6)

where


+
1
1
2 + 4B 2
B
=
+
B
±
(B
−
B
)
11
22
11
22
12 , (D.7)
2χ±b(C.15)
where Bij are evaluated at q = 0. The thermodynamic−
stability conditions at q = 0 are given by χ−a , χb > 0+−−(note that χ+
a , χb > 0 if χa , χb > 0). These inequalities
are ensured if and only if eqs. (54)–(57) are satisﬁed.

and
*
k)
δF
2Λ1 β ± + Λ2 α± + Λ4 β ∓ + Λ5 α∓
=
±
δφ
2

c
˜ 2 φ± + λ 2 ∇˜ 2 ρ± ,−2λ1 ∇2(1 − λ1 )2 + λ22 .

The thermodynamic stability in large wave numbers is
ensured by M ± > 0, which is equivalent to eq. (53).
For small wave numbers, the eigenvalues of A are givenby
σd2 q 2 + κd2 q 4 + O(q 6 ),
1
=
(D.4)
2
χa (q)(1/χ±a ) + O(q ),Functional derivatives of eqs. (11) and (12) are givenby+(C.16)

respectively. Using eqs. (42), (44), (C.15) and (C.16), we
can calculate the right-hand sides of eqs. (C.11)–(C.14) toobtain eqs. (47)–(50).

Appendix D. Thermodynamic stability
Appendix D.1. Stability at q → 0 and q → ∞The static ﬂuctuations and stability of the system is characterized by the eigenvalues of the matrices A(q) and B(q)Appendix D.2. Stability at intermediate wave numbers
Assuming that the modes at q = 0 and ∞ arethermodynamically stable, we have A22 (q), B11 (q) >0 from eqs. (53) and (54). Then, ĥ, ρ and ρ̄
can be integrated out from the Boltzmann weightexp(−F [ĥ, ρ, φ, ρ̄, φ̄]/kB T ) to obtain an eﬀective free energy Feﬀ [φ, φ̄] = −kB T ln[ Dĥ Dρ Dρ̄ e−F/kB T ], where
DĥDρ Dρ̄ denotes the functional integral and kB is theBoltzmann constant. Since our free energy F is quadratic,
Feﬀ can be obtained by minimizing F with respect to ĥ,ρ and ρ̄. By equating the right-hand sides of eqs. (44)Eur. Phys. J. E (2016) 39: 52

Page 19 of 21

10 4

Table 3. Values of dimensionless quantities in eqs. (D.13),(D.16) and (D.18) evaluated for the parameter values in sect. 4
and Δa = Δλ = ν = 1.
s3

τa∗

10

−2

2.86 × 10

5

7 × 10106

2.86 × 10−4

10−8

2.86 × 10−2

7 × 1091010

2.86 × 10−62

σ (erg/cm )
−4

s1

s2

(a)

10 3
10 2
10 -8
10
10 -6 10 -4

and (C.15) to zero, we can eliminate ĥ, ρ and ρ̄ fromeq. (18) to obtain


|φ̄|2
d2 q 1 |φ|2
+
.
(D.8)
Feﬀ =
2
(2π) 2 χφ (q) χφ̄ (q)

10 -2

(b)

1

Here the susceptibilities χφ (q) and χφ̄ (q) for φ and φ̄ aregiven by
det A
1
=
,
(D.9)
χφ (q)A11 A22 − A212

10 2

1

0.9
0.8

and
det B
1
=
,
χφ̄ (q)B11

respectively.
We expand 1/χφ and 1/χφ̄ in powers of q to have)*1
= τa + s1 Δ2a + Δλ − s2 (ντa )2 (qd)2 + O(q 4 ),kχφ (q)(D.11)
and)*1= τb + s1 Δ2b + Δλ (qd)2 + O(q 4 ),kχφ̄ (q)2ck,s2 = ,2kd
σΛ2 − Λ5λ2Δa =,−2 − Λ3
2

Cφ  s2 (s3 + s2 ν 2 τa )(ντa )2 ,s3 =(D.15)

1/2

10 2(D.17)

where we have deﬁned

Notice that the stability conditions at q = 0 in eqs. (54)(57) are equivalent to the conditions τa , τb , A22 (0),B11 (0) > 0.
In eq. (D.12), the term quadratic in q is always positive, which indicates that φ̄ and ρ̄ do not exhibit any
instability at intermediate wave numbers when they arestable at q = 0. In eq. (D.11), however, the quadraticterm is negative if
s1 (Δ2a + Δλ )τa >s2 ν 2

1

In this case, we need to add a quartic term Cφ (qd)4 toeq. (D.11) in order to examine the stability at intermediate wave numbers. The coeﬃcient of the quartic term isobtained as(D.14)

Λ2 + Λ5
λ2.
−
2 + Λ32



10 -2

Fig. 14. Plots of (a) qmin deﬁned by eq. (D.20) as a functionof τa − τa∗ , and (b) χφ (0)/χφ (qmin ) as a function of τa . The
parameters are the same as in table 3. The reciprocal of χφ
exhibits a minima at q = qmin (> 0) for τa > τa∗ . We see in (b)that the minimum value 1/χφ (qmin ) does not become zero as
far as 1/χφ (0) = kτa > 0. The two curves in (b) for diﬀerent σvalues coincide almost completely for τa > τa∗ = 2.86 × 10−4 .

(D.13)

s1 =

10 -4

(D.12)

where the reduced temperatures τa and τb were deﬁned
in eqs. (71) and (84), respectively. In the above, we havedeﬁned dimensionless combinations

Δb =0.7 -6
10(D.10)≡ τa∗ .

(D.16)

κ.
σd2(D.18)In deriving eq. (D.17), we have assumed s1  s2 , s3 . In
table 3, we list the values of s1 , s2 , s3 and τa∗ for theparameter values chosen in sect. 4 and Δa = Δλ = ν = 1.
For τa > τa∗ , the reciprocal of the susceptibility 1/χφhas a minima at an intermediate wave number q = qmin ,s2 ν 4 [τa2 − (τa∗ )2 ]2
1= τa − 2,kχφ (qmin )4Cφ(D.19)

1/2
1 s2 ν 2 {τa2 − (τa∗ )2 }qmin =.
d2Cφ(D.20)

wherePage 20 of 21

Eur. Phys. J. E (2016) 39: 52

From table 3, we can assume τa∗ ∼s1 /s2  1 ands2 /s3 = kd2 /κ to be of order of unity so that1kχφ (qmin )⎧⎪τ(τa − τa∗  s3 /s2 ),
⎪
⎨ a


1
τa
⎪
⎪
3+
(τa  τa∗ ),⎩4
1 + (s2 /s3 )ν 2 τa(D.21)
andqmin 

⎧⎪⎪⎨

τa − τa∗(s3 + s2 ν 2 τa∗ )τa∗

1/2

1
×
d ⎪
⎪
⎩[2(s + s ν 2 τ )]−1/232a(τa − τa∗  τa∗ ),(τa  τa∗ ).

(D.22)
Since we see 1/χφ (qmin ) > 0 in eq. (D.21), the instability
at q = qmin does not takes place in both regimes.
In ﬁg. 14(a), we plot qmin as a function of τa − τa∗ .
For both σ = 10−8 and 10−4 erg/cm2 , we have√ qmin ∼(τa − τa∗ )1/2 for τa − τa∗  τa∗ , qmin  1/(d 2s3 ) for−1/2for τa τa∗  τa  s3 /(s2 ν 2 ), and qmin ∼ τa2s3 /(s2 ν ). These behaviors are in good agreement witheq. (D.22). In ﬁg. 14(b), we plot 1/χa (qmin ) multiplied byχa (0) = 1/(kτa ) as a function of τa , where the parameters are the same as in (a). We see that the curves fordiﬀerent σ values almost coincide and in agreement witheq. (D.21). The ratio χφ (0)/χφ (qmin ) monotonically decreases from unity to the lower bound 3/4 as τa increases.
Therefore we conclude that instability does not occur atintermediate wave numbers when the modes are stable at
q = 0. Therefore the overall thermodynamic stability isensured by the conditions in eqs. (54)–(57).
Leibler and Andelman discussed the instability in
two-component membranes at intermediate wave numbers [56]. Although they did not explicitly take into account the bilayer structure, their model is similar to ours
because the composition-bending coupling is explicitly
taken into account. However, they treated the coeﬃcientsof the power series of the susceptibility as independent
parameters. In our study, the coeﬃcient in eq. (D.11)
and eq. (D.17) vary simultaneously when Λi values arechanged.

References
1. B. Alberts, A. Johnson, P. Walter, J. Lewis, M. Raﬀ,
Molecular Biology of the Cell (Garland Science, New York,
2008).
2. R. Lipowsky, E. Sackmann (Editors), Structure and Dynamics of Membranes – From Cells to Vesicles (Elsevier,
Amsterdam, 1995).
3. L. Kramer, J. Chem. Phys. 55, 2097 (1971).
4. F. Brochard, J.F. Lennon, J. Phys. (Paris) 36, 1035 (1975).
5. U. Seifert, S.A. Langer, Europhys. Lett. 23, 71 (1993).

6. A. Yeung, E. Evans, J. Phys. II 5, 1501 (1995).
7. W. Pfeiﬀer, S. König, J.F. Legrand, T. Bayerl, D. Richter,
E. Sackmann, Europhys. Lett. 23, 457 (1993).
8. T. Pott, P. Méléard, Europhys. Lett. 59, 87 (2002).
9. S.A. Shkulipa, W.K. den Otter, W.J. Briels, Phys. Rev.
Lett. 96, 178302 (2006).
10. J.-B. Fournier, N. Khalifat, N. Puﬀ, M.I. Angelova, Phys.
Rev. Lett. 102, 018102 (2009).
11. A.-F. Bitbol, N. Puﬀ, Y. Sakuma, M. Imai, J.-B. Fournier,
M.I. Angelova, Soft Matter 8, 6073 (2012).
12. S.L. Veatch, S.L. Keller, Biochim. Biophys. Acta 1746,
172 (2005).
13. A.R. Honerkamp-Smith, S.L. Veatch, S.L. Keller, Biochim.
Biophys. Acta 1788, 53 (2009).
14. S. Komura, D. Andelman, Adv. Colloid Interface Sci. 208,
34 (2014).
15. A.R. Honerkamp-Smith, B.B. Machta, S.L. Keller, Phys.
Rev. Lett. 108, 265702 (2012).
16. K. Seki, S. Komura, M. Imai, J. Phys.: Condens. Matter
19, 072101 (2007).
17. K. Inaura, Y. Fujitani, J. Phys. Soc. Jpn. 77, 114603(2008).
18. S. Ramachandran, S. Komura, K. Seki, M. Imai, Soft Matter 7, 1524 (2011).
19. S. May, Soft Matter 5, 3148 (2009).
20. M.D. Collins, Biophys. J. 94, L32 (2008).
21. M.D. Collins, S.L. Keller, Proc. Natl. Acad. Sci. U.S.A.
105, 124 (2008).
22. J. Zhang, B. Jing, N. Tokutake, S.L. Regen, J. Am. Chem.
Soc. 126, 10856 (2004).
23. J. Zhang, B. Jing, V. Janout, S.L. Regen, Langmuir 23,
8709 (2007).
24. W.-C. Lin, C.D. Blanchette, T.V. Ratto, M.L. Longo, Biophys. J. 90, 228 (2006).
25. M.J. Stevens, J. Am. Chem. Soc. 127, 15330 (2005).
26. J.D. Perlmutter, J.N. Sachs, J. Am. Chem. Soc. 133, 6563
(2011).
27. A.J. Wagner, S. Loew, S. May, Biophys. J. 93, 4268 (2007).
28. G.G. Putzel, M. Schick, Biophys. J. 94, 869 (2008).
29. Y. Hirose, S. Komura, D. Andelman, Chem. Phys. Chem.
10, 2839 (2009).
30. Y. Hirose, S. Komura, D. Andelman, Phys. Rev. E 86,
021916 (2012).
31. S.R. de Groot, P. Mazur, Non-Equilibrium Thermodynamics (Dover, New York, 1984).
32. L.D. Landau, E.M. Lifshitz, Fluid Mechanics, 2nd edition(Pergamon, Oxford, 1987).
33. L. Onsager, Phys. Rev. 37, 405 (1931).
34. L. Onsager, Phys. Rev. 38, 2265 (1931).
35. M. Doi, J. Phys.: Condens. Matter 23, 284118 (2011).
36. A.-F. Bitbol, L. Peliti, J.-B. Fournier, Eur. Phys. J. E 34,
53 (2011).
37. J. Song, R.E. Waugh, Biophys. J. 64, 1967 (1993).
38. G. Niggemann, M. Kummrow, W. Helfrich, J. Phys. II 5,
413 (1995).
39. W. Rawicz, K.C. Olbrich, T. McIntosh, D. Needham, E.
Evans, Biophys. J. 79, 328 (2000).
40. A. Tian, C. Johnson, W. Wang, T. Baumgart, Phys. Rev.
Lett. 98, 208102 (2007).
41. E.A. Evans, R. Waugh, L. Melnik, Biophys. J. 16, 585(1976).
42. E. Evans, W. Rawicz, Phys. Rev. Lett. 64, 2094 (1990).

Eur. Phys. J. E (2016) 39: 52
43. P.F.F. Almeida, W.L.C. Vaz, Structure and Dynamics ofMembranes – From Cells to Vesicles, edited by R. Lipowsky, E. Sackmann (Elsevier, Amsterdam, 1995) Chapt. 6.
44. R. Merkel, E. Sackmann, E. Evans, J. Phys. France 50,
1535 (1989).
45. A. Horner, S.A. Akimov, P. Pohl, Phys. Rev. Lett. 110,
268101 (2013).
46. A. Onuki, Phase Transition Dynamics (Cambridge University Press, Cambridge, 2002).
47. M. Mell, L. Moleiro, Y. Hertle, I. Lôpez-Montero, F.J. Cao,
P. Fouquet, T. Hellweg, F. Monroy, Chem. Phys. Lipids
185, 61 (2015).

Page 21 of 21
48. J.-B. Fournier, J. Non-Linear Mech. 75, 67 (2015).
49. J.J. Williamson, P.D. Olmsted, Biophys. J. 108, 1963
(2015).
50. A. Callan-Jones, M. Durand, J.-B. Fournier, Soft Matter
12, 1791 (2016).
51. R. Granek, S. Pierrat, Phys. Rev. Lett. 83, 872 (1999).
52. D. Lacoste, A.W.C. Lau, Europhys. Lett. 70, 418 (2005).
53. S. Komura, K. Yasuda, R. Okamoto, J. Phys.: Condens.
Matter 27, 432001 (2015).
54. M. Arroyo, A. DeSimone, Phys. Rev. E 79, 031915 (2009).
55. M. Rahimi, M. Arroyo, Phys. Rev. E 86, 011932 (2012).
56. S. Leibler, D. Andelman, J. Phys. (Paris) 48, 2013 (1987).