1. Introduction
Flow around biological membranes holds a great deal of importance in a multitude of systems (Seifert Reference Seifert1999; Herzenberg et al. Reference Herzenberg, Tung, Moore, Herzenberg and Parks2006; Shi et al. Reference Shi, Graber, Baumgart, Stone and Cohen2018; Salmond et al. Reference Salmond, Khanna, Owen and Williams2021). These biological membranes are modelled using surrogate structures known as vesicles. Vesicles are complex droplets having a lipid bilayer on the surface instead of simple fluid boundaries. The presence of this lipid bilayer imparts an elastic resistance to bending (Helfrich Reference Helfrich1973). These lipid bilayer membranes act like two-dimensional (2-D) fluid sheets that resist any stretching or compression (Campelo et al. Reference Campelo, Arnarez, Marrink and Kozlov2014). Apart from being a marvellous surrogate system for understanding the biophysics of cellular membranes, lipid vesicles are often used as carriers for drug delivery processes (Needham Reference Needham1999; Guo & Szoka Reference Guo and Szoka2003).
Lipid bilayers are made up for long chain compounds known as phospholipids that contain a hydrophilic head and a hydrophobic tail (Alberts Reference Alberts2017). When the bilayer consists of a single type of phospholipid, we call them ‘single component’ vesicles. These vesicles have been the centre of a plethora of studies over the past five decades. The physical properties and thermodynamics of such vesicles can be characterized by a lipid bilayer exhibiting a uniform bending resistance (Lipowsky & Sackmann Reference Lipowsky and Sackmann1995). Previous studies have explored the dynamics of such single component vesicles in great detail by studying multiple modes of vesicle motion in different systems: (i) shear flow – tank-treading (Keller & Skalak Reference Keller and Skalak1982; Seifert Reference Seifert1999; Abkarian, Lartigue & Viallat Reference Abkarian, Lartigue and Viallat2002), swinging (Noguchi & Gompper Reference Noguchi and Gompper2007), tumbling (Biben & Misbah Reference Biben and Misbah2003; Rioual, Biben & Misbah Reference Rioual, Biben and Misbah2004; Noguchi & Gompper Reference Noguchi and Gompper2005; Kantsler & Steinberg Reference Kantsler and Steinberg2006) and vacillated breathing (Misbah Reference Misbah2006); (ii) extensional flow (Zhao & Shaqfeh Reference Zhao and Shaqfeh2013; Narsimhan et al. Reference Narsimhan, Spann and Shaqfeh2014, Reference Narsimhan, Spann and Shaqfeh2015; Dahl et al. Reference Dahl, Narsimhan, Gouveia, Kumar, Shaqfeh and Muller2016; Kumar, Richter & Schroeder Reference Kumar, Richter and Schroeder2020); (iii) general linear flows (Vlahovska & Gracia Reference Vlahovska and Gracia2007; Lin & Narsimhan Reference Lin and Narsimhan2019); (iv) oscillatory flows (Lin et al. Reference Lin, Kumar, Richter, Wang, Schroeder and Narsimhan2021).
The deep understanding of single component vesicles has provided a base for further explorations into systems that are closer to reality. Often, these lipid bilayer membranes contain multiple phospholipids along with cholesterol molecules interspersed between them (John et al. Reference John, Schreiber, Kubelt, Herrmann and Müller2002). The existence of multiple molecules in these bilayers makes for an incredible amalgamation of phase equilibrium thermodynamics and fluid mechanics (Safran Reference Safran2018). Some of these phospholipids (e.g. dipalmitoylphosphatidylcholine (DPPC)) have a larger affinity towards cholesterol molecules than others (e.g. 1,2-Dioleoyl-sn-glycero-3-phosphocholine (DOPC)) (Veatch & Keller Reference Veatch and Keller2005a; Davis, Clair & Juhasz Reference Davis, Clair and Juhasz2009; Uppamoochikkal, Tristram-Nagle & Nagle Reference Uppamoochikkal, Tristram-Nagle and Nagle2010). This preferential separation into phases leads to the formation of ordered and disordered liquid phases on the membrane surface (Veatch & Keller Reference Veatch and Keller2005b), a term known as ‘lipid rafts’ (Simons & Ikonen Reference Simons and Ikonen1997). These rafts have a relevance in signal transduction (Simons & Toomre Reference Simons and Toomre2000) and protein transfer across membranes, thereby having implications for health and diseases (Michel & Bakovic Reference Michel and Bakovic2007). These factors underscore the importance of studying such systems.
From a mechanical viewpoint, the existence of multiple phases imparts inhomogeneous properties like the bending rigidity to the vesicle due to the differences in bending stiffnesses of each constituent phospholipid (Baumgart, Hess & Webb Reference Baumgart, Hess and Webb2003). Moreover, the resultant lateral phases fall prey to a tussle between convective motion due to the background fluid and surface diffusive motion due to the inherent molecular properties of the phases (Yanagisawa et al. Reference Yanagisawa, Imai, Masui, Komura and Ohta2007; Taniguchi, Yanagisawa & Imai Reference Taniguchi, Yanagisawa and Imai2011; Arnold, Gubbala & Takatori Reference Arnold, Gubbala and Takatori2023).
While creating medical diagnostic devices, often, the measurement of mechanical properties of such vesicles is important (Kollmannsberger & Fabry Reference Kollmannsberger and Fabry2011; Lei et al. Reference Lei2021). These measurements help in improving the control and precision of devices. Previously, multiple numerical simulations have been performed to understand the motion of multicomponent vesicles under shear flow (Sohn et al. Reference Sohn, Tseng, Li, Voigt and Lowengrub2010; Liu et al. Reference Liu, Marple, Allard, Li, Veerapaneni and Lowengrub2017; Gera & Salac Reference Gera and Salac2018b). These studies have highlighted the influence of line tension and bending rigidity, along with the membrane tension of the vesicle, on the modes of motion that the vesicle undergoes – tank-treading, tumbling, swing/phase-treading, vertical banding, among others. More recently, under the limit of dominant advective forces, researchers came up with an analytical treatment of 2-D multicomponent vesicles and their swinging-to-tumbling transition (Gera, Salac & Spagnolie Reference Gera, Salac and Spagnolie2022). The transition primarily depended on a ratio of the bending stiffness of the two phases and the capillary number of the vesicle. While informative and thoughtful, the study lacked an analytical treatment of the case when diffusive time scales are comparable to that of the convective time scales. Secondly, the study treated a multicomponent vesicle as a 2-D inextensible membrane, thus leaving room for out-of-plane deformations. To the authors’ knowledge, such an analytical treatment has not been provided yet.
We aim to bridge this gap of knowledge through this study that focuses on the semianalytical prediction for a three-dimensional (3-D), nearly spherical, multicomponent vesicle placed under a background shear flow. We leverage the vector spherical harmonics-based techniques previously used for single component vesicles and drops (Vlahovska, Loewenberg & Blawzdziewicz Reference Vlahovska, Loewenberg and Blawzdziewicz2005; Vlahovska & Gracia Reference Vlahovska and Gracia2007) to solve the underlying nonlinear dynamical equations up to leading order, while ensuring conservation of the composition and vesicle surface area. This helps us arrive at reduced-order equations governing the shape and composition of the vesicle and the phospholipid–cholesterol phases, respectively. We delineate multiple motions exhibited by the vesicle and discuss their dependence on material properties. The aim of this study is to equip an experimental researcher with a theory that could a priori predict the vesicle dynamics based on the material specifications and control variables.
We discuss problem set-up, solution strategy and numerical methods in §§ 2 and 3. We discuss results for multicomponent vesicles in the high Péclet number limit in § 4, where membrane coarsening effects are unimportant. We then discuss results in the intermediate Péclet number limit in § 5, where coarsening effects come into play. We discuss implications of asymmetric phase distributions in § 6. We then conclude the study with a discussion in § 7.
2. Problem formulation
We consider a lipid membrane vesicle containing a ternary mixture of phospholipids and cholesterol undergoing phase separation on the membrane surface (see figure 1). This vesicle is placed in an unbounded domain containing a background shear flow
$\boldsymbol{u}^{\infty } = \dot {\gamma } y \boldsymbol{\hat {x}}$
. This membrane encloses a Newtonian fluid of viscosity
$(\lambda -1) \eta$
with a surrounding fluid of viscosity
$\eta$
, where
$\lambda -1$
is the viscosity ratio between the inner and outer fluids. The bilayer has a bending stiffness of
$\kappa _{c}$
that is dependent on the phospholipid distribution characterized by an order parameter
$q$
. This dependence will be explained in § 2.1. Furthermore, the lipid molecules impose a membrane tension to preserve the membrane surface area
$A$
and there exists a line tension between the separated phases (indicated by black and white coloured lipid heads). The vesicle has a characteristic size
$R$
which is the radius of a sphere of the same volume.
Schematic of system. The inset shows a zoomed-in version of the lipid bilayer and its properties (see § 2.1).
2.1. Free energy of membrane
The lipid membrane energy landscape has been a problem of interest for the past few decades. Historically, a lipid bilayer membrane is treated as an elastic surface. Multiple models have been used to describe such materials – the Keller Skalak model (Keller & Skalak Reference Keller and Skalak1982), area difference model (Seifert Reference Seifert1997) and Helfrich model (Helfrich Reference Helfrich1973; Campelo et al. Reference Campelo, Arnarez, Marrink and Kozlov2014) to name a few. In our study, we use a Canham–Helfrich model.
The free energy of a multicomponent membrane is governed by three factors: bending, mixing and surface tension energies. The bending energy is given by
\begin{equation} W_{\textit{bend}} = \int {\frac {1}{2}\kappa _{c} (2H-c_{0})^{2}{\textrm d}S}, \end{equation}
where
$\kappa _{c} = k_{0} + k_{1}q$
is a bending modulus that depends on the lipid phase behaviour characterized by an order parameter
$q$
. The order parameter
$q$
is a variable whose values can range from
$-1$
corresponding to pure
$L_{d}$
(cholesterol-lacking) and
$+1$
corresponding to pure
$L_{o}$
(cholesterol-rich). Thermodynamically speaking, this is the local coordinate along a tie line in the two-phase coexistence region that is present in the ternary phase diagram for a system like DOPC, DPPC and cholesterol. A tie line represents the line that connects the compositions of the ordered and disordered phases at equilibrium, as commonly seen in liquid–liquid extraction systems (Geankoplis Reference Geankoplis2003). The values of
$k_0$
and
$k_1$
are
$k_0 = {1}/{2} ( \kappa _{lo} + \kappa _{ld} )$
and
$k_1 = {1}/{2} ( \kappa _{lo} - \kappa _{ld} )$
, where
$\kappa _{l0}$
and
$\kappa _{ld}$
are the bending moduli of the
$L_0$
and
$L_d$
phases, respectively. Lastly, the parameters
$H$
and
$c_{0}$
are the mean and spontaneous curvatures of the bilayer leaflet. We treat
$c_{0}=0$
for a symmetric leaflet.
The mixing energy is given by the Landau–Ginzberg equation (Gompper & Klein Reference Gompper and Klein1992):
\begin{equation} W_{\textit{mix}} = \int {\left (\frac {\tilde {a}}{2}q^{2}+\frac {\tilde {b}}{4}q^{4} + \frac {\gamma ^{2}}{2}|\boldsymbol{\nabla} _{\!s}q|^{2}\right ){\textrm d}S}. \end{equation}
The first two terms represent a double-well potential that gives rise to phase separation (for
$\tilde {a} \lt 0$
), and the last term is an energy penalty for creating phases that is related to line tension. This free energy has been used in many studies of lipid systems undergoing phase separation. These parameters are related to the experimentally measured phase split concentrations (
$q^{\pm }$
), line tension
$(\xi ^{\textit{line}})$
and the interface width
$(\varepsilon ^{w\textit{idth}})$
as follows:
\begin{equation} q^{\pm } = \pm \sqrt {|\tilde {a}|/\tilde {b}}, \qquad \xi ^{\textit{line}} = \frac {2\sqrt {2}}{3\tilde {b}}|\tilde {a}|^{3/2}\gamma, \qquad \varepsilon ^{w\textit{idth}} = \sqrt {\frac {2\gamma ^{2}}{|\tilde {a}|}}. \end{equation}
The Landau–Ginzberg equation has been used to qualitatively model bilayer membranes (Gera & Salac Reference Gera and Salac2018b) – see the appendix of Camley & Brown (Reference Camley and Brown2014) for estimated values of
$\tilde {a},\tilde {b}$
and
$\gamma$
for a specific system (
$R\sim O(nm)$
). In their study, they mention that an approximate value for
$q^+ \approx 0.3$
would be appropriate for phase-split concentrations.
Lastly, the surface tension energy is given by
\begin{equation} W_{\sigma } = \int {\sigma {\textrm d}S} = \text{constant}. \end{equation}
The above contribution is constant because the lipid area per molecule is conserved and thus the membrane is incompressible. The surface tension
$\sigma$
thus acts as a Lagrange multiplier to ensure area conservation.
2.2. Fluid motion and membrane interface dynamics
We solve the Stokes equations inside and outside the vesicle that govern the velocity, stress and pressure distributions:
\begin{align}(\lambda -1)\eta {\nabla} ^{2}\boldsymbol{u}^{\textit{in}} &= \boldsymbol{\nabla }p^{\textit{in}}; \qquad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^{\textit{in}} = 0;\end{align}
\begin{align}\eta {\nabla} ^{2}\boldsymbol{u}^{\textit{out}} &= \boldsymbol{\nabla }p^{\textit{out}}; \qquad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^{\textit{out}} = 0.\end{align}
These equations are subject to the following boundary conditions: (i) far-field,
$\boldsymbol{u}^{\textit{out}} \rightarrow \boldsymbol{u}^{\infty }$
as
$|\boldsymbol{x}| \rightarrow \infty$
; (ii) continuity of velocity,
$\boldsymbol{u}^{\textit{in}} = \boldsymbol{u}^{\textit{out}}$
for
$\boldsymbol{x} \in S$
; (iii) membrane incompressibility,
$\boldsymbol{\nabla}_{\kern-1pt s} \boldsymbol{\cdot }\boldsymbol{u}^{\textit{in}} = 0$
for
$\boldsymbol{x} \in S$
; (iv) traction balance,
$\boldsymbol{f}^{\textit{hyd}} = \boldsymbol{f}^{\textit{mem}}$
for
$\boldsymbol{x} \in S$
.
In the traction balance,
$\boldsymbol{f}^{\textit{hyd}} = \boldsymbol{n} \boldsymbol{\cdot }(\boldsymbol{T}^{\textit{out}} - \boldsymbol{T}^{\textit{in}} )$
is the hydrodynamic traction from viscous stresses, where
$\boldsymbol{n}$
is the outward-pointing normal vector and
$\boldsymbol{T} = \mu ( \boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{u}^T ) - p \boldsymbol{I}$
is the stress tensor for a fluid with viscosity
$\mu$
(
$\mu = \eta$
outside the vesicle,
$\mu = (\lambda -1) \eta$
inside the vesicle). The membrane traction
$\boldsymbol{f}^{\textit{mem}}$
is equal to the first variation of the membrane free energy with respect to position:
$\boldsymbol{f}^{\textit{mem}} = \delta W/\delta \boldsymbol{x}$
. This can be broken into bending, mixing and surface tension contributions
$\boldsymbol{f}^{\textit{mem}} = \boldsymbol{f}^{\textit{bend}}+\boldsymbol{f}^{\textit{mix}}+ \boldsymbol{f}^{\sigma }$
, with expressions for each shown as follows:
\begin{align}\boldsymbol{f}^{\textit{bend}} &= \frac {\delta W_{\textit{bend}}}{\delta \boldsymbol{x}} = -\boldsymbol{n}{\nabla} ^{2}_{s}(\kappa _{c}(2H)) + \kappa _{c}(2H)\left [\boldsymbol{\nabla} _{\kern-2pt s}(2H) + (2K-4H^{2}) \right ]{\boldsymbol{n}}\nonumber \\ &\quad +\, 2H\left (\frac {1}{2}\kappa _{c}(2H)^{2}\right )\boldsymbol{n} - \boldsymbol{\nabla} _{\kern-2pt s}\left (\frac {1}{2}\kappa _{c}(2H)^{2}\right)\!,\end{align}
\begin{align}\boldsymbol{f}^{\textit{mix}} &= \frac {\delta W_{\textit{mix}}}{\delta \boldsymbol{x}}= \gamma ^{2} \boldsymbol{\nabla} _{\kern-2pt s}\boldsymbol{\cdot }(\boldsymbol{\nabla} _{\kern-2pt s}q {\otimes } \boldsymbol{\nabla} _{\kern-2pt s}q) + 2H\left (\frac {1}{2}\gamma ^{2}|\boldsymbol{\nabla} _{\kern-2pt s}q|^{2} + g(q)\right )\boldsymbol{n} \nonumber \\ & \quad-\, \boldsymbol{\nabla} _{\kern-2pt s}\left (\frac {1}{2}\gamma ^{2}|\boldsymbol{\nabla} _{\kern-2pt s}q|^{2} + g(q)\right)\!,\end{align}
\begin{align}\boldsymbol{f}^{\sigma } &= \frac {\delta W_{\sigma }}{\delta \boldsymbol{x}} = 2\sigma H\boldsymbol{n} - \boldsymbol{\nabla} _{\kern-2pt s}\sigma.\end{align}
The reader is directed to the following publications for details on how these equations are derived: Napoli & Vergori Reference Napoli and Vergori2010 and Gera Reference Gera2017. In the above expressions,
$\otimes$
represents a dyadic product,
$\boldsymbol{\nabla}_{\kern-1pt s} = ( \boldsymbol{I} - \boldsymbol{nn} ) \boldsymbol{\cdot }\boldsymbol{\nabla}$
is the surface gradient and
$g = ({\tilde {a}}/{2}) q^2 + ({\tilde {b}}/{4}) q^4$
is the double well potential in the mixing free energy (2.2). The quantities
$K = \text{det}(\boldsymbol{L}) = C_{1}C_{2}$
and
$H = ({1}/{2}) \text{tr}(\boldsymbol{L}) = ({C_{1}+C_{2}})/{2}$
are the Gaussian and mean curvatures of the interface, respectively, where
$\boldsymbol{L} = \boldsymbol{\nabla} _{\!s}\boldsymbol{n}$
is the surface curvature tensor. The surface tension
$\sigma$
is a Lagrange multiplier to enforce membrane incompressibility (
$\boldsymbol{\nabla}_{\kern-1pt s} \boldsymbol{\cdot }\boldsymbol{u} = 0$
on the interface), and hence must be solved in addition to the velocity and pressure fields.
In addition to the flow field, one must also solve for order parameter
$q$
on the membrane surface. The flow and coarsening behaviour of the order parameter satisfy a Cahn–Hilliard equation (Gera Reference Gera2017),
\begin{align}\frac {\partial q}{\partial t} + \boldsymbol{\nabla} _{\kern-2pt s}\boldsymbol{\cdot }(\boldsymbol{u} q) & = \frac {\nu }{\zeta _0} {\nabla} ^{2}_{s} \zeta,\end{align}
\begin{align}\zeta & = \frac {\delta W}{\delta q} = \left ( \tilde {a}q + \tilde {b}q^{3} - \gamma ^{2}{\nabla} ^{2}_{s}q + \frac {k_{1}}{2}(2H)^{2} \right )\!,\end{align}
where in the above expression,
$\nu$
is the mobility of the phospholipids (units of metres squared per second) and
$\zeta = \delta W/\delta q$
is the surface chemical potential (units of energy per area). The above expression states that lipids are either convected along the surface or move via gradients in chemical potential. The value reference chemical potential
$\zeta _0$
is given in Gera Reference Gera2017.
Lastly, we impose the kinematic boundary condition to avoid any slip on the membrane surface. If we parameterize the vesicle radius as
$r = r_s(\boldsymbol{x}, t)$
, the kinematic boundary condition is below, where
${\textrm D}/{{\textrm D}t} = ({\partial }/{\partial t}) + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla}$
is the material derivative:
\begin{equation} \frac {\textrm D}{{\textrm D}t} \left (r - r_s \right )= 0 \qquad \boldsymbol{x} \in S. \end{equation}
2.3. Time scales and dimensionless quantities
The vesicle dynamics occur over four time scales. The first one is a bending time scale
$t_{\textit{bend}} = R^{3}\eta /k_0$
that denotes the time taken by a vesicle to restore its equilibrium configuration under the action of bending forces. The second time scale is that of the flow
$t_{\dot {\gamma }} = \dot {\gamma }^{-1}$
.The third represents the time scale for coarsening
$t_{q} = R^{2}/\nu$
and the fourth represents the rotational time scale
$t_r = \lambda {\dot {\gamma }}^{-1}$
. We pick the same characteristic scales for physical quantities as previously used in flow studies for single component vesicles (Vlahovska & Gracia Reference Vlahovska and Gracia2007). All lengths are non-dimensionalized by the equivalent radius
$R$
, all times by the flow time scale
$t_{\dot {\gamma }}=\dot {\gamma }^{-1}$
, and all velocities by
$U_{flow} = R/t_{\dot {\gamma }} = R\dot {\gamma }$
. All pressures and viscous stresses are scaled by
$\eta \dot {\gamma }$
, whereas the membrane tractions
$\boldsymbol{f}^{\textit{mem}}$
are scaled by
$k_{0}/R^{2}$
.
Table 1 lists the set of physical parameters for this problem and their typical experimental values, while table 2 lists the dimensionless numbers for this problem. The first three dimensionless numbers are ones that are also found in the single component vesicle literature – the viscosity ratio parameter
$\lambda = (\mu ^{\textit{in}} + \mu ^{\textit{out}})/\mu ^{\textit{out}}$
between the inner and outer fluids, the capillary number
$\chi = \eta R^3{\dot {\gamma }}/k_0$
relating the viscous to bending forces on the membrane and the dimensionless excess area
$\varDelta = A/R^2 - 4\pi$
representing the floppiness of the vesicle. The remaining dimensionless numbers occur only in multicomponent vesicles. Of these, the most important ones are the dimensionless bending stiffness difference between the two phases
$\beta = k_1/k_0 = (\kappa _{lo} - \kappa _{ld})/(\kappa _{lo}+\kappa _{ld})$
, the Cahn number
$Cn = \gamma / (R\sqrt {\zeta _0} )$
(i.e. ratio of line tension energy to the energy scale of phase separation), the surface Péclet number
$\textit{Pe} = R^{2}\dot {\gamma }/\nu$
(i.e. ratio of coarsening time from diffusion to flow time) and the line tension parameter
$\alpha = k_0/\gamma ^2$
(ratio between bending energy to line tension energy). Note that the Cahn number can be re-expressed as the ratio between the interface thickness and vesicle size
$Cn = \epsilon ^{w\textit{idth}}/(\sqrt {2}R)$
. Furthermore, the average concentration
$q_{0}$
affects the position along the local energy landscape. It also affects stresses arising from the phase energy (2.6b).
Physical parameter ranges and orders of magnitude.
Dimensionless parameter ranges and orders of magnitude.
3. Solution for a nearly spherical vesicle
3.1. Perturbation expansion
We solve for the vesicle shape and composition as a function of time in the limits of small changes in deformation and concentration. We use
$\epsilon _1=\varDelta ^{1/2}$
and
$\epsilon _2 = q^+$
as perturbation variables with
$\epsilon _1 \sim \epsilon _2$
, where
$\varDelta$
is the excess area and
$q^+ = \sqrt {-{a}/{b}}$
is the phase split concentration (see table 2). This approximation (i.e. weak segregation approximation), has been used in multiple studies related to copolymer systems (Leibler Reference Leibler1980; Fredrickson & Helfand Reference Fredrickson and Helfand1987; Seul & Andelman Reference Seul and Andelman1995), and has been applied to multicomponent vesicles in the past (Taniguchi et al. Reference Taniguchi, Kawasaki, Andelman and Kawakatsu1994; Kumar, Gompper & Lipowsky Reference Kumar, Gompper and Lipowsky1999). We solve the dynamical equations in § 2.2 to leading order, keeping volume (
$V = 4\pi /3$
) and surface area (
$A = 4\pi + \varDelta$
) conserved to
$O(\epsilon _1^{2})$
and average order parameter
$(q_0 = A^{-1}\int q {\textrm d}A)$
conserved to
$O(\epsilon _2^2)$
.
Briefly, we expand the vesicle shape, surface tension and composition in terms of spherical harmonics. The non-dimensional vesicle quantities (radius
$r_{s}$
, membrane tension
$\sigma$
, and order parameter
$q$
) are written as
\begin{equation} Z(\theta , \phi ) = Z_0 + \sum _{l=1}^{\infty } \sum _{m=-l}^{+l} Z_{lm}Y_{lm}(\theta , \phi ) \end{equation}
where
$Z \equiv (r_{s},\sigma ,q)$
. In the above equation,
$Y_{lm}(\theta , \phi )$
are spherical harmonics (Appendix A), while
$Z_{lm}= (f_{lm},\sigma _{lm}, q_{lm})$
are coefficients for shape, surface tension and composition to be solved, with
$f_{lm}, \sigma _{lm} \sim O(\epsilon _1)$
and
$q_{lm} \sim O(\epsilon _2)$
, with
$\epsilon _1 \sim \epsilon _2$
. The isotropic quantities
$Z_0 = (r_0, \sigma _0, Q_0)$
are determined by applying the constraints of constant volume, constant area, and average order parameter. We obtain
$r_0 = 1 - {1}/{4\pi }\sum _{lm}f_{lm}f_{lm}^{*}$
and
$Q_0 = q_0 - {1}/{2\pi } \sum _{lm} q_{lm}f_{lm}^*$
. The isotropic surface tension
$\sigma _0$
is obtained implicitly from the expression
$\sum _{lm} (l + 2)(l-1) f_{lm}f^{*}_{lm} = \varDelta$
, which arises from conserving area to
$O(\epsilon _1^2)$
.
3.2. Solving Stokes equations
We solve the Stokes equations around the vesicle to determine how the shape and surface tension coefficients (
$f_{lm}, \sigma _{lm}$
) evolve over time. This section provides a high level overview of the steps, with algebraic details provided in Appendix B.
We Taylor expand the membrane tractions
$\boldsymbol{f}^{\textit{mem}}$
in (2.6) to
$O(\Delta ^{1/2})$
on the surface of a sphere. We then solve Stokes equations with this traction boundary condition (
$\chi \boldsymbol{f}^{\textit{hydro}} = \boldsymbol{f}^{\textit{mem}}$
at
$r = 1$
), along with boundary conditions of continuity of velocity across the interface (
$\boldsymbol{u}^{\textit{out}} = \boldsymbol{u}^{\textit{in}}$
at
$r = 1$
), membrane incompressibility (
$\boldsymbol{\nabla}_{\kern-1pt s} \boldsymbol{\cdot }\boldsymbol{u}^{\textit{in}} = 0$
at
$r = 1$
) and far-field velocity (
$\boldsymbol{u}^{\textit{out}} \rightarrow y \boldsymbol{\hat {x}}$
as
$r \rightarrow \infty$
). We note that this idea is commonplace in studies of single component vesicles (Misbah Reference Misbah2006; Vlahovska & Gracia Reference Vlahovska and Gracia2007) – the difference here is that we are examining a situation where
$\boldsymbol{f}^{\textit{mem}}$
takes a more complicated form (see § 2.2). Stokes flow is solved using fundamental basis sets for Lamb’s solution – i.e. a series solution using vector spherical harmonics. We use the notation in Vlahovska et al. (Reference Vlahovska, Loewenberg and Blawzdziewicz2005) and Bławzdziewicz et al. (Reference Bławzdziewicz, Vlahovska and Loewenberg2000), details are found in Appendix B.
Once one performs this procedure, one applies the kinematic boundary condition (2.8). Doing so yields the differential equation for the shape mode
$f_{lm}$
of the vesicle:
\begin{equation} \frac {{\textrm d} f_{lm}}{{\textrm d}t} = \frac {im}{2} f_{lm} + c_{lm2}^{+}. \end{equation}
The first term on the right-hand side comes from the rigid body rotation from shear flow, while the next term
$(c_{lm2}^{+})$
comes from the extensional deformation, which depends on the shape modes
$f_{lm}$
, concentration modes
$q_{lm}$
and isotropic surface tension
$\sigma _0$
. The final results are presented here,
\begin{equation} c_{lm2}^{+} = \lambda ^{-1}C_{lm} + \lambda ^{-1}\chi ^{-1} \left [ ( A_l + B_l \sigma _0) f_{lm} + D_l q_{lm} \right ] \end{equation}
where the coefficients (
$A_l, B_l, C_{lm}, D_l$
) are listed in (B10)–(B11) in Appendix B. These depend on the the mode number
$(l,m)$
as well as the dimensionless quantities discussed in table 2.
We note that the above expression (3.3) depends on the isotropic surface tension
$\sigma _0$
, which one obtains by applying constraint of constant excess area given by
${{\textrm d} \varDelta }/{{\textrm d}t} = 0$
, where the excess area is given by
$\varDelta = \sum _{lm} (l + 2)(l-1) f_{lm}f^{*}_{lm}$
. This constraint couples all modes in the differential equation for the shape in (3.2).
3.3. Solving the Cahn–Hilliard equation
To determine how the concentration modes
$q_{lm}$
evolve over time, we solve the Cahn–Hilliard equation (2.7). In the perturbative limits
$\varDelta \ll 1,$
$q^+ \ll 1$
, we Taylor expand all quantities to
$O(\varDelta ^{1/2})$
and
$O(q^+)$
on the unit sphere except the terms coming from the double-well potential (i.e.
${1}/{2}*aq^{2} + {1}/{4}*bq^{4}$
), where we keep all higher-order terms in the chemical potential. The reason why we perform this task is that such higher-order terms are needed to have a quartic free energy expression, which is necessary to have two-phase coexistence with a tie line. We note that similar approximations have been applied for equilibrium studies of multicomponent vesicles (Luo & Maibaum Reference Luo and Maibaum2020). In the cited study, the free energies are Taylor expanded to quadratic order in the shape and concentration modes except the double well term, where modes are kept to quartic order. This yields tractions and chemical potentials to the same order of accuracy as the current study.
The Cahn–Hilliard equation becomes
\begin{eqnarray} \frac {{\textrm d}q_{lm}}{{\textrm d}t} &= \dfrac {im}{2}q_{lm} - \dfrac {1}{\textit{Pe}} \left[ \varLambda _{lm} l(l+1) + Cn^{2}l^{2}(l+1)^{2}q_{lm} \right. \nonumber\\ &\quad \left. +\, 2\alpha Cn^{2}\beta (l-1)l(l+1)(l+2)f_{lm} \right], \end{eqnarray}
where
$\varLambda _{lm} = \int _{\varOmega } (aq+bq^{3} )Y_{lm}{\textrm d}\varOmega$
is the chemical potential from the double-well potential, projected onto spherical harmonics by integrating over a unit sphere. This term is weakly nonlinear as discussed above.
3.4. Numerical procedure
The final equations we solve are (3.2) for the vesicle shape along with conservation of area, as well as (3.4) for the concentration. These are coupled, nonlinear ordinary differential equations (ODEs).
When solving these equations, we decompose the shape and concentration modes into real and imaginary parts
\begin{equation} f_{lm} = f^{\prime}_{lm}+if^{\prime \prime}_{lm} ,\qquad q_{lm} = q^{\prime}_{lm}+iq^{\prime \prime}_{lm}, \end{equation}
and perform an ODE solve for the components (
$f^{\prime}_{lm}, f^{\prime \prime}_{lm}, q^{\prime}_{lm}, q^{\prime \prime}_{lm}$
). We only solve for components with
$m \geq 0$
since for the radius to be real,
$f^{\prime}_{l-m} = (-1)^m f^{\prime}_{lm}$
and
$f^{\prime \prime}_{l-m} = (-1)^{m+1} f^{\prime \prime}_{lm}$
, with similar relations holding for the concentration. Since we need to ensure that the constraint of area conservation is satisfied
$A = 4\pi + \varDelta$
, utilizing regular solvers is not feasible. We have developed a predictor–corrector scheme to solve the ODEs as mentioned in Appendix C. After solving for the spectral coefficients, we reconstruct the vesicle shape and concentration fields using the spherical harmonic series in (3.1). A benchmark is provided in Appendix D for examining coarsening on a rigid sphere in the absence of flow. We have also benchmarked our results against those of previous studies pertaining to single component vesicles (Vlahovska & Gracia Reference Vlahovska and Gracia2007).
For most of the simulations we perform, we choose a cutoff mode number between
$l_{\textit{max}} = 6$
and
$l_{\textit{max}} = 10$
when computing the order parameter. We find that beyond this number, the dominant modes
$0 \leq l \leq 4$
for shape and order parameter change less than 2 %, and thus the trends discussed in the results section do not change appreciably. Figure 2 shows an example of a convergence test. Here, we computed the order parameter
$q(\theta ,\phi , t)$
using different cutoff modes, and compared the results with a simulation using a high cutoff mode. The relative error we computed was
$e = \int_S (\bar {q}_{\textit{approx}} - \bar {q}_{\textit{true}}) {\textrm{d}}A / \int_S \bar {q}_{\textit{true}}{\textrm{d}}A$
, where
$\bar {q}_{\textit{true}}$
is the time-averaged value at
$(\theta ,\phi )$
using the highest cutoff mode and
$\bar {q}_{\textit{approx}}$
is the time-averaged value at
$(\theta , \phi )$
using the lower cutoff mode. For small values of
$Cn$
(
$Cn \ll 1$
), we expect that a spectral method may lead to inaccuracies as one will encounter the sharp interface limit that is hard to resolve owing to the Gibbs phenomenon.
Error convergence plots for
$Cn=0.5$
. The parameters for the simulations are
$\alpha =1,\beta =0.1,\varDelta =0.01,\chi =0.5,a=-0.1,b=1,q_0=0,\textit{Pe}=5$
,
$\lambda = 2$
. The
$y$
-axis represents the ratio of the difference in the order parameter over the vesicle surface divided by the average magnitude of the order parameter of the case with most modes
$l_{\textit{max}} = 11$
.
4. Multicomponent vesicles,
$\boldsymbol{Pe} \boldsymbol{\gg 1}$
4.1. Overview and initial conditions
This section will examine multicomponent vesicles in the limit of large Péclet number (
$\textit{Pe} \gg 1$
), where the flow time scale is much faster than the coarsening time scale. In this regime, coarsening effects are unimportant, and the dynamics is determined by the inhomogeneous bending rigidity set by the initial condition. This regime has previously been studied for 2-D vesicles (Gera et al. Reference Gera, Salac and Spagnolie2022) where motions such as tumbling and swinging have been observed. However, we are unaware of similar studies for 3-D systems.
One question that naturally arises is whether the 3-D geometry will alter the dynamics significantly compared with the 2-D case. To answer this question, we will structure this section as follows. We will first examine initial conditions for concentration (
$q_{lm}$
) and shape (
$f_{lm}$
) where the
$l = 2$
modes are excited, giving rise to a prolate-like shape with two lobes of a stiff (
$L_o$
) phase. These were argued to be the dominant modes in the 2-D study. We will first excite the
$l=2$
modes that are symmetric about the flow-shear plane (
$l = 2, m = 0, \pm 2$
), replicating conditions similar to a 2-D study, and then excite the
$l=2$
modes that do not obey this symmetry – coined out-of-plane modes (see figures 3a and 3b). We will conclude by examining initial conditions in which many
$l\neq 2$
modes are excited for concentration (figure 3
c).
Initial conditions for
$\textit{Pe} \gg 1$
simulations. (a) In-plane
$l = 2$
modes,
$f^{\prime}_{22}=f^{\prime \prime}_{22}=f_{20}=\sqrt {0.1\varDelta },q^{\prime}_{22}=\sqrt {\varDelta /8}$
; (b) out-of-plane
$l=2$
modes,
$f^{\prime}_{22}=f^{\prime \prime}_{22}=f_{20}=\sqrt {0.1\varDelta };q^{\prime}_{22}=q^{\prime}_{20}=,q^{\prime}_{21}=\sqrt {\varDelta /8}$
. (c) One lobe initial condition with
$l \neq 2$
concentration modes,
$f^{\prime}_{22}(0)=f^{\prime \prime}_{22}(0)=f_{20}(0) = \sqrt {{\varDelta }/{10}}$
, and
$q^{\prime}_{lm}(0) = \sqrt {{\varDelta }/{8}},q^{\prime \prime}_{lm}(0)=0;l=1,2,3,m=0,\ldots,l$
.
Overall, when the initial shape and concentration profiles obey in-plane symmetry, one sees in-plane tumbling and swinging behaviour very similar to the 2-D studies. When out-of-plane modes are excited, one will also see tumbling and swinging, but the swinging can exhibit weak, out-of-plane oscillations, while the tumbling can occur out of plane. If the initial concentration profile is not symmetric with respect to the flow-gradient plane, the profile will remain non-symmetric and will retain the same number of lobes as the initial condition (thus, the concentration solution is not limited to the
$l=2$
subspace). Details are provided below.
4.2. Exciting
$l = 2$
modes in-plane
This section explores the effects of exciting the
$l=2$
modes in-plane for shape and concentration. This initial condition corresponds to the situation mentioned in figure 3(a). In this case, an analytical theory for in-plane dynamics can be developed, similar to previously developed 2-D theories (Gera et al. Reference Gera, Salac and Spagnolie2022).
When
$\textit{Pe} \to \infty$
, the diffusive coarsening effects are negligible and the phospholipids are convected by the background shear flow. Setting
$\textit{Pe}^{-1} = 0$
in (3.4) yields
\begin{equation} \frac {{\textrm d} q_{lm}}{{\textrm d}t} = \frac {im}{2} q_{lm}, \end{equation}
and hence the solution is oscillatory:
$q_{lm}(t) = q_{lm}(0) \exp (imt/2)$
. The shape of the vesicle is then given by (3.2) along with conservation of area constraint, using the above expression for
$q_{lm}$
. We will solve for the shape and concentration modes for
$l=2, m = \pm 2$
. We find that the
$l = 2, m=0$
shape mode (
$f_{20}$
) quickly decays to zero.
Since the dynamics are in plane, we will examine the equatorial plane of the vesicle (
$\theta = \pi /2$
) and decompose the shape and concentration as follows (using the notation of Gera et al. Reference Gera, Salac and Spagnolie2022):
\begin{align}r_s(\theta = \pi /2, \phi , t) &= 1 + {\epsilon _1} a_2(t) \cos (2\phi ) + {\epsilon _1} b_2(t) \sin (2\phi ),\end{align}
\begin{align}q(\theta = \pi /2, \phi , t) &= {\epsilon _2} \cos (2\phi + t),\end{align}
where
$a_2$
and
$b_2$
are time dependent coefficients to be determined, and
$\epsilon _1 = \varDelta ^{1/2}$
and
$\epsilon _2 = q^+$
are the small parameters in this problem. In terms of the coefficients (
$f_{lm}, q_{lm}$
) discussed in the paper,
$f_{2\,\pm \,2} = 2 {\epsilon _1} \sqrt {{2\pi }/{15}} ( a_2 \mp i b_2 )$
and
$q_{2\,\pm\, 2} = 2 {\epsilon _2} \sqrt {{2\pi }/{15}} \exp (\pm it)$
. If we define a modified capillary number as
$C = \chi /{\epsilon _1}$
and modified time
$\tau = t/{\epsilon _1}$
, we find the leading-order shape dynamics in (3.2) for
$C \sim O(1)$
,
$\tau \sim O(1)$
and
$\epsilon _2/\epsilon _1 \sim O(1)$
to be
\begin{align} \frac {{\textrm d} a_2}{{\textrm d} \tau } &= -\varTheta \left (a_2, b_2, {\epsilon _1 \tau , \frac {\epsilon _2}{\epsilon _1}}, C, \lambda , \beta \right ) b_2, \nonumber\\ \frac {{\textrm d} b_2}{{\textrm d} \tau } &= \varTheta \left (a_2, b_2, {\epsilon _1 \tau , \frac {\epsilon _2}{\epsilon _1}}, C, \lambda , \beta \right ) a_2 . \end{align}
These sets of ODEs take the same form as the 2-D case, except that the expression
$\varTheta$
is different due to the geometry of the problem. Appendix E shows the algebra to obtain
$\varTheta$
– here we state the final results. For the 3-D system, we obtain
\begin{equation} \varTheta _{3D} = \frac {512 \pi }{5(9 + 23\lambda )} \left [ \frac {5}{8} a_2 + {\frac {\epsilon _2}{\epsilon _1}\beta C^{-1} \left ( a_2 \sin (\epsilon _1 \tau ) + b_2 \cos (\epsilon _1 \tau ) \right )} \right ]\!, \end{equation}
while for the 2-D system, we obtain
\begin{equation} \varTheta _{2D} = \frac {3 \pi }{2\lambda } \left [ a_2 + {\frac {\epsilon _2}{\epsilon _1} \beta C^{-1} \left ( a_2 \sin (\epsilon _1 \tau ) + b_2 \cos (\epsilon _1 \tau ) \right )} \right ]\!, \end{equation}
where in the 2-D case the small parameter
$\epsilon _1$
is related to the excess length as
$\epsilon _1 = \varDelta _{2D}^{1/2}$
(see Appendix E). The expressions for
$\varTheta$
depend on the viscosity ratio
$\lambda$
, reduced area
$\epsilon _1^2 = \varDelta$
and a lumped bending stiffness parameter
$S=({\epsilon _2}/{\epsilon _1})\beta C^{-1} = \beta q^{+}/\chi$
. Results are discussed below.
The system of equations (4.3) admit a periodic solution for the shape modes
$a_2(t)$
and
$b_2(t)$
. Two different behaviours are seen:
-
(i) swinging – after an initial transient, the vesicle undergoes weak oscillations in its inclination angle
$\phi _i$
about a fixed value, while the concentration field performs a full rotation along the membrane (for visualization, see figure 4); -
(ii) tumbling – the vesicle shape and concentration field undergo in-plane rigid body rotation (for visualization, see figure 5).
Swinging visualization in flow-gradient plane for
$\textit{Pe} \gg 1$
. The initial condition is figure 3(a), and the parameters are
$\textit{Pe}=10^4,a=-0.1,b=1,Cn=0.4,\alpha =1,\beta =0.5,\varDelta =0.1,\chi =0.6,q_0=0, \lambda =2$
. The blue phase represents the softer phase and the yellow phase represents the stiffer phase. The black double arrow represents oscillations about a fixed axis (dashed line) after an initial transient. The blue arrow represents the motion of the phases around the vesicle. See Supplementary movie (Video1.mp4) for colourbar details.
Tumbling visualization in the flow-gradient plane for
$\textit{Pe} \gg 1$
. The initial condition is figure 3(a), and the parameters are
$\textit{Pe}=10^4,a=-0.1,b=1,Cn=0.4,\alpha =1,\beta =0.5,\varDelta =0.1,\chi =0.01,q_0\!=\!0,\lambda\! =\! 2$
. The blue phase represents the softer phase and the yellow phase represents the stiffer phase. The black and blue arrows represent rigid body motion of the shape and concentration field, while the dashed line is the initial orientation. See Supplementary movie (Video2.mp4) for colourbar details.
For a fixed viscosity ratio
$\lambda$
and excess area
$\varDelta$
, the transition between the two behaviours is determined by the dimensionless quantity
$S=({\epsilon _2}/{\epsilon _1}) \beta C^{-1} = \beta q^+/\chi$
. Below a critical value of this parameter, swinging occurs, while tumbling occurs above the critical value. Figure 6 illustrates a typical phase diagram. We see that the critical conditions exhibit similar trends for 3-D and 2-D vesicles, although their quantitative values might be different (Gera et al. Reference Gera, Salac and Spagnolie2022). The explanation for these trends is given by the authors as follows. (i) At low background shear rates, the vesicle elongates in the direction of the softer phase in a quasisteady manner. These regions have a higher curvature. In order to match the phase to the curvature, the vesicle undergoes rigid body tumbling. (ii) On the other hand, when the shear rates are high, the viscous stresses dominate over the elastic stresses thereby pushing the stiffer phase past the high curvature region. In order to account for this phase-shape mismatch, the vesicle undergoes a swinging motion.
Phase diagram for
$\textit{Pe} \to \infty$
limit at
$\varDelta =0.1$
and
$\lambda = 2,q^{+}=0.1$
.
Figure 7 shows a plot for
$a_2(t)$
and
$b_2(t)$
for swinging and tumbling vesicles, comparing the 3-D and 2-D theories in the limit
$\textit{Pe} \rightarrow \infty$
. As mentioned, the agreement is not exact due to multiple prefactors entering the 3-D analysis compared with the 2-D case. However, the agreement for tumbling is nearly exact in both cases. The discrepancy primarily arises due to the differences in expressions for
$\varTheta _{2D}$
and
$\varTheta _{3D}$
.
Here
$\textit{Pe} \rightarrow \infty$
results comparison with Gera et al. (Reference Gera, Salac and Spagnolie2022) at
$\varDelta =0.0001$
by solving (4.3) and assuming
$\epsilon _2 =\epsilon _1=\varDelta ^{1/2}$
.
Lastly, in Appendix F, we show a comparison between the
$\textit{Pe} \rightarrow \infty$
theory (4.3) and full numerical simulations (solving (3.2) and (3.4)). We see that the analytical theory matches well with numerics for
$\textit{Pe} \sim O(100)$
or above for both tumbling and swinging. We have also seen that shape deformations
$f_{20}$
in simulations rapidly decay to zero, indicating that deformations in this mode do not appear to affect the theory quantitatively.
4.3. Effects of higher modes being excited
In this section, we look at effects of exciting
$q_{lm}$
modes that are not constrained to
$l=2$
. The initial condition considered in the following simulations is
$f^{\prime}_{22}(0)=f^{\prime \prime}_{22}(0)=f_{20}(0) = \sqrt {{\varDelta }/{10}}$
, and
$q^{\prime}_{lm}(0) =\sqrt {{\varDelta }/{8}},q^{\prime \prime}_{lm}(0)=0;l=1,2,3,m=0,\ldots,l$
. On superimposing these modes, the resultant structure is one lobe containing the stiffer phase and the remaining being primarily the softer phase (see figure 3
c).
Similar to the previous case, swinging is seen where the concentration performs a full rotation along the membrane while the inclination angles exhibit weak oscillations in the shear-flow and the flow–vorticity planes (figure 8). The modes are plotted in figure 9 to understand the numerical details. We calculate the inclination angles (
$\theta _i, \phi _i$
) with respect to the shear-flow plane and shear-vorticity plane using the moment of inertia tensor (Ramanujan & Pozrikidis Reference Ramanujan and Pozrikidis1998). The vesicle maintains its one lobe concentration profile during its motion, indicating that the spatial concentration pattern is set by the initial condition.
Swinging visualization when higher-order
$q_{lm}$
modes are excited for
$\textit{Pe} \gg 1$
. The initial condition is figure 3(c), and the parameters for simulations are
$\beta =0.5,\alpha =1,a=-0.1,b=1,\textit{Pe}=10^4,\varDelta =0.1,\chi =0.6,Cn=0.4,\lambda =2,q_0=0$
. The blue phase represents the softer phase and the yellow phase represents the stiffer phase. The black double arrow represents oscillations about a fixed axis (dashed line) after an initial transient. The blue arrow represents the motion of the phases around the vesicle. See Supplementary movie (Video7.mp4,Video8.mp4) for colourbar details.
Modes and inclination angles for figure 8 (swinging for
$\textit{Pe} \gg 1$
).
In the second type of motion (tumbling), as seen in figure 10, the oscillations in the shear-flow plane are significant (rigid-body rotations), whereas the oscillations in the vorticity-flow plane are small. We plot the modes in figure 11. We can observe that the
$f_{3m}$
deformations are an order of magnitude smaller than
$f_{2m}$
modes which is expected due to the nature of the background flow. This vesicle behaves almost like an oblate-spheroid rotating in flow.
Tumbling visualization when higher-order
$q_{lm}$
modes are excited for
$\textit{Pe} \gg 1$
. The initial condition is figure 3(c), and the parameters for simulations are
$\beta =0.5,\alpha =1,a=-0.1,b=1,\textit{Pe}=10^4,\varDelta =0.1,\chi =0.01,Cn=0.4,q_0=0,\lambda =2$
. The blue phase represents the softer phase and the yellow phase represents the stiffer phase. The black and blue arrows represent rigid body motion of the shape and concentration field, and the dashed line is the initial orientation. See Supplementary movies (Video9.mp4,Video10.mp4) for colourbar details.
Modes and inclination angles for figure 10 (tumbling for
$\textit{Pe} \gg 1$
).
We also performed multiple simulations with initial conditions same as those in figure 3(b). The primary observation is that the nature of motion remains similar to figures 8 and 10, albeit with a different concentration topology. In these simulations, the order parameter fields look like figure 3(b). The movies for such motions can be seen in Video3.mp4 (swinging), Video4.mp4 (swinging), Video5.mp4 (tumbling) and Video6.mp4 (tumbling).
5. Multicomponent vesicles, intermediate
$\boldsymbol{Pe}$
In this section, we focus on vesicles undergoing motion for a
$\textit{Pe}\sim O(1{-}10)$
. Under this regime, the flow and coarsening time scales are comparable, which give rise to a wider range of vesicle shape behaviours compared with the previous section. Additional parameters to consider (beyond the ones discussed previously) are (i) line tension parameter
$\alpha = k_0/\gamma ^2$
, which indicates the relative bending to line tension energy, and (ii) Cahn number
$Cn = \epsilon ^{w\textit{idth}}/(\sqrt {2}R)$
, which is related to the interface width compared with the vesicle size. This is also related to line tension. Just like in the previous section, we will divide our analysis based on the initial conditions of the simulations – (i) in plane dynamics (flow-shear), and (ii) out-of-plane dynamics.
5.1. In-plane dynamics
Below summarizes the behaviour observed in the long time limit when one excites the
$l=2$
shape and concentration modes in-plane. The initial condition is given in figure 3(a).
-
(i) Tank-treading. Here, both the vesicle’s shape and concentration field are stationary. The vesicle is ellipsodal and fixed at a steady inclination angle
$\phi _i$
, and the concentration field is frozen (figure 12). This behaviour can only be found when coarsening effects are present, since the concentration field in the
$\textit{Pe} \gg 1$
limit only admits periodic solutions (see (4.1)). -
(ii) Tumbling. The vesicle’s shape and concentration field exhibit rigid body rotation (figure 13). One can observe that the
$f_{20}$
mode decays rapidly, thereby leaving purely in-plane rotations. -
(iii) Vacillated – breathing/trembling. The vesicle performs rotations in the flow-gradient plane while there are significant deformations in the flow–vorticity plane represented by the
$f_{20}$
mode (figure 14). A detailed understanding of this motion is given in Misbah (Reference Misbah2006).
Tank-treading for
$\textit{Pe} \sim O(1)$
. The initial condition is figure 3(a), and the parameters are
$\chi = 0.1,\alpha =1,Cn=1,\beta = 0.1,\lambda =2,\textit{Pe}=5,\varDelta =0.1$
,
$a =-0.1, b=1, q_0=0$
. The blue phase represents the softer phase, and the yellow phase represents the stiffer phase. The inset snapshot represents the frozen state and inclination angle of the vesicle. See the Supplementary movie (Video11.mp4) for colourbar information.
Tumbling for
$\textit{Pe} \sim O(1)$
. The initial condition is figure 3(a), and the parameters are
$\alpha =1,\beta =0.5,Cn=0.4,\chi =0.01,\textit{Pe}=5,\lambda =2,\varDelta =0.1$
,
$a=-0.1, b=1,q_0=0$
. The colour bar represents the order parameter
$q$
. The blue phase represents the softer phase and the yellow phase represents the stiffer phase. The blue and black arrows represent rigid-body motion of lipids and shape, and the dashed line is the initial orientation. See Supplementary movie (Video12.mp4) for visualization.
Breathing for
$\textit{Pe} \sim O(1)$
. The initial condition is figure 3(a), and the parameters are
$\chi = 0.01,\alpha =1,Cn=0.8,\beta = 0.1,\lambda =10,\textit{Pe}=5,\varDelta =0.1$
,
$a=-0.1, q_0=0, b=1$
. The inset snapshot represents the visualization of the vesicle undergoing breathing motion. The blue phase represents the softer phase and the yellow phase represents the stiffer phase. See Supplementary movies (Video13.mp4,Video14.mp4) for visualization and colourbar information.
Figures 12–14 show snapshots of these behaviours, as well as typical plots for the shape modes
$f_{lm}$
, concentration modes
$q_{lm}$
and in-plane inclination angle
$\phi _i$
. Movies are also provided in the Supplementary material.
Occasionally, we have observed situations where the vesicle undergoes tumbling motion before the motion dampens out at very long times to tank-treading behaviour. This has been noted in figure 15. We believe this is not a numerical instability but arises from the nonlinear mode-mixing term
$bq^{3}$
in the Cahn–Hilliard equation. It is important therefore to understand what is the time scale at which these motions are observed experimentally. In figure 15, we see that the tumbling motion persists up to
$t\sim 10$
, which translates to a physical time of
$10{\dot {\gamma }}^{-1}$
that could range from
$O(1)-O(10)$
${\textrm{s}}$
. Shear flow experiments are usually limited by issues like photobleaching or the escape of the vesicle from the field-of-view. Hence, we believe in order to characterize the motion, the time scale of observation should be reported for a phase diagram. For the other dynamical regimes listed above, we do not observe a change of behaviour at very long times
$(t \sim O(100))$
.
An example of tumbling motion dampening to give rise to tank-treading. The parameters are
$\lambda =2,\beta =0.25,Cn=0.5,\textit{Pe}=5,\varDelta =0.1,\chi =0.01,a=-0.1,b=1,q_0=0, \alpha =1$
.
5.2. Out-of-plane dynamics
In this section, we inspect the effects of exciting out-of-plane, higher-order modes not constrained to the
$l=2$
space. The initial condition is in figure 3(c). In the following simulations, unless specified, we used a cutoff mode
$l_{\textit{max}} = 8$
when solving the shape and order parameter. The following behaviours were observed.
5.2.1. Tank-treading (stationary)
In some simulations, we see the vesicle is frozen with a steady shape and concentration profile. This typically occurs when the interface width between phases is diffuse (
$Cn\sim 1$
), where the vesicle behaves similarly to a single component vesicle with matched viscosity
$\lambda =2$
– stationary tank-treading motion. The steady inclination in-plane angle
$\phi _{i}$
depends on the excess area and has been discussed previously (Vlahovska & Gracia Reference Vlahovska and Gracia2007). The out-of-plane inclination angle
$\theta _{i}$
is
$\pi /2$
, which indicates that the vesicle remains in the flow-shear (gradient) plane. We also note that the initial condition that primarily consisted of one yellow (stiff) lobe develops two lobes and stays stationary (indicated by
$q^{\prime}_{22},q^{\prime \prime}_{22}$
). Essentially, this motion is very similar to the in-plane counterpart discussed in § 5.1. See Supplementary movie Video15.mp4 for more details.
5.2.2. Tumbling
When the capillary number
$\chi$
is low, the vesicle is stiffer and tumbles under the influence of the background shear flow. The vesicle preserves its one-yellow-lobe configuration. The
$f_{2m}$
modes exhibit oscillatory motion. This motion is similar to the high
$\textit{Pe}$
limit tumbling (see figure 10). See Supplementary movies Video16.mp4 and Video17.mp4 for a visualization.
5.2.3. Swinging/phase-treading
When the capillary number
$\chi$
is higher, the vesicle undergoes swinging motion where the concentration rotates along the membrane and the shape exhibits weak oscillations about a fixed orientation. Similar to the high
$\textit{Pe}$
limit, the shape oscillations are weakly out of plane. However, unlike the high
$\textit{Pe}$
limit where the topology of the concentration profile is set by the initial condition, here the topology of the profile can change due to coarsening effects. For example, for the initial condition discussed in figure 3(c) where the vesicle has one single stiff (
$L_o$
) domain, coarsening during swinging gives rise to two stiff phases treading along the surface (figures 16 and 17). However, this effect is highly dependent on flow strength, bending stiffness and line tension. If for example we keep the same conditions as figure 16, except we increase the capillary number (increase
$\chi$
) and increase line tension (decrease
$\alpha$
), one can see a situation where the one lobe initial condition persists and gives rise to
$l=1$
mode swinging in the subsequent results (figures 18 and 19). In this motion, the
$q_{11},q_{1-1}$
modes persist while all other higher modes
$q_{lm}$
decay. The vesicle behaves like a Janus particle with divided phases, similar to what has been noted in previous studies (Gera & Salac Reference Gera and Salac2018a).
Phase-treading/swinging snapshots when higher-order
$q_{lm}$
modes are excited at
$\textit{Pe} \sim O(1)$
. The initial condition is figure 3(c), and parameters are
$\beta =0.5,\alpha =3,a=-0.1,b=1,\textit{Pe}=5,\varDelta =0.1,\chi =0.3,Cn=0.2,\lambda =2,q_0=0$
. The blue arrow represents the lipid motion whereas the double-headed black arrow represents the vesicle shape oscillations about a fixed angle (black dashed line). The concentration profile coarsens from one lobe to two lobes. See Supplementary movies (Video18.mp4, Video19.mp4) for visualization and colourbar information.
Modes and inclination angles for figure 16 (phase-treading at
$\textit{Pe} \sim O(1)$
with topology change).
Phase-treading/swinging snapshots at
$\textit{Pe} \sim O(1)$
when higher-order
$q_{lm}$
modes are excited with
$l=1$
modes driving the lipid motion. The initial condition is figure 3(c), and the parameters are
$\beta =0.5,\alpha =1,a=-0.1,b=1,\textit{Pe}=5,\varDelta =0.1,\chi =0.6,Cn=0.2,\lambda =2$
and
$q_0=0$
. The blue arrow represents the lipid motion whereas the double-headed black arrow represents oscillations about a fixed angle (black dashed line). See Supplementary movies (Video20.mp4, Video21.mp4) for visualization and colourbar information.
Modes and inclination angles for figure 18 (phase-treading at
$\textit{Pe} \sim O(1)$
without topology change).
Phase diagrams for
$\alpha$
versus
$\chi$
at
$\varDelta =0.1,\textit{Pe}=5,\beta =0.5,\lambda =2,a=-0.1,b=1,q_0=0$
. In these simulations, the initial condition is figure 3(c).
5.2.4. Phase diagrams
This section provides phase diagrams to quantify general trends for vesicle dynamics. Since there is a large region of phase space to explore, we will isolate the effect of a few variables mentioned in § 2.3.
The first set of variables we will explore are the capillary number (
$\chi$
) and line tension parameter (
$\alpha$
). We will choose a moderately deflated vesicle (
$\varDelta = 0.1$
) with matched interior and exterior fluid viscosity (
$\lambda = 2$
) and zero average order parameter (
$q_0 = 0$
), and choose
$O(1)$
values for the other parameters:
$\textit{Pe} = 5, \beta = 0.5$
. Figure 20 shows the types of motion observed by the vesicle for different values of capillary number
$\chi$
and line tension parameter
$\alpha$
, when the initial condition is the same as that of figure 3(c). Two different interface widths are examined:
$Cn = 0.2$
(more sharp) and
$Cn = 0.55$
(less sharp). The regimes are reported over time window
$t \approx 50$
.
We will now explain the reasons behind these phase diagrams. In the first case where we see the relatively sharper interface
$Cn=0.2$
(figure 20
a), we see that at small capillary numbers (
$\chi$
), the vesicle exhibits tumbling motion. The rationale for this is that a low value of
$\chi$
increases
$f_{lm}/\chi$
in the shape equations (3.2)–(3.3). This effect causes the vesicle to elongate in the softer phase direction, after which it is rotated by the vorticity of the background flow, leading to tumbling motion, very similar to what is seen in the high
$\textit{Pe}$
limit (see § 4). When we increase the
$\chi$
value, depending on how stiff the vesicle is (given by
$\alpha$
) we can get swinging/phase-treading motion. The
$({im}/{2})q_{lm}$
term in the Cahn–Hilliard equation (3.4) drives the oscillatory behaviour of
$q_{lm}$
. This translates to the shape equation (3.2)–(3.3) showing oscillatory motion due to the
$q_{lm}$
term, but the
$\chi ^{-1}$
prefactor in (3.3) dampens these oscillations leading to swinging/phase-treading motion. We have further subdivided the phase-treading/swinging motion into
$l=1$
and
$l=2$
driven, akin to figures 16 and 18.
When the vesicle interface becomes more diffuse (
$Cn$
increases), we see that the vesicle shows behaviour more similar to single-component vesicles. This is seen in figure 20(b) where the
$Cn=0.55$
. At small
$\chi$
values, the term
$\chi /\varDelta ^{1/2}$
is small thereby driving tumbling motion as stated before. As
$\chi$
increases, the
$f_{lm}$
oscillations are dampened out in the shape equations (3.2)–(3.3) and the
$q_{lm}$
oscillations in the Cahn–Hilliard equation (3.4) are dampened due to the line tension penalty term
$Cn^{2}l^{2}(l+1)^{2}q_{lm}$
coming into play. As this
$Cn$
increases further, the vesicle will largely show tank-treading behaviour, which is expected for single component vesicles with a matched viscosity (
$\lambda =2)$
.
6. Effect of asymmetric distributions
In the previous sections, we presented results based on
$q_{0}=0$
which indicates an effective
$50:50$
phase distribution of the two separated phases. However, it is possible that we might encounter cases where
$q_0 \neq 0$
. In such cases, our theory can help predict the dynamics in limits where the
$q_{0}$
value is not large enough to lie outside the phase-split region. On inspecting the shape equations (B10), we see that the shape-concentration coupling term
$D_l$
could make a drift caused by
$q_0$
. One can observe from (B10) that for
$l=1$
modes, the shape coefficients take the following form:
\begin{equation} \frac {{\textrm d}f_{1m}}{{\textrm d}t} = \frac {imf_{1m}}{2} -\chi ^{-1}\left [\frac {4\mu _0}{9\alpha Cn^{2}}\right ]q_{1m}, \end{equation}
where
$\mu _{0} = aq_{0}+bq_{0}^3$
This indicates that any non-zero initial condition in
$q_{1m}$
would induce a drift in
$f_{1m}$
, which would lead to translational motion – i.e. centre of mass tracking.
Lateral drift caused by asymmetric phase split (
$q \neq 0$
). In both cases, the initial condition is
$f^{\prime}_{22}(0)=f^{\prime \prime}_{22}(0)=f_{20}(0) = \sqrt {{\varDelta }/{10}}$
, and
$q^{\prime}_{10}(0) = q^{\prime}_{11}(0) = q^{\prime}_{22}=\sqrt {{\varDelta }/{8}}$
. The parameters for the simulations are
$\textit{Pe}=5,\lambda =2,a=-0.1,b=1,\varDelta =0.1,\beta =0.5,\chi =0.6,\alpha =1,Cn=0.4$
.
This phenomenon can be observed in figure 21. When one performs simulations at the same initial conditions but examines two different average concentrations
$q_{0}=0.01$
and
$q_0 = 0$
, one observes that a non-zero
$q_0$
generates a oscillatory drift in the
$f^{\prime}_{11},f^{\prime \prime}_{11}$
modes and a non-zero drift in
$f_{10}$
. This indicates that the centre of mass of the vesicle performs sinusoidal oscillations about a fixed point in space. This lateral drift has been observed in previous studies as well (Gera & Salac Reference Gera and Salac2018b). Lastly, we note that although the mechanisms are different, a similar oscillatory translational motion has been observed in single component vesicles in quadratic flows with confinement, and has been coined ‘snaking’ (Lyu et al. Reference Lyu, Chen, Farutin, Jaeger, Misbah and Leonetti2023).
7. Conclusion
In this study, we explored the effect of shear flow around a multicomponent vesicle containing a ternary mixture of phospholipids and cholesterol molecules within the bilayer. We analysed the small-deformation, weakly segregated limit and found that multicomponent vesicles exhibit very different dynamics compared with their single component counterparts. Results were examined in two distinct limits: one where the coarsening dynamics on the surface are unimportant (
$\textit{Pe} \gg 1$
), and another where the coarsening time scale is comparable to the flow time scale (
$\textit{Pe} \sim O(1)$
).
When
$\textit{Pe} \gg 1$
, the phases are advected by the shear flow over the surface. The strength of the coupling with the shape deformations, given by
$S = \beta q^{+}/\chi$
, determines the nature of this motion – swinging or tumbling. The initial conditions heavily dictate the dynamics in this system. In cases where purely the
$l=2$
in-plane modes are excited, the results show qualitative agreement with previous 2-D studies on such vesicles at the
$\textit{Pe} \rightarrow \infty$
limit (Gera et al. Reference Gera, Salac and Spagnolie2022). However, when out-of-plane modes are excited, it is possible to observe weakly out-of-plane tumbling and swinging motion. The concentration profile will perform periodic dynamics and retain its initial topology over time. Thus, if it starts asymmetric with respect to the flow-gradient plane, it will remain asymmetric over time.
In the intermediate
$\textit{Pe}$
limit (
$\textit{Pe} \sim O(1)-O(10)$
), we observe variety of motions based on the initial conditions and physical parameters such as tumbling, breathing, swinging/phase-treading and tank-treading. The biggest difference between
$\textit{Pe} \sim O(1)$
and
$\textit{Pe} \gg 1$
dynamics is that coarsening can occur on the membrane, which can give rise to stationary shapes and concentration profiles like tank-treading (which is not admissible in the
$\textit{Pe} \rightarrow \infty$
limit), and can alter the topology (i.e. number of lobes) in the concentration profile during swinging or tank-treading. We provided a general overview of the dynamics and illustrated phase diagrams for the different dynamical regimes. Generally, tumbling is favoured at small capillary numbers. At higher capillary numbers, phase-treading or tank-treading can occur, with the latter favoured at lower line tensions (e.g. more diffuse interface between phases). Lastly, we discussed the role of asymmetric phase split on the dynamics (
$q \neq 0$
), which can give rise to effects like oscillatory translational motion.
While the weak-segregation, small-deformation limit eases up the mathematical rigour largely due to the exclusion of nonlinear stresses, we believe it provides a great starting point for further semianalytical treatments of these systems using higher-order theories.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.10632.
Funding
The authors would like to acknowledge support from National Science Foundation (grant 2147559-CBET).
Declaration of interests
The authors report no conflicts of interest.
Appendix A. Spherical harmonics
We define spherical harmonics
$Y_{lm}$
using the following convention (Jones Reference Jones1985):
\begin{equation} Y_{lm}(\theta , \phi ) = \sqrt {\frac {(2l+1)(l-m)!}{4\pi (l+m)!}}(-1)^{m}P_{l}^{m}(\cos \theta )e^{im\phi }, \end{equation}
where
$(\theta ,\phi )$
are the polar and azimuthal angles, and
$P_l^m$
are associated Legendre polynomials. The inner products of these spherical harmonics are orthonormal – i.e.
$\int Y_{lm} Y_{l^{\prime}m^{\prime}}^* {\textrm{d}} \varOmega = \varDelta _{ll^{\prime}} \varDelta _{mm^{\prime}}$
, where the integral is over the unit sphere and
$*$
represents complex conjugate. Furthermore, the following relationship holds:
$Y_{l -m} = (-1)^m Y_{lm}^*$
.
We can also define vector spherical harmonics as follows:
\begin{equation} \boldsymbol{y}_{lm0} = \left [l(l+1)\right ]^{-1/2}r\boldsymbol{\nabla }Y_{lm}, \qquad \boldsymbol{y}_{lm1} = i\boldsymbol{\hat {r}}\times \boldsymbol{y}_{lm0}, \qquad \boldsymbol{y}_{lm2} = \boldsymbol{\hat {r}}Y_{lm}. \end{equation}
These vector spherical harmonics follow the properties of orthogonality by dot product and inner product.
Appendix B. Algebraic details for Stokes flow solution
B.1. Lamb’s solution
Using vector spherical harmonics in (A2), one can define the velocity basis sets
$\boldsymbol{u^{\pm }}_{lmq}$
(
$q = 0, 1, 2$
) where
$\pm$
represents growing or decaying harmonics, respectively,
\begin{align}\boldsymbol{u^{-}}_{lm0}& = \frac {1}{2}r^{-l}(2-l+lr^{-2})\boldsymbol{y}_{lm0} + \frac {1}{2}r^{-l}[l(l+1)]^{1/2}(1-r^{-2})\boldsymbol{y}_{lm2},\end{align}
\begin{align}\boldsymbol{u^{-}}_{lm1} &= r^{-l-1}\boldsymbol{y}_{lm1},\end{align}
\begin{align}\boldsymbol{u^{-}}_{lm2} & = \frac {1}{2}r^{-l}(2-l)\left (\frac {l}{l+1}\right )^{1/2}(1-r^{-2})\boldsymbol{y}_{lm0} + \frac {1}{2}r^{-l}(l+(2-l)r^{-2})\boldsymbol{y}_{lm2},\end{align}
\begin{align}\boldsymbol{u^{+}}_{lm0}& = \frac {1}{2}r^{l-1}(-(l+1)+(l+3)r^{2})\boldsymbol{y}_{lm0} - \frac {1}{2}r^{l-1}[l(l+1)]^{1/2}(1-r^{2})\boldsymbol{y}_{lm2},\end{align}
\begin{align}\boldsymbol{u^{+}}_{lm1} &= r^{l}\boldsymbol{y}_{lm1},\end{align}
\begin{align}\boldsymbol{u^{+}}_{lm2} &= \frac {1}{2}r^{l-1}(3+l)\left (\frac {l+1}{l}\right )^{1/2}(1-r^{2})\boldsymbol{y}_{lm0} + \frac {1}{2}r^{l-1}((l+3)-(l+1)r^{2})\boldsymbol{y}_{lm2}.\end{align}
In terms of these basis sets, the velocity fields inside and outside the vesicle are
\begin{align}\boldsymbol{u}^{\textit{out}}& = \sum _{lmq}c^{\infty }_{lmq} \boldsymbol{u}^{+}_{lmq}(\boldsymbol{x}) + \sum _{lmq}( c_{lmq}^{-} - c^{\infty }_{lmq} ) \boldsymbol{u}^{-}_{lmq}(\boldsymbol{x}),\end{align}
\begin{align}\boldsymbol{u}^{\textit{in}} &= \sum _{lmq}c_{lmq}^{+}\boldsymbol{u}^{+}_{lmq}(\boldsymbol{x}).\end{align}
In the above equation,
$c_{lmq}^{\infty }$
are the coefficients associated with the far field
$\boldsymbol{u}^{\infty } = y \hat {\boldsymbol{x}}$
, which correspond to
$c_{2 \,\pm\, 2 0}^{\infty } = \mp i \sqrt {\pi /5}$
,
$c_{2\, \pm \,2 2}^{\infty } = \mp i \sqrt {2\pi /15}$
$c_{101}^{\infty } = i \sqrt {2\pi /3}$
. The coefficients
$c_{lmq}^{-}$
and
$c_{lmq}^{+}$
are associated with the disturbance fields outside and inside vesicle, which one must solve.
From the velocity fields in (B3), one can compute the hydrodynamic stresses on the unit sphere
$r = 1$
. These are written as follows:
\begin{equation} \boldsymbol{n} \boldsymbol{\cdot }\boldsymbol{T}^{\textit{in}} = (\lambda - 1) \sum _{lmq} \tau _{lmq}^{\textit{in}} \boldsymbol{y}_{lmq}, \qquad \boldsymbol{n} \boldsymbol{\cdot }\boldsymbol{T}^{\textit{out}} = \sum _{lmq} \tau _{lmq}^{\textit{out}} \boldsymbol{y}_{lmq}. \end{equation}
Using the notation and results from Bławzdziewicz et al. (Reference Bławzdziewicz, Vlahovska and Loewenberg2000), the coefficients
$\tau _{lmq}$
are related to the velocity coefficients
$(c_{lmq}^{\infty }, c_{lmq}^{+}, c_{lmq}^{-})$
as follows (details for
$\varTheta _{q^{\prime}q}$
are given in Vlahovska & Gracia Reference Vlahovska and Gracia2007):
\begin{equation} \tau ^{\textit{out}}_{lmq} = \sum _{q^{\prime}=0}^{2}c^{\infty }_{lmq^{\prime}}\big(\varTheta ^{+}_{q^{\prime}q}-\varTheta ^{-}_{q^{\prime}q}\big) + \sum _{q^{\prime}=0}^{2}c_{lmq^{\prime}}^{-}\varTheta ^{-}_{q^{\prime}q} ,\qquad \tau ^{\textit{in}}_{lmq} =\sum _{q^{\prime}=0}^{2}c_{lmq^{\prime}}^{+}\varTheta ^{+}_{q^{\prime}q}. \end{equation}
B.2. Solving equations
Solving the Stokes equations reduces to solving the unknown coefficients
$(c_{lmq}^+, c_{lmq}^-)$
associated with the disturbance velocity fields, as well as the surface tension
$\sigma _{lm}$
on the interface. For each mode
$(l, m)$
, there are seven unknowns to solve:
\begin{equation} \left(c_{lm0}^{+}, c_{lm1}^{+}, c_{lm2}^{+}, c_{lm0}^{-}, c_{lm1}^{-}, c_{lm2}^{-}, \sigma _{lm}\right)\!. \end{equation}
These are found by applying continuity of velocity across the surface
$\boldsymbol{u}^{\textit{out}} = \boldsymbol{u}^{\textit{in}}$
at
$r = 1$
(three equations), membrane incompressibility
$\boldsymbol{\nabla}_{\kern-1pt s} \boldsymbol{\cdot }\boldsymbol{u}^{\textit{in}} = 0$
at
$r = 1$
(one equation), and traction balance
$\chi \boldsymbol{f}^{hydro} = \boldsymbol{f}^{\textit{mem}}$
at
$r = 1$
(three equations), where
$\boldsymbol{f}^{\textit{mem}}$
is Taylor expanded to
$O(\varDelta ^{1/2})$
on the unit sphere.
Applying the boundary conditions is straightforward except the traction balance. Below, we write the membrane traction,
\begin{equation} \boldsymbol{f}^{\textit{mem}} = {t}^{\textit{mem}}_{lm0}\boldsymbol{y}_{lm0}+{t}^{\textit{mem}}_{lm2}\boldsymbol{y}_{lm2} \end{equation}
where
\begin{align}{t}^{\textit{mem}}_{lm0} &= \left (\left [-2\beta - \frac {\mu _{0}}{\alpha Cn^{2}}\right ]q_{lm} - \sigma _{lm}\right )\sqrt {l(l+1)},\end{align}
\begin{align}{t}^{\textit{mem}}_{lm2} &= 2\beta l(l+1) q_{lm} \nonumber \\ &\quad + \left (\frac {(l+2)(l-1)g_{0}}{\alpha Cn^{2}} + (1+\beta q_{0})(l-1)(l)(l+1)(l+2)+ \sigma _{0}(l-1)(l+2)\right )f_{lm} \nonumber \\ &\quad +\, 2\sigma _{lm}.\end{align}
In the above equations, the quantities
\begin{equation} g_{0} = \frac {1}{2}{a}q_{0}^{2} + \frac {1}{4}{b}q_{0}^{4}; \qquad \mu _{0} = aq_{0} + bq_{0}^{3} \end{equation}
represent the double well potential and chemical potential evaluated at the base state
$q_0$
. The dimensionless quantities
$\alpha , Cn, \beta$
are defined in table 2, which correspond to the line tension parameter, Cahn number and dimensionless bending stiffness difference.
In the dynamical equation for
$f_{lm}$
in the main text (3.2), one needs an expression for
$c_{lm2}^{+}$
. After considerable algebra, one obtains expression (3.3), with coefficients
$C_{lm}$
,
$A_l$
,
$B_l$
and
$D_l$
given as follows:
\begin{align}C_{lm} &= \frac {\left[c_{lm0}^{\infty }\sqrt {l(l+1)}(2l+1) + c^{\infty }_{lm2}(4l^{3}+6l^{2}-4l-3)\right]}{d(\lambda ,l)},\end{align}
\begin{align}A_{l} &= -\frac {(l-1)l(l+1)(l+2)}{d(\lambda ,l)}\left [ (1+\beta q_{0})l(l+1) + \frac {g_{0}}{\alpha Cn^{2}}\right ],\end{align}
\begin{align}B_{l} &= -\frac {(l-1)l(l+1)(l+2)}{d(\lambda ,l)},\end{align}
\begin{align}D_l &= -\frac {2(l-1)l(l+1)(l+2)\beta }{d(\lambda ,l)} - \frac {2l(l+1)\mu _{0}}{\alpha Cn^{2}d(\lambda ,l)},\end{align}
where
\begin{equation} d(\lambda ,l) = \frac {9}{\lambda } + (2l^{3}+3l^{2}-5). \end{equation}
Appendix C. Numerical scheme to solve ODEs
We start with the following ODEs for the deformation and order parameter dynamics (for more details of the equations, see Venkatesh Reference Venkatesh2025):
\begin{align}\frac {{\textrm{d}} f_{lm}}{{\textrm{d}}t} &= \frac {im}{2} f_{lm} + \lambda ^{-1}C_{lm} + \lambda ^{-1}\chi ^{-1} \left [ ( A_l + B_l \sigma _0) f_{lm} + D_l q_{lm} \right ],\end{align}
\begin{align}\frac {{\textrm{d}}q_{lm}}{{\textrm{d}}t} &= \frac {im}{2}q_{lm} - \frac {1}{\textit{Pe}} \big[ \varLambda _{lm} l(l+1) + Cn^{2}l^{2}(l+1)^{2}q_{lm} \big. \nonumber\\ & \quad \big. +\, 2\alpha Cn^{2}\beta (l-1)l(l+1)(l+2)f_{lm} \big].\end{align}
In the above equations,
$f_{lm},q_{lm}$
and
$\sigma _{0}$
are unknown. We start by splitting
$f_{lm},q_{lm}$
into their real and imaginary parts – e.g.
$f_{lm}=f^{\prime}_{lm}+if^{\prime \prime}_{lm}$
– while keeping in mind that
$f_{l-m} = (-1)^{|m|}f^{*}_{l|m|}$
. For the deformation, this gives us equations of the form
$(m\geq 0)$
,
\begin{align}{\frac {{\textrm{d}}f^{\prime}_{lm}}{{\textrm{d}}t} = g\left(f^{\prime}_{lm},f^{\prime \prime}_{lm},q^{\prime}_{lm},q^{\prime \prime}_{lm},\lambda ,l,\chi \right) + \sigma _{0} \mathscr{M}\left(f^{\prime}_{lm},\lambda ,l,\chi \right)},\end{align}
\begin{align}{\frac {{\textrm{d}}f^{\prime \prime}_{lm}}{{\textrm{d}}t} = h\left(f^{\prime}_{lm},f^{\prime \prime}_{lm},q^{\prime}_{lm},q^{\prime \prime}_{lm},\lambda ,l,\chi \right) + \sigma _{0}\mathscr{P}}\left(f^{\prime}_{lm},\lambda ,l,\chi \right),\end{align}
where
$g,h,\mathscr{M}$
and
$\mathscr{P}$
are functions of shape and deformation.
At time
$t=t_{i}$
, we perform a first-order Euler predictor–corrector scheme to march the vesicle shape (C3a)–(C3b). We first take a predictor step while neglecting the
$\sigma _{0}$
terms
\begin{align}{\frac {\tilde {f}^{\prime }_{lm}-f^{{\prime }_{i}}_{lm}}{\Delta t} = g\left(f^{{\prime }_{i}}_{lm},f^{{\prime {\prime }}_{i}}_{lm},q^{{\prime }_{i}}_{lm},q^{{\prime {\prime }}_{i}}_{lm},\lambda ,l,\chi \right)},\end{align}
\begin{align}{\frac {\tilde {f}^{\prime {\prime }}_{lm}-f^{{{\prime }{\prime }}_{i}}_{lm}}{\Delta t} = h\left(f^{{\prime }_{i}}_{lm},f^{{{\prime }{\prime }}_{i}}_{lm},q^{{\prime }_{i}}_{lm},q^{{{\prime }{\prime }}_{i}}_{lm},\lambda ,l,\chi \right)},\end{align}
where
$\tilde {f}^{\prime}_{lm},\tilde {f}^{\prime \prime}_{lm}$
are predictor values. We now proceed to take a corrector step to time
$t=t_{i+1}$
, now considering only the
$\sigma _{0}$
portion of (C3a)–(C3b),
\begin{align}{\frac {{f^{\prime}_{lm}}^{i+1}-\tilde {f}^{\prime}_{lm}}{\Delta t} = \sigma _{0}\mathscr{M}(\tilde {f^{\prime}}_{lm}, \lambda , l, \chi )},\end{align}
\begin{align}{\frac {{f^{\prime \prime}_{lm}}^{i+1}-\tilde {f}^{\prime \prime}_{lm}}{\Delta t} = \sigma _{0}\mathscr{P}(\tilde {f^{\prime \prime}}_{lm}, \lambda , l, \chi )}.\end{align}
We choose
$\sigma _{0}$
such that the shape at
$t = t_{i\,+\,1}$
satisfies the area constraint
$A = 4\pi + \varDelta$
, in other words
\begin{equation} \sum _{l,m\geq 0}{\frac {(l-1)(l+2)}{2} \Big ( \big ({f^{\prime}_{lm}}^{i\,+\,1} \big )^{2} +\big({f^{\prime \prime}_{lm}}^{i\,+\,1}\big )^{2} \Big )} = \varDelta .\end{equation}
The latter step is done with root finding (using fzero in MATLAB with tolerance
$10^{-8}$
). The reason for performing this step is to ensure that the nonlinear equations do not drift the differential–algebraic equation system away from the constraint (Ascher, Chin & Reich Reference Ascher, Chin and Reich1994). This method seems to perform better than simply imposing
$\sum _{lm}(l-1)(l+2) {{\textrm{d}}f_{lm}}/{{\textrm{d}}t}f^{*}_{lm} = 0$
to calculate
$\sigma _{0}$
analytically.
After solving for shape, we perform time stepping for the Cahn–Hilliard equation in (C2). We perform the procedure from (Yoon et al. Reference Yoon, Jeong, Lee, Kim, Kim, Lee and Kim2020). We treat the equation implicitly except the nonlinear term
$aq+bq^{3}$
, which we split into an explicit term
$(a-2)q + bq^3$
and an implicit term
$2q$
. The procedure is shown below. After converting
$q_{lm}$
into real and imaginary components, we obtain
\begin{align}\frac {{q^{\prime}_{lm}}^{i+1} - {q^{\prime}_{lm}}^{i+1} }{\Delta t} &= \frac {-m}{2}{q^{\prime \prime}_{lm}}^{i+1} \nonumber\\ & \quad -\, \frac {1}{\textit{Pe}} \left[{\textrm{Re}} \left({\varLambda ^{\prime}_{lm}}^{i}\right)l(l+1) + 2{q^{\prime}_{lm}}^{i+1}l(l+1) + Cn^{2}l^{2}(l+1)^{2}{q^{\prime}_{lm}}^{i+1} \right. \nonumber\\ & \quad \left. +\, 2\alpha Cn^{2}\beta (l-1)l(l+1)(l+2){f^{\prime}_{lm}}^{i+1} \right],\end{align}
\begin{align}\frac {{q^{\prime \prime}_{lm}}^{i+1} - {q^{\prime}_{lm}}^{i+1} }{\Delta t} &= \frac {m}{2}{q^{\prime}_{lm}}^{i+1} \nonumber\\ & \quad -\, \frac {1}{\textit{Pe}} \left[{\textrm{Im}} \left({\varLambda ^{\prime}_{lm}}^{i}\right)l(l+1) + 2{q^{\prime \prime}_{lm}}^{i+1}l(l+1) + Cn^{2}l^{2}(l+1)^{2}{q^{\prime \prime}_{lm}}^{i+1} \right. \nonumber\\ & \quad \left. +\, 2\alpha Cn^{2}\beta (l-1)l(l+1)(l+2){f^{\prime \prime}_{lm}}^{i+1}\right],\end{align}
where
$\varLambda ^{\prime}_{lm} = \int _{\varOmega } (({a}-2)q+{b}q^{3} )Y_{lm}d\varOmega$
, and
$\textrm{Re} (.)$
and
$\textrm{Im} (.)$
represent real and imaginary parts of an expression.
Appendix D. Convergence of the spectral method in the absence of flow
For the Cahn–Hilliard equation without flow, we benchmarked our simulation against a sharp-interface limit theory developed by Rätz Reference Rätz2016, where the authors developed an axisymmetric solution to Cahn–Hilliard dynamics on a sphere. The authors solve the following equation:
\begin{equation} \frac {\partial q}{\partial t} = D{\nabla} ^{2}_{s}\big (aq+bq^{3} - \gamma ^{2}\boldsymbol{\nabla} _{\kern-2pt s}^{2}q\big ) \end{equation}
where
$t = \gamma t^*$
, with
$t^*$
being a modified time. They solve this equation for an initial condition being a striped domain where
$q=1$
from
$\theta = 0$
to
$\theta = \theta _1 = 0.43$
,
$q = -1$
from
$\theta = \theta _1 = 0.43$
to
$\theta = \theta _2 = 0.92$
, and
$q = 1$
from
$\theta =\theta _2 = 0.92$
to
$\theta = \pi$
. In the sharp interface limit as
$\gamma \rightarrow 0$
, the three domains eventually merge into two domains, with the steady-state solution being
$\theta _1 = 0$
and
$\theta _2 = 0.812$
for
$a = 1$
,
$b = -1$
and
$D = 1$
.
For different values of
$\gamma$
, we plot
$\theta _1$
versus
$t^*$
, where
$\theta _1$
represents the first point at which the order parameter changes from a positive to a negative value. We then calculate the sum of squared errors between the numerical solution
$\theta _{\gamma }$
and the semianalytical solution
$\theta _s$
(valid at
$\gamma = 0$
) for
$\theta _1$
. The semianalytical solution is given by the following differential equation for
$\theta _1$
:
Sum of squared errors compared with analytical solution of surface Cahn–Hilliard equation Rätz (Reference Rätz2016). The values of
$\gamma$
are
$\gamma =0.15,0.1,0.08,0.035$
. In the simulations,
$D=1$
. The sum is calculated over all points in time. The time is simulated until
$t=0.1$
.
\begin{equation} \frac {{\textrm{d}}\theta _{i}}{{\textrm{d}}t^*} = -\frac {c^{+}}{\sin \theta _{i}} \end{equation}
where
$i=1,2$
and
$c^{+} = -({\sqrt {2}}/{3}) ( ( {\cot {\theta _{1}}+\cot {\theta _{2}}}) / ((\log {|\tan {({\theta _{1}}/{2})}|} )+ (\log |\tan ({\theta _{2}}/{2})| )))$
. Figure 22 shows the error for different values of
$\gamma$
and shows that the numerical solution converges to the analytical solution in the sharp interface limit
$\gamma \rightarrow 0$
. Note that we mainly deal with the diffuse interface limit in our study and do not need the
$\gamma$
values to be extremely small.
Appendix E. Algebra,
$\boldsymbol{Pe} \boldsymbol{\rightarrow \infty}$
theory (
$\boldsymbol{l\,=\,2}$
modes excited in-plane)
We examine (3.2)–(3.3) and use the expressions
$f_{22} = 2 \epsilon _1 \sqrt {{2\pi }/{15}} ( a_2 - i b_2 )$
and
$q_{22} = 2 \epsilon _2 \sqrt {{2\pi }/{15}} \exp (it)$
. If we define a modified capillary number as
$C = \chi /\epsilon _1$
and modified time
$\tau = t/\epsilon _1$
, the system of equations at leading order for
$C \sim O(1), \tau \sim O(1)$
, and
$\epsilon _2/\epsilon _1 \sim O(1)$
are
\begin{eqnarray} \frac {{\textrm{d}} a_2}{{\textrm{d}} \tau } &=& \lambda ^{-1} C^{-1} \left [ (A_2 + B_2 \sigma _0) a_2 + {\frac {\epsilon _2}{\epsilon _1} D_2 \cos (\epsilon _1 \tau )} \right ]\!, \nonumber\\ \frac{{\textrm{d}} b_2}{{\textrm{d}} \tau } &= &\lambda ^{-1} C^{-1} \left [ (A_2 + B_2 \sigma _0) b_2 - {\frac {\epsilon _2}{\epsilon _1} D_2 \sin (\epsilon _1 \tau )} \right ] - \lambda ^{-1} \frac {1}{2i} \sqrt {\frac {15}{2\pi }} C_{22}. \end{eqnarray}
In the above equation,
$A_2, B_2, C_{22}$
and
$D_2$
are coefficients defined in (B10)–(B11) for
$l = 2, m = 2$
. The surface tension
$\sigma _0$
is a Lagrange multiplier to enforce the surface area constraint – i.e.
$a_2^2 + b_2^2 = {15}/{32\pi }$
. We solve for
$\sigma _0$
using the expression
$a_2 ({{\textrm{d}} a_2}/{{\textrm{d}} \tau }) + b_2 ({{\textrm{d}} b_2}/{{\textrm{d}} \tau }) = 0$
, and plug it back into the above ODE (E1) to yield
\begin{eqnarray} \frac {{\textrm{d}} a_2}{{\textrm{d}} \tau } &=& -\varTheta b_2, \nonumber\\ \frac {{\textrm{d}} b_2}{{\textrm{d}} \tau } &=& \varTheta a_2, \end{eqnarray}
where
$\varTheta$
is
\begin{equation} \varTheta = -\frac {1}{a_2^2 + b_2^2} \left [ \frac {1}{2i} \sqrt {\frac {15}{2\pi }} \lambda ^{-1} C_{22} a_2 + {\frac {\epsilon _2}{\epsilon _1} \lambda ^{-1} C^{-1} D_2 \left ( a_2 \sin (\epsilon _1 \tau ) + b_2 \cos (\epsilon _1 \tau ) \right )} \right ]\! .\end{equation}
Substituting the expressions for
$C_{22}$
and
$D_2$
(i.e.
$C_{22} = -60 \lambda i \sqrt {{2\pi}/{15}} ( 9 + 23 \lambda )^{-1}$
,
$D_2 = -48 \lambda \beta ( 9 + 23 \lambda )^{-1}$
for
$q_0 = 0$
), as well as the constraint
$a_2^2 + b_2^2 = {15}/{32\pi }$
yields (4.4).
In Gera et al. Reference Gera, Salac and Spagnolie2022, the authors obtained the same form of the differential equation as (E2) for
$a_2$
and
$b_2$
for a 2-D vesicle, when the radius is written as
$r = 1 + \epsilon a_2(t) \cos (2\phi ) + \epsilon b_2 \sin (2\phi )$
. However, the expression for
$\varTheta$
is different due the fact the vesicle is 2-D rather than 3-D. Their expression is (using the notation in our paper)
\begin{equation} \varTheta _{2D} = \frac {1}{\lambda } \frac {1}{Q^2} \left [ a_2 + {\frac {\epsilon _2}{\epsilon _1}\beta C^{-1} \left ( a_2 \sin (\epsilon _1 \tau ) + b_2 \cos (\epsilon _1 \tau ) \right )} \right ] \end{equation}
where
$Q$
is a constant equal to
$Q^2 = {1}/{3}\sum _{n=1}^{\infty }(n^2 - 1) (a_n^2 + b_n^2)$
. If one defines a 2-D excess length parameter as
$\varDelta _{2D} = L/R - 2\pi$
, where
$L$
is the membrane length and area
$A = \pi R^2$
, one obtains
\begin{equation} \varDelta _{2D} = \frac {3\pi }{2} \epsilon _1^2 Q^2. \end{equation}
Keeping consistent with the notation used in our paper, if we define the small parameter
$\epsilon = \varDelta _{2D}^{1/2}$
, this yields
$Q^2 = {2}/{3\pi }$
, which gives the final expression written in the main text (4.5).
High
$\textit{Pe}$
swinging motion for
$\alpha =0.1,Cn=0.1,\beta =0.5,\textit{Pe}=1000,\chi =2,q_{0}=0, a=-0.1,b=1,\varDelta =0.0001,\lambda =2$
. The dashed lines represent the full numerical solution of (3.2), (3.4) whereas the circles represent
$f^{\prime}_{22},f^{\prime \prime}_{22},q_{22}^{\prime},q^{\prime \prime}_{22}$
for the
$\textit{Pe} \rightarrow \infty$
theory equations (4.2), (4.3). The blue circles/dashed curves represent the
$f^{\prime \prime}_{22},q^{\prime \prime}_{22}$
modes. The black circles/dashed curves represent the
$f^{\prime}_{22},q^{\prime}_{22}$
modes.
High
$\textit{Pe}$
tumbling motion for
$\alpha =0.1,Cn=0.1,\beta =0.5,\textit{Pe}=1000,\chi =0.05,q_{0}=0, a=-0.1,b=1,\varDelta =0.0001,\lambda =2$
. We use
$l_{\textit{max}}=8$
in these simulations. The dashed lines represent the full numerical solution of (3.2), (3.4) whereas the circles represent
$f^{\prime}_{22},f^{\prime \prime}_{22},q_{22}^{\prime},q^{\prime \prime}_{22}$
for the
$\textit{Pe} \rightarrow \infty$
theory equations (4.2), (4.3). The blue circles/dashed curves represent the
$f^{\prime \prime}_{22},q^{\prime \prime}_{22}$
modes. The black circles/dashed curves represent the
$f^{\prime}_{22},q^{\prime}_{22}$
modes.
Appendix F. Comparison between theory and simulation,
$\boldsymbol{Pe} \boldsymbol{\rightarrow \infty}$
, in-plane deformations and concentrations
This section performs full numerical simulations to check the validity of
$\textit{Pe} \rightarrow \infty$
theory discussed in § 4.2. Figures 23 and 24 show cases where swinging and tumbling occur at
$\textit{Pe}=1000$
. The circles indicate the theory (4.3) for the
$q_{lm}$
and
$f_{lm}$
modes, while the dashed lines are the full numerical simulations (solving (3.2) and (3.4)). In both simulations, we start with an initial condition
$f^{\prime}_{22}(0)=\sqrt {\varDelta /10},f^{\prime \prime}_{22}(0)=-\sqrt {\varDelta /10},q^{\prime}_{22}(0)=2q^{+}\sqrt {2\pi /15}$
. We observe that simulations closely match the theory, and this behaviour generally holds for
$\textit{Pe} \sim O(100)$
and above. Another interesting observation is the validity of the expressions for the
$\textit{Pe} \to \infty$
(4.3). These expressions hold only for
$\varDelta \ll 1$
. As
$\varDelta$
increases, the rotational term in (3.2) becomes significant and starts influencing the dynamics as well. Hence, the approximation given by expressions (4.3) are not accurate as
$\varDelta$
increases. Lastly, we generally observe in our simulations that the out-of-plane
$f_{20}$
deformations rapidly decay to zero, and thus this mode does not seem to affect the theory.
Open access
