2.1. Governing equations for vesicle membrane interacting with background fluid

We consider a vesicle immersed in a background fluid within a rectangular domain $\varOmega$ , whose membrane is represented by a closed surface $\varGamma$ . The system is modelled as a two-phase flow, with $\varGamma$ serving as the interface between the inner and outer fluids. To capture complex shape deformations, we adopt a phase-field formulation in which the sharp interface $\varGamma$ is replaced by a diffuse interface characterized by a scalar field $\phi (\boldsymbol{x}, t)$ . The phase-field variable $\phi$ follows a hyperbolic tangent profile across the membrane interface, ensuring a smooth transition from the interior ( $\phi = 1$ ) to the exterior ( $\phi = -1$ ) of the vesicle. Specifically, $\phi$ can be defined as

(2.1) \begin{align} \phi (\boldsymbol{x}, t) = \tanh \left (\frac { r(\boldsymbol{x}, t)}{2\sqrt {2}\varepsilon _\phi }\right )\!, \quad \boldsymbol{x} \in \bar {\varOmega }, \end{align}

where $r(\boldsymbol{x}, t)$ denotes the signed distance from $\boldsymbol{x}$ to the nearest point on the membrane interface $\varGamma$ , taken positive inside and negative outside, and $\varepsilon _\phi$ controls the thickness of the diffuse interface.

Both the inner and outer fluids are assumed to be incompressible, with equal density ( $\rho _{{in}} = \rho _{\textit{out}}$ ) but distinct viscosities, denoted by $\mu _{{in}}$ and $\mu _{\textit{out}}$ . The governing equations for the coupled two-phase flow and membrane evolution are formulated as follows:

(2.2) \begin{align} &(\rho \boldsymbol{u})_t + \boldsymbol{\nabla } \boldsymbol{\cdot }(\rho \boldsymbol{u} \otimes \boldsymbol{u}) = \boldsymbol{\nabla } \boldsymbol{\cdot }(\mu \unicode{x1D63F}(\boldsymbol{u})) - \boldsymbol{\nabla }\! p + \boldsymbol{F} +\boldsymbol{\nabla } \boldsymbol{\cdot }(\zeta (\phi ) \unicode{x1D64B} \lambda _{\textit{local}}), \\[-12pt]\nonumber\end{align}

(2.3) \begin{align} &\boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u}=0, \\[-12pt]\nonumber\end{align}

(2.4) \begin{align} &\phi _t+\boldsymbol{\nabla } \boldsymbol{\cdot }(\phi \boldsymbol{u} )= \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{q_{\phi }}, \\[-12pt]\nonumber\end{align}

(2.5) \begin{align} &\xi \varepsilon _\phi ^2 \boldsymbol{\nabla } \boldsymbol{\cdot }(\phi ^2 \boldsymbol{\nabla } \lambda _{\textit{local}}) + \zeta (\phi ) \unicode{x1D64B} : \boldsymbol{\nabla } \boldsymbol{u} = 0,\end{align}

where (2.2) and (2.3) represent the conservation of momentum and fluid incompressibility via the Navier–Stokes equations, and (2.4) describes the evolution of the vesicle membrane using a phase-field Cahn–Hilliard (CH) formulation. The membrane inextensibility condition is enforced through (2.5). Here, $\boldsymbol{u}$ and $p$ denote the bulk velocity and pressure, respectively. The rate of strain tensor is given by $ \unicode{x1D63F}(\boldsymbol{u}) = (\boldsymbol{\nabla } \boldsymbol{u} + \boldsymbol{\nabla } \boldsymbol{u}^{{T}})/2$ . The spatially varying viscosity is interpolated as $\mu (\phi ) = {(1-{\phi })\mu _{in}}/{2} + {(1+{\phi })\mu _{out}}/{2}$ . The vesicle–fluid interaction force $\boldsymbol{F}$ and the interfacial flux $\boldsymbol{q}_\phi$ will be defined later.

In addition, membrane inextensibility is modelled using the diffuse domain approach of Aland et al. (Reference Aland, Egerer, Lowengrub and Voigt2014), which extends the local surface constraint $\boldsymbol{\nabla }_{\varGamma } \boldsymbol{\cdot }\boldsymbol{u} = \unicode{x1D64B} : \boldsymbol{\nabla } \boldsymbol{u} = 0$ from the vesicle membrane $\varGamma$ to the entire computational domain $\varOmega$ . Here, the surface gradient operator is defined as

(2.6) \begin{align} \boldsymbol{\nabla }_{\varGamma } = \unicode{x1D64B}\, \boldsymbol{\nabla } = \boldsymbol{\nabla } - \boldsymbol{n} (\boldsymbol{n} \boldsymbol{\cdot }\boldsymbol{\nabla }), \end{align}

where $ \unicode{x1D64B} = \boldsymbol{I} - \boldsymbol{n} \otimes \boldsymbol{n}$ denotes the projection tensor onto the tangential plane of the membrane, with the unit normal vector $\boldsymbol{n} = -\boldsymbol{\nabla } \phi / |\boldsymbol{\nabla } \phi |$ pointing from the interior to the exterior of the vesicle. The inextensibility constraint is enforced by introducing a local Lagrange multiplier $\lambda _{local}$ together with a harmonic regularization term, leading to (2.5), where the parameter $\xi \gt 0$ controls the strength of the regularization. The regularized surface delta function is then defined as

(2.7) \begin{align} \zeta (\phi ) = \eta _\phi \left ( \frac {f(\phi )}{\varepsilon _\phi } + \frac {\varepsilon _\phi }{2} |\boldsymbol{\nabla }\phi |^2 \right )\!, \end{align}

where $f(\phi ) = {(1 - \phi ^2)^2}/{4}$ is a double-well potential, $\varepsilon _\phi$ is the membrane thickness and $\eta _\phi$ is a scaling parameter that relates the diffuse interface to the sharp-interface model. The resulting constraint term contributes additional membrane tension forces in the momentum equation, completing the formulation of (2.2).

Furthermore, the velocity $\boldsymbol{u}$ is subject to Dirichlet boundary conditions, while the phase-field variable $\phi$ and its associated flux $\boldsymbol{q}_{\phi }$ satisfy no-flux boundary conditions,

(2.8) \begin{align} &\boldsymbol{u} =\boldsymbol{u}_{\boldsymbol{bc}} \,\,\,\,\textrm {on}\,\,\,\,\partial \varOmega ,\end{align}

(2.9) \begin{align} &\boldsymbol{n_b} \boldsymbol{\cdot }\boldsymbol{\nabla } \phi =\boldsymbol{n_b} \boldsymbol{\cdot }\boldsymbol{q_{\phi }}=0\,\,\,\,\,\,\,\,\textrm {on}\,\,\,\,\partial \varOmega ,\end{align}

where $\boldsymbol{n_b}$ denotes the unit outward normal vector to the boundary $\partial \varOmega$ of the rectangular computational domain $\varOmega$ .

2.2. Governing equation for multicomponent on the vesicle membrane

We model the vesicle membrane $\varGamma$ as a two-phase fluid consisting of liquid-ordered ( $l_o$ ) and liquid-disordered ( $l_d$ ) phases with distinct physical properties. The membrane composition is described by a second scalar phase-field variable $c_{\varGamma }$ defined on the membrane $\varGamma$ , where $c_{\varGamma } = 1$ corresponds to the $l_d$ phase, and $c_{\varGamma } = -1$ to the $l_o$ phase. The two components are separated by a diffuse interface of thickness $\varepsilon _c$ , across which $c_{\varGamma }$ transitions smoothly. The phase separation dynamics on $\varGamma$ are governed by the surface CH equation (Zhiliakov et al. Reference Zhiliakov, Wang, Quaini, Olshanskii and Majd2021),

(2.10) \begin{align} \left (c_{\varGamma }\right )_t + \boldsymbol{{\boldsymbol{\nabla} }_{\varGamma }} \boldsymbol{\cdot }\left ( c_{\varGamma } \boldsymbol{u} \right ) &= \boldsymbol{{\boldsymbol{\nabla} }_{\varGamma }} \boldsymbol{\cdot }\boldsymbol{q_c} , \end{align}

where the diffusive flux of the membrane component, $\boldsymbol{q}_c$ , will be specified later.

To couple the surface phase dynamics (2.10) with the bulk phase-field system (2.2)–(2.5) in the rectangular domain $\varOmega$ , we employ the diffuse domain method (Rätz & Voigt Reference Rätz and Voigt2006; Li et al. Reference Li, Lowengrub, Rätz and Voigt2009; Teigen et al. Reference Teigen, Li, Lowengrub, Wang and Voigt2009), which is an accurate and efficient numerical approach for solving partial differential equations (PDEs) on complex surfaces. The vesicle membrane $\varGamma$ is embedded into the larger, rectangular domain $\varOmega$ using a regularized surface delta function $\zeta (\phi )$ derived from the phase-field variable $\phi$ . This embedding eliminates the need for explicit interface tracking. Moreover, the method establishes a unified, energy-based framework by extending surface PDEs into the bulk domain $\varOmega$ through the regularized delta function $\zeta (\phi )$ (2.7). Accordingly, the surface CH (2.10) is reformulated in $\varOmega$ as

(2.11) \begin{align} \left (\zeta (\phi ) c_{\varGamma }\right )_t + \boldsymbol{\nabla } \boldsymbol{\cdot }\left (\zeta (\phi ) c_{\varGamma } \boldsymbol{u} \right ) &= \boldsymbol{\nabla } \boldsymbol{\cdot }\left ( \zeta (\phi ) \boldsymbol{q_c} \right ). \end{align}

Details of this extension are provided in Appendix A. Notably, previous studies have demonstrated that in the diffuse domain method, the extended PDE (2.11) asymptotically converges to the original surface PDE (2.10) as the diffuse interface thickness approaches zero ( $\varepsilon _{\phi } \to 0$ ) (Li et al. Reference Li, Lowengrub, Rätz and Voigt2009; Yu, Chen & Thornton Reference Yu, Chen and Thornton2012; Lervåg & Lowengrub Reference Lervåg and Lowengrub2015; Chadwick et al. Reference Chadwick, Stewart, Enrique, Du and Thornton2018; Poulsen & Voorhees Reference Poulsen and Voorhees2018; Yu, Guo & Lowengrub Reference Yu, Guo and Lowengrub2020; Guo et al. Reference Guo, Yu, Lin, Wise and Lowengrub2021; Guo Reference Guo2023).

In addition, the membrane component $c_\varGamma$ and its flux $\boldsymbol{q_{c}}$ are subject to no-flux boundary conditions on the boundary of the rectangular domain $\varOmega$ ,

(2.12) \begin{align} \boldsymbol{n_b} \boldsymbol{\cdot }\boldsymbol{\nabla c_{\varGamma }} =\boldsymbol{n_b} \boldsymbol{\cdot }\boldsymbol{q_{c}}=0\,\,\,\,\textrm {on}\,\,\,\,\partial \varOmega . \end{align}

Thanks to the diffuse domain method, the membrane component equation for $c_{\varGamma }$ (2.11) can be naturally coupled with the vesicle–fluid interaction system (2.2)–(2.5) and solved simultaneously within the same rectangular domain $\varOmega$ . The unknown terms $\boldsymbol{F}$ , $\boldsymbol{q_{\phi }}$ and $\boldsymbol{q_{c}}$ in (2.2)–(2.5) and (2.11) will be derived via an energy-variation approach based on the principle of free energy dissipation, thereby ensuring thermodynamic consistency.

2.3. Free energy of the multicomponent vesicle–fluid system

The free energy $\mathcal{F}$ of a multicomponent vesicle immersed in a background fluid comprises contributions from both the fluid and the membrane. These include the kinetic energy of the background fluid $ \mathcal{F}_{\textit{K}}$ , the membrane bending energy $ \mathcal{F}_{\textit{B}}$ , the surface energy $ \mathcal{F}_{\textit{S}}$ , the line energy $ \mathcal{F}_{\textit{L}}$ and a penalty term $ \mathcal{F}_{\textit{A}}$ enforcing surface area conservation. The free energy is expressed as

(2.13) \begin{align} \mathcal{F} &= \mathcal{F}_{{K}} + \mathcal{F}_{{B}} + \mathcal{F}_{{S}} + \mathcal{F}_{{L}} + \mathcal{F}_{{A}} \nonumber\\ &= \int _{\varOmega } \frac {1}{2} \rho |\boldsymbol{u}|^2 \, \mathrm{d}V + \int _{\varOmega } \frac {1}{2 \varepsilon _\phi } \kappa _B(c_\varGamma ) \omega _\phi ^2 \, \mathrm{d}V + \int _{\varOmega } \sigma _S(c_\varGamma ) \zeta (\phi ) \, \mathrm{d}V \nonumber\\ &\quad + \int _{\varOmega } \sigma _L \zeta (c_\varGamma ) \zeta (\phi ) \, \mathrm{d}V \,\,+ \frac {M_s}{2} \frac {(S(\phi ) - S(\phi _0))^2}{S(\phi _0)}. \end{align}

In the bending energy $\mathcal{F}_{{B}}$ , the scalar field

(2.14) \begin{align} \omega _\phi = \eta _\phi \left ( \frac {f'(\phi )}{\varepsilon _\phi } - \varepsilon _\phi \varDelta \phi \right )\!, \end{align}

serves as a phase-field approximation of the membrane mean curvature. Here, zero spontaneous curvature is assumed and the contribution from Gaussian curvature is neglected. The bending rigidity $\kappa _B(c_\varGamma )$ , associated with the bending energy $\mathcal{F}_{{B}}$ , and the surface tension $\sigma _S(c_\varGamma )$ , defined in the surface energy $\mathcal{F}_{{S}}$ , are both modelled as functions of the membrane component field $c_\varGamma$ ,

(2.15) \begin{align} \kappa _B(c_\varGamma ) = \frac {1 + c_\varGamma }{2} \kappa _{B_{l_d}} + \frac {1 - c_\varGamma }{2} \kappa _{B_{l_o}}, \quad \sigma _S(c_\varGamma ) = \frac {1 + c_\varGamma }{2} \sigma _{S_{l_d}} + \frac {1 - c_\varGamma }{2} \sigma _{S_{l_o}}. \end{align}

In the line energy $ \mathcal{F}_{{L}}$ , the function

(2.16) \begin{align} \zeta (c_\varGamma ) = \eta _c \left ( \frac {f(c_\varGamma )}{\varepsilon _c} + \frac {\varepsilon _c}{2} |\boldsymbol{\nabla }c_\varGamma |^2 \right ) \end{align}

serves as a regularized delta function that captures the diffuse interface between membrane components. Here, the double-well potential $ f(c_\varGamma ) = (1 - c_\varGamma ^2)^2/4$ promotes phase separation, $\varepsilon _c$ represents the interfacial thickness between the components, $\eta _c$ is a scaling parameter and $ \sigma _L$ denotes the line tension. In the penalty energy $\mathcal{F}_{{A}}$ , the vesicle surface area is defined as $S(\phi ) = \int _\varOmega \zeta (\phi ) \, \mathrm{d}V$ . The reference area $S(\phi _0)$ is measured at the initial time $t = 0$ and assumed constant during the simulation. Surface area conservation is enforced by this penalty term, with its strength controlled by the penalty coefficient $M_S$ , which has been found to help improve the accuracy of the method (Aland et al. Reference Aland, Egerer, Lowengrub and Voigt2014).

By applying an energy-variation approach, we derive explicit expressions for the terms $\boldsymbol{F}, \boldsymbol{q_\phi }$ and $\boldsymbol{q_c}$ that were previously left unspecified in the governing equations. Detailed derivations are provided in Appendix B.

2.4. Non-dimensionalization

To non-dimensionalize the model, we choose the characteristic length, velocity and time scales as $L$ , $U$ and $L/U$ , respectively. Fluid properties are normalized with respect to the outer fluid (where $\phi = -1$ ), using characteristic viscosity $\mu _{\textit{out}}$ and density $\rho _{\textit{out}}$ , while membrane properties are normalized by those of the $l_o$ phase ( $c_\varGamma = -1$ ). The non-dimensional form of the multicomponent vesicles–fluid model is then given by

(2.17) \begin{align} &(\rho \boldsymbol{u})_t + \boldsymbol{\nabla } \boldsymbol{\cdot }(\rho \boldsymbol{u} \otimes \boldsymbol{u}) = \frac {1}{\textit{Re}}\boldsymbol{\nabla } \boldsymbol{\cdot }(\mu \unicode{x1D63F}(\boldsymbol{u})) -\boldsymbol{\nabla }\! p +\frac {1}{\textit{Re}} \boldsymbol{F}\! +\!\boldsymbol{\nabla } \boldsymbol{\cdot }\left (\zeta (\phi ) \unicode{x1D64B} \lambda _{\textit{local}}\right )\! , \\[-12pt]\nonumber\end{align}

(2.18) \begin{align} &\boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u}=0, \\[-12pt]\nonumber\end{align}

(2.19) \begin{align} &\phi _t+\boldsymbol{\nabla } \boldsymbol{\cdot }(\phi \boldsymbol{u} ) = \frac {1}{\textit{Pe}} \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{q_{\phi }}, \\[-12pt]\nonumber\end{align}

(2.20) \begin{align} &\xi {\varepsilon _\phi }^{2} \boldsymbol{\nabla } \boldsymbol{\cdot }\left (\phi ^{2} \boldsymbol{\nabla } \lambda _{\textit{local}}\right )+\zeta (\phi ) \unicode{x1D64B}: \boldsymbol{\nabla } \boldsymbol{u}=0, \\[-12pt]\nonumber\end{align}

(2.21) \begin{align}&\left (\zeta (\phi ) c_{\varGamma }\right )_t + \boldsymbol{\nabla } \boldsymbol{\cdot }\left (\zeta (\phi ) c_{\varGamma } \boldsymbol{u} \right )= \frac {1}{{\textit{Pe}}_{\varGamma }} \boldsymbol{\nabla } \boldsymbol{\cdot }\left ( \zeta (\phi ) \boldsymbol{q_c}\right )\!,\end{align}

where

(2.22) \begin{align} \boldsymbol{F}&= g \boldsymbol{\nabla } \phi + {\omega _c \boldsymbol{\nabla } (\zeta (\phi ) c_{\varGamma }) } , \end{align}

(2.23) \begin{align} \boldsymbol{q_{\phi }} &= m(\phi ) \boldsymbol{\nabla } {g}, \end{align}

(2.24) \begin{align} \boldsymbol{q_c} &= m(c_{\varGamma }) \boldsymbol{\nabla } {\omega _c}, \end{align}

with

(2.25) \begin{align} g&= b\omega _{\phi } +\frac {1}{\textit{Ca}} \left ({\kappa _B(c_{\varGamma })}{\eta _\phi } \frac {f^{\prime \prime }(\phi )}{{\varepsilon _\phi }^2} \omega _{\phi } - {\eta _\phi } \varDelta ({\kappa _B(c_{\varGamma })} \omega _{\phi } ) \right ) - {\eta _\phi } {\varepsilon _\phi } \boldsymbol{\nabla } b \boldsymbol{\cdot }\boldsymbol{\nabla } \phi , \end{align}

(2.26) \begin{align} b &= -\omega _{c} c_{\varGamma } +{\textit{Cn}}_{\varGamma } \delta (c_{\varGamma }) + \mathcal{M}_{\text{s}} \frac {S(\phi )-S(\phi _0)}{S(\phi _0)} + Cs_{\varGamma } {\sigma _S(c_{\varGamma })} . \\[9pt] \nonumber \end{align}

The chemical potential associated with the membrane component field $c_\varGamma$ is given by

(2.27) \begin{align} \zeta (\phi )\omega _c &=\frac {1}{\textit{Ca}} \frac {{\varepsilon _\phi }^{-1}}{2}{\kappa _B}^{\prime }(c_{\varGamma }) {\omega _{\phi }}^{2} + Cs {\sigma _S}^{\prime }(c_{\varGamma }) {\zeta (\phi )} \nonumber \\ &\quad + {\textit{Cn}}_{\varGamma } \zeta (\phi )\left ({\eta _c} \frac {f^{\prime }(c_{\varGamma })}{\varepsilon _c}-{\eta _c} {\varepsilon _c} \varDelta c_{\varGamma }\right )- {\textit{Cn}}_{\varGamma } {\eta _c} {\varepsilon _c} \boldsymbol{\nabla } \zeta (\phi ) \boldsymbol{\cdot }\boldsymbol{\nabla }{c_{\varGamma }}. \end{align}

Equations (2.17) and (2.18) describe the dynamics of the incompressible background fluid, while the evolution of the vesicle membrane is governed by (2.19). The local inextensibility of the vesicle is enforced via (2.20), and the lateral phase separation of membrane components is governed by (2.21). The boundary conditions are applied as specified in (2.8), (2.9) and (2.12).

Notably, we emphasize that $\sigma (c_\varGamma )$ and $\lambda _{\textit{local}}$ play fundamentally different roles in our model. The field $\lambda _{\textit{local}}$ arises solely from the geometric inextensibility constraint and therefore cannot represent composition-dependent surface tension. In contrast, $\sigma (c_\varGamma )$ is a genuine material property that encodes the membrane’s intrinsic energetics, and its spatial gradients can actively drive Marangoni-type tangential flows. Thus, the two quantities operate at distinct physical levels: $\lambda _{\textit{local}}$ enforces kinematic inextensibility, whereas $\sigma (c_\varGamma )$ governs the membrane’s thermodynamic response to compositional heterogeneity. Although their associated tangential stresses may locally reinforce or partially counteract each other under certain flow conditions, the constraint-induced tension cannot generally replace or compensate for material surface tension since it is uniquely determined by the inextensibility requirement. A systematic investigation of how constraint-induced tension redistributes in the presence of heterogeneous material tension is beyond the scope of the present work and will be considered in future studies.

The vesicle dynamics described by (2.17)–(2.27) are governed by a set of dimensionless parameters that characterize the coupling between the background fluid and membrane dynamics, the lateral phase separation and the influence of line tension and component-dependent physical properties. These dimensionless parameters are summarized in table 1.

Table 1.

Physical and dimensionless parameters used in the simulations.

Our model integrates key physical features, including membrane component-dependent bending rigidity heterogeneity, surface tension heterogeneity, line tension, inextensibility and viscosity contrasts between the internal and external fluids. Its thermodynamically consistent formulation ensures free energy dissipation and physical consistency, enabling the simulation of complex dynamic processes under realistic biological conditions. This comprehensive framework provides a powerful tool for investigating vesicle behaviour under shear flow and for elucidating the role of membrane heterogeneity in governing non-equilibrium dynamics.

3.1. Phase separation on stationary multicomponent vesicles

We first investigate phase separation dynamics on the membrane of multicomponent vesicles in comparison with the experimental results reported in Wang et al. (Reference Wang, Palzhanov, Quaini, Olshanskii and Majd2022). That study examined the phase behaviour of GUVs and provided detailed observations of domain evolution. To facilitate a direct comparison, we simplify our set-up by assuming a stationary vesicle with a fixed spherical shape, focusing on capturing the essential features of membrane phase separation.

Here $\bar {\varOmega } = [-1.5, 1.5]^3$ and is discretized using a uniform Cartesian grid of $128^3$ . The fixed spherical vesicle is represented by a phase-field with radius $r_0=1$ , corresponding to a physical vesicle radius $R = 10 \,{\unicode{x03BC}}$ m (Wang et al. Reference Wang, Palzhanov, Quaini, Olshanskii and Majd2022), and membrane thickness $\varepsilon _\phi = 0.02$ . The initial phase-field variable $\phi (\boldsymbol{x}, 0)$ is prescribed as

(3.1) \begin{align} \phi (\boldsymbol{x}, 0) = \tanh \left (\frac {r_0 - r(\boldsymbol{x})}{2\sqrt {2}\varepsilon _\phi }\right )\!, \quad r(\boldsymbol{x}) = \sqrt {x^2 + y^2 + z^2}, \quad \boldsymbol{x} =(x,y,z) \in \bar {\varOmega }, \end{align}

representing a smooth diffuse interface centred around the vesicle surface. To model an initially homogeneous liposome, we first generate random numbers $ c_{\textit{rand}}(\boldsymbol{x})$ with a uniform distribution between 0 and 1. Then, we apply the following cutoff for the initial surface fraction $ c_\varGamma (\boldsymbol{x}, 0)$ :

(3.2) \begin{align} c_{\varGamma }(\boldsymbol{x}, 0) = \begin{cases} \,\,1 \quad (l_{d}), & \text{if } c_{\textit{rand}}(\boldsymbol{x}) \in [0, a_{l_d}], \\ -1 \quad (l_{o}), & \text{if } c_{\textit{rand}}(\boldsymbol{x}) \in (a_{l_d}, 1], \end{cases} \end{align}

where $a_{l_d}$ and $a_{l_o} = 1 - a_{l_d}$ represent the prescribed area fractions of the $l_d$ and $l_o$ phases, respectively.

No-flux boundary conditions are applied to the phase-field variables $\phi$ , $c_{\varGamma }$ and $\omega _c$ at the domain boundary $\partial \varOmega$ .

The model parameters are chosen to match those reported in Wang et al. (Reference Wang, Palzhanov, Quaini, Olshanskii and Majd2022): line tension $\sigma _L = 1.2 \,\text{pN}$ and diffusion coefficient $D_c = 10^{-5}\,\text{cm}^2\,\text{s}^{-1}$ . The lipid mobility is adjusted such that ${\textit{Cn}}_\varGamma /{\textit{Pe}}_\varGamma = 0.12$ . Bending rigidity and surface tension are assumed to be uniform across the $l_o$ and $l_d$ phases, with $\kappa _B = \sigma _S = 1$ . The interfacial thickness of membrane components is set to $\varepsilon _c = 0.02$ , fitted to match experimental GUVs.

Two area fractions, $a_{l_d} = 0.3$ and $a_{l_d} = 0.7$ , are considered to match the experimental conditions in Wang et al. (Reference Wang, Palzhanov, Quaini, Olshanskii and Majd2022). To account for stochastic variability, we perform 10 independent simulations with different random seeds of $c_{\varGamma }$ for each area fraction.

Figures 2 and 3 show the phase separation process on the membrane for area fractions $ a_{l_d} = 0.3$ and $ a_{l_d} = 0.7$ , respectively, compared with the experimental results from Wang et al. (Reference Wang, Palzhanov, Quaini, Olshanskii and Majd2022). During phase separation, membrane components coarsen, forming distinct, stable domains, as illustrated in the rightmost images of both figures 2 and 3 that closely match the experimental snapshots.

Figure 2.

Phase separation on GUVs at $a_{l_d} = 0.3$ : experimental snapshots (a) from Wang et al. (Reference Wang, Palzhanov, Quaini, Olshanskii and Majd2022) and numerical results (b), showing $l_o$ domains (blue) forming within the $l_d$ phase (red).

Figure 3.

Phase separation on GUVs at $a_{l_d} = 0.7$ : experimental snapshots (a) from Wang et al. (Reference Wang, Palzhanov, Quaini, Olshanskii and Majd2022) and numerical results (b), showing $l_d$ domains (red) forming within the $l_o$ phase (blue).

To further validate the accuracy of our model, we compare the time evolution of the number of lipid domains obtained from our numerical simulations with the experimental measurements reported in Wang et al. (Reference Wang, Palzhanov, Quaini, Olshanskii and Majd2022) (digitized using WebPlotDigitizer https://automeris.io/). In the simulations, lipid domains are identified directly on the membrane surface by extracting the isosurface $\phi =0$ , interpolating $c_\varGamma$ onto this surface, thresholding at $c_\varGamma =0$ and applying MATLAB’s standard surface-connected-component labelling routine bwconncomp (https://www.mathworks.com/help/images/ref/bwconncomp.html.) As shown in figure 4, for both area fractions under investigation ( $a_{l_d} = 0.3$ and $a_{l_d} = 0.7$ ), the experimental data (markers) are generally bounded by the minimum (green) and maximum (red) values from the ensemble of our simulations. The average numerical result (blue) also lies in close agreement with the experimental trend.

Figure 4.

Time evolution of number of lipid rafts for $a_{l_d} = 0.3$ (a) and $a_{l_d} = 0.7$ (b). Markers represent experimental data; solid lines represent numerical results (blue, average; green, minimum; red, maximum).

Figure 5.

Time evolution of the total perimeter of lipid domains at $a_{l_d} = 0.3$ (a) and $a_{l_d} = 0.7$ (b). Markers represent experimental data; solid lines represent numerical results.

Figure 6.

Experimental images (a) from Wang et al. (Reference Wang, Palzhanov, Quaini, Olshanskii and Majd2022) and numerical results (b) showing the time evolution of $l_o$ domains (blue) within the $l_d$ phase (red). The simulation is initialized with 12 symmetrically placed $l_o$ domains and an area fraction $a_{l_o} = 0.64$ .

In addition, we further compare the evolution of the total perimeter of lipid domains, quantified by $\int _\varOmega \zeta (\phi ) \zeta (c_\varGamma )\, \mathrm{d}V$ in our simulations, where $\zeta (\phi )$ localizes the membrane and $\zeta (c_\varGamma )$ localizes the phase boundary, with the experimental data (digitized using WebPlotDigitizer), as shown in figure 5. The simulation results (solid lines), averaged over all numerical simulations, exhibit excellent agreement with the experimental measurements (markers) for both $a_{l_d} = 0.3$ and $a_{l_d} = 0.7$ .

As reported in Wang et al. (Reference Wang, Palzhanov, Quaini, Olshanskii and Majd2022), three distinct structural patterns were observed experimentally at an area fraction of $a_{l_d} = 0.3$ prior to equilibrium: (i) the minority $l_d$ phase domains (red) embedded within the majority $l_o$ phase (blue), as shown in figure 2; (ii) majority $l_o$ phase domains (blue) enclosed within the minority $l_d$ phase (red), as illustrated in figure 6; and (iii) a mixed pattern featuring both types of domains, as depicted in figure 7.

Patterns (ii) and (iii) were attributed to the persistence of a metastable third phase, which influences domain coarsening and leads to fluctuations in domain area fractions before equilibrium (Van Santen Reference Van Santen1984; Anderson & Lekkerkerker Reference Anderson and Lekkerkerker2002). For instance, the measured $l_o$ area fraction in case (ii) was approximately 6 % lower than its equilibrium value.

Following Wang et al. (Reference Wang, Palzhanov, Quaini, Olshanskii and Majd2022), we reproduce the experimentally observed patterns by adjusting the initial conditions: (i) prescribing prepositioned, symmetrically arranged $l_o$ domains within an $l_d$ background, and (ii) reducing the $l_o$ area fraction by approximately $6\,\%$ . Accordingly, the second pattern – $l_o$ domains fully enclosed within the $l_d$ phase (figure 6) – is reproduced by initializing 12 symmetrically placed $l_o$ domains with an area fraction $a_{l_o}=0.64$ . The third pattern, featuring coexisting $l_o$ and $l_d$ domains (figure 7), is obtained by initializing six symmetrically distributed $l_o$ domains and setting $a_{l_o}=0.65$ .

The quantitative and qualitative agreement between our reproduced simulations and the experimental observations demonstrates the effectiveness and accuracy of our model in capturing the phase dynamics associated on the vesicle membrane.

Figure 7.

Experimental images (a) from Wang et al. (Reference Wang, Palzhanov, Quaini, Olshanskii and Majd2022) and numerical results (b) with front and rear views over time, showing a combination of $l_o$ phase (blue) and $l_d$ phase (red). The simulation is initialized with six symmetrically placed $l_o$ domains and an area fraction $a_{l_o} = 0.65$ .

3.2. Dynamics of a two-component vesicle under shear flow

To assess the capability of the model in capturing membrane–fluid interactions, we reproduce the experimental observations of multicomponent GUVs under shear flow reported in Tusch et al. (Reference Tusch, Loiseau, Al-Halifa, Khelloufi, Helfer and Viallat2018). These experiments revealed dynamic behaviours such as tank treading, swinging and tumbling, driven by asymmetric phase distributions and periodic elastic energy storage along the contact line.

In our simulations, $\bar {\varOmega } = [-2, 2]^3$ , discretized on a uniform grid with $128^3$ . To match the experimental set-up, the initial vesicle shape is modelled as a prolate ellipsoid with a specified reduced volume. The initial phase variable is given by

(3.3) \begin{align} \phi (\boldsymbol{x}, 0) = \tanh \left (\frac {r_0 - r(\boldsymbol{x})}{2\sqrt {2}\varepsilon _\phi }\right )\!, \quad r (\boldsymbol{x})= \sqrt {\frac {x^2}{d_1^2} + \frac {y^2}{d_2^2} + \frac {z^2}{d_3^2}}, \quad \boldsymbol{x} = (x, y, z) \in \bar {\varOmega }, \end{align}

where $d_1 \gt d_2 = d_3$ , resulting in elongation along the $X$ -axis (see figure 8). The membrane component $c_{\varGamma }$ is initially prepositioned at the tips of the vesicle along the $X$ -axis, with a prescribed area fraction $a_{l_d}$ representing the fraction of the $l_d$ phase.

Figure 8.

(a) Geometric parameters used to characterize vesicle dynamics under shear flow. The vesicle shape is described by its principal semiaxes $(d_1,d_2,d_3)$ , and the inclination angle $\theta$ is defined as the angle between the major axis $d_1$ and the flow direction (X-axis). (b) Initial phase distribution used in the simulations of § 4, shown in the X–Z view.

Figure 9.

Experimental images (a) from Tusch et al. (Reference Tusch, Loiseau, Al-Halifa, Khelloufi, Helfer and Viallat2018) and numerical results (b) showing tank treading and swinging at $ {\textit{Cn}}_\varGamma = 8$ , for a vesicle with reduced volume $v = 0.95$ , flattening $f = 0.3$ and area fraction $a_{l_d} = 0.7$ .

The shear flow is introduced by imposing boundary conditions. Specifically, Dirichlet boundary conditions $\boldsymbol{u}_{bc} = (\gamma z, 0, 0)$ are imposed at the top and bottom boundaries in the $z$ -direction, where $\gamma$ denotes the dimensionless shear rate (see table 1). No-flux boundary conditions are enforced for all phase-field variables ( $\phi$ , $\omega _\phi$ , $c_\varGamma$ , $\omega _c$ ) in the $z$ -direction, while periodic boundary conditions are applied in the $x$ - and $y$ -directions for all variables.

Physical parameters are mainly taken from Tusch et al. (Reference Tusch, Loiseau, Al-Halifa, Khelloufi, Helfer and Viallat2018): mean vesicle radius $R = 10 \,{\unicode{x03BC}} \text{m}$ , fluid viscosity $\mu _{in} = \mu _{out} = 1.2 \times 10^{-3} \,\text{Pa}\boldsymbol{\, }\text{s}$ , bending rigidity $\kappa _{B_{l_o}} = \kappa _{B_{l_d}} = 10^{-19}\,\text{J}$ and shear rate $\dot {\gamma } = 3\,\text{s}^{-1}$ . Fluid densities are taken as $\rho _{in} = \rho _{out} = 1000 \,\text{kg}\,\text{m}^{-3}$ (Luby-Phelps Reference Luby-Phelps1999b ). Surface tensions are fixed at $\sigma _{S_{l_d}} = \sigma _{S_{l_o}} = 10^{-6}\,\text{N}\,\text{m}^{-1}$ to ensure uniform membrane properties. The line tension $\sigma _L$ is varied as $4,1$ and $0.1\,\text{pN}$ to reproduce the three dynamical regimes observed experimentally.

Using $R$ as the characteristic length and $\dot {\gamma }R$ as the velocity scale, the resulting dimensionless parameters are Reynolds number ${\textit{Re}} = 3 \times 10^{-4}$ , bending capillary number ${\textit{Ca}} = 36$ , Péclet numbers ${\textit{Pe}} = 1000$ and ${\textit{Pe}}_\varGamma = 10$ , dimensionless surface tension strength $Cs_\varGamma = 27$ and surface area penalty coefficient $\mathcal{M}_{{S}} = 2000$ . The Cahn number is set as ${\textit{Cn}}_\varGamma = 8,\,2,\,0.2$ for the three regimes (figures 9–11). Interfacial thickness parameters $\varepsilon _c = 0.05$ for membrane components and $\varepsilon _\phi = 0.05$ for the membrane are fitted to experimental GUVs.

Figures 9–11 illustrate vesicle dynamics under different reduced volumes, showing close agreement with the experimental results from Tusch et al. (Reference Tusch, Loiseau, Al-Halifa, Khelloufi, Helfer and Viallat2018). In figure 9, the initial ellipsoid with $ d_1 = 1$ , $ d_2 = d_3 = 0.7$ , yields $ v \approx 0.95$ , flattening parameter $ f = ({d_1 - d_2})/{d_1} = 0.3$ and area fraction $ a_{l_d} = 0.7$ . At ${\textit{Cn}}_\varGamma =8$ , the vesicle exhibits tank treading and swinging characterized by a steady prolate shape with oscillatory inclination angle.

At a moderate reduced volume $ v \approx 0.8$ , shown in figure 10, with $ d_1=1, d_2 = d_3 = 0.35$ , $ f = 0.65$ and $ a_{l_d} = 0.6$ , the vesicle exhibits pronounced deformation at ${\textit{Cn}}_\varGamma =2$ . It appears elliptical when the $ l_o$ – $ l_d$ interface aligns with the short axis and resembles two adhered vesicles when aligned with the long axis, exhibiting tank treading and swinging accompanied by large deformations.

At lower volume $v \approx 0.7$ (figure 11), with $d_1=1, d_2 = d_3 = 0.25$ , $f = 0.75$ and $a_{l_d} = 0.5$ , the vesicle transitions to a tumbling regime at ${\textit{Cn}}_\varGamma =0.2$ . The interface length between $l_o$ and $l_d$ is markedly reduced along the long axis due to line tension, which drives the system to minimize line energy, facilitating large deformations and promoting tumbling under shear.

Figure 10.

Experimental images (a) from Tusch et al. (Reference Tusch, Loiseau, Al-Halifa, Khelloufi, Helfer and Viallat2018) and numerical results (b) showing strong deformation during tank treading and swinging at $ {\textit{Cn}}_\varGamma = 2$ , for a vesicle with reduced volume $ v = 0.8$ , flattening parameter $ f = 0.65$ and area fraction $ a_{l_d} = 0.6$ .

Figure 11.

Experimental images (a) from Tusch et al. (Reference Tusch, Loiseau, Al-Halifa, Khelloufi, Helfer and Viallat2018) and numerical results (b) showing the tumbling at $ {\textit{Cn}}_\varGamma = 0.2$ , for a vesicle with reduced volume $ v = 0.7$ , flattening parameter $ f = 0.75$ and area fraction $ a_{l_d} = 0.5$ .

To further validate our model, we track the time evolution of the inclination angle $\theta$ of a vesicle undergoing tank treading and swinging behaviours at reduced volume $v = 0.95$ and area fraction $a_{l_d} = 0.7$ (see figure 9). The inclination angle $\theta$ is defined as the angle between the vesicle’s long axis and the flow direction, constrained within $\theta \in [-{\pi }/{2}, {\pi }/{2} ]$ (see figure 8). As shown in figure 12, the simulated results (solid line) closely match the experimental measurements (circles) reported in Tusch et al. (Reference Tusch, Loiseau, Al-Halifa, Khelloufi, Helfer and Viallat2018) (digitized using WebPlotDigitizer), capturing the characteristic oscillations of the tank treading and swinging regime with a mean angle of approximately $30^\circ$ . The numerical results remain consistently within the experimental variability, thereby confirming not only the accuracy and robustness of the model but also its predictive capability in reproducing flow-induced vesicle dynamics.

Figure 12.

Time evolution of inclination angle $ \theta$ for a vesicle with reduced volume $ v = 0.95$ , flattening parameter $ f = 0.3$ and area fraction $ a_{l_d} = 0.7$ at ${\textit{Cn}}_\varGamma =8$ , comparing numerical results (solid line) with experimental data (circles) from Tusch et al. (Reference Tusch, Loiseau, Al-Halifa, Khelloufi, Helfer and Viallat2018).

4.1. Bending heterogeneity-driven dynamics

We systematically explore the coupling between shear flow and bending heterogeneity in multicomponent vesicles by varying the bending rigidity contrast $\kappa _B = \kappa _{B_{l_d}} / \kappa _{B_{l_o}}$ and the bending capillary number ${\textit{Ca}} = \mu _{\textit{out}} U L^2 / \kappa _{B_{l_o}}$ , both readily tuneable in experiments.

The numerical set-up closely follows § 3.2. The vesicle is initialized as a prolate ellipsoid with principal axes $d_1 = 1$ and $d_2 = d_3 = 0.565$ , as defined in (3.3), corresponding to a reduced volume $v \approx 0.91$ . Two symmetric $l_d$ -phase domains (red) are prepositioned at the vesicle tips along the $X$ -axis (see figure 8), with an initial area fraction $a_{l_d} = 0.4$ for the membrane component field $c_{\varGamma }$ . The dimensionless shear rate is set to $\gamma = 1$ .

Moreover, the dimensionless parameters are specified as follows: Reynolds number ${\textit{Re}} = 10^{-3}$ ; Cahn number ${\textit{Cn}}_{\varGamma } = 0.05$ ; Péclet numbers ${\textit{Pe}} = 1000$ and ${\textit{Pe}}_\varGamma = 1$ ; dimensionless surface-tension strength $Cs_{\varGamma } = 10$ ; surface-area penalty coefficient $\mathcal{M}_{\mathrm{S}} = 2000$ ; surface tension contrast $\sigma _S = 1$ . In all cases, the $l_d$ phase is assumed to be softer than the $l_o$ phase, i.e. $\kappa _{B_{l_d}} \lt \kappa _{B_{l_o}}$ .

Figure 13 presents a phase diagram spanned by $(\kappa _B, Ca)$ , capturing diverse dynamics of multicomponent vesicles and identifying six distinct dynamical regimes.

1. (i) Lateral ring banding. The vesicle adopts an oblate, tilted shape, with $ l_d$ phase domains stretching along the edges to form a lateral ring band composed of a single slender domain.

2. (ii) Lateral broken-ring banding. Similar to (I), but the $ l_d$ phase domains form a broken lateral ring composed of two slender domains.

3. (iii) Tank treading. The vesicle assumes a prolate, tilted shape, with $ l_d$ phase domains stably located at the high-curvature tips (reported in Liu et al. (Reference Liu, Marple, Allard, Li, Veerapaneni and Lowengrub2017) and Venkatesh & Narsimhan (Reference Venkatesh and Narsimhan2024)).

4. (iv) Vertical broken-ring banding. The $l_d$ phase domains stretch vertically along the flow direction on the membrane of a prolate, tilted vesicle. This results in a vertical broken ring composed of two slender domains.

5. (v) Vertical ring banding. Similar to (IV), but the $ l_d$ phase domains merge into a complete vertical ring band composed of a single slender domain (reported in Gera & Salac (Reference Gera and Salac2018b )).

6. (vi) Phase treading. The $ l_d$ phase domains migrate periodically along the membrane of a prolate, tilted vesicle from one tip to the other tip in the flow direction, while maintaining the initial domain symmetry (reported in Liu et al. (Reference Liu, Marple, Allard, Li, Veerapaneni and Lowengrub2017), Gera & Salac (Reference Gera and Salac2018b ) and Venkatesh & Narsimhan (Reference Venkatesh and Narsimhan2024)).

Note that, regimes (I), (II) and (IV) have not been previously reported and represent novel dynamical regimes driven by membrane bending heterogeneity under shear flow.

Figure 13.

Phase diagram of vesicle dynamics under shear flow as a function of bending rigidity contrast $\kappa _B$ and capillary number ${\textit{Ca}}$ . Regimes (I)–(VI) are distinguished by colours and markers.

The phase diagram also reveals two prominent ‘dynamical belts’: (i) a ‘low-rigidity full-dynamics belt’ ( $ 0.25 \leq \kappa _B \leq 0.65$ ), where the membrane rigidity contrast is small enough for all six regimes to appear across a wide range of ${\textit{Ca}}$ ; and (ii) a ‘low-capillary full-dynamics belt’ ( $0.027 \lesssim Ca \lesssim 0.043$ ), where relatively weak shear allows all six regimes to emerge over a range of $\kappa _B$ .

Within the ‘low-rigidity full-dynamics belt’ ( $0.25 \leq \kappa _B \leq 0.65$ ), increasing ${\textit{Ca}}$ drives a flow-mediated transition pathway through the dynamical regimes: at low ${\textit{Ca}}$ , weak shear and dominant bending energy lead the vesicle to adopt an oblate shape, with low-rigidity $l_d$ domains accumulating along high-curvature lateral edges, forming (I) lateral ring banding. As ${\textit{Ca}}$ increases moderately, enhanced shear stresses partially rupture the $l_d$ domains at the lateral edges, producing (II) lateral broken-ring banding. With further increase in ${\textit{Ca}}$ , the vesicle elongates into a prolate shape, and the $l_d$ domains stabilize at the high-curvature tips, resulting in (III) tank treading. At higher ${\textit{Ca}}$ , the $l_d$ domains stretch vertically along the membrane, but insufficient membrane fluidity prevents full closure into a ring, resulting in (IV) vertical broken-ring banding. At sufficiently high ${\textit{Ca}}$ , the $l_d$ domains are further stretched and coalesce into (V) vertical ring banding. Under even stronger ${\textit{Ca}}$ , persistent and periodic domain migration emerges, characterized by symmetric phase migration, known as (VI) phase treading.

Outside the ‘low-rigidity full-dynamics belt’ ( $\kappa _B \gt 0.65$ ), the diversity of accessible dynamical regimes progressively narrows. For $0.65 \leq \kappa _B \lt 0.75$ , only regimes (III)–(VI) remain, while for $\kappa _B \gt 0.75$ , regime (III) also vanishes. This reduction in dynamical regimes results from the diminishing bending-rigidity contrast between the $l_d$ and $l_o$ phases, which weakens the tendency of low rigidity $l_d$ domains to localize at high-curvature regions and suppresses bending-driven behaviours such as (I) lateral ring banding.

Notably, within the ‘low-capillary full-dynamics belt’ ( $0.027 \leq Ca \lesssim 0.0435$ ) shown in figure 13, we observe a complete sequence of transitions from regime (I) to regime (VI) as $\kappa _B$ increases. This highlights the strong sensitivity of both domain and shape evolution to bending rigidity contrast under moderate shear, and identifies this belt as a critical parameter space for capturing diverse vesicle dynamical regimes.

We further observe that the critical capillary numbers separating successive regimes-for example, the transition from (IV) to (V) – generally decrease with increasing $\kappa _B$ , except for the (V) to (VI) transition at $\kappa _B \leq 0.45$ . This trend indicates that stronger flow is typically required to induce domain migration and evolution under pronounced bending heterogeneity, reflecting a competition between shear-induced forces and membrane bending rigidity contrast.

Moreover, within the (V) vertical ring banding regime in figure 13, we identify three subregimes based on ring dynamics.

1. (V _a ) Temporary ring. The $l_d$ domains briefly merge into a vertical ring that ruptures and stabilizes into a broken-ring configuration.

2. (V _b ) Cyclic broken and closing ring. The $l_d$ domains repeatedly merge into a vertical ring that undergoes persistent cycles of rupture and reconnection.

3. (V _c ) Stable ring. The $l_d$ domains experiences several rupture-reconnect cycles but eventually stabilizes and remains intact over long times.

These subregimes emerge from transitions driven by the bending heterogeneity $\kappa _B$ and capillary number ${\textit{Ca}}$ . At $\kappa _B \geq 0.65$ , a weak bending heterogeneity reduces the preference of $l_d$ domains for high-curvature tips, restricting dynamics to (V $_{{a}}$ ). For $0.45 \lt \kappa _B \lt 0.65$ , increasing ${\textit{Ca}}$ triggers sequential transitions from (V $_{{a}}$ ) to (V $_{{b}}$ ) and ultimately to (V ${c}$ ). At $\kappa _B \leq 0.45$ , a strong heterogeneity destabilizes the ring, suppressing (V $_{{c}}$ ). Notably, (V $_{{a}}$ ) is distinct from regime (IV): in (V $_{{a}}$ ), a transient ring forms and then ruptures, whereas in (IV), the $l_d$ domains never coalesce into a ring.

Figure 14 illustrates six distinct vesicle dynamic regimes (I)–(VI) at $\kappa _B = 0.65$ within the ‘low-rigidity full-dynamics belt’. Figure 14(a) shows their temporal evolution (see supplementary Movie-1, Movie-2, Movie-3, Movie-4, Movie-5 and Movie-6, corresponding to regimes (I)–(VI), respectively), are available at https://doi.org/10.1017/jfm.2026.11276, while figure 14(b) presents the corresponding bending and line energy profiles from $t = 2$ to $t = 38$ . To capture the key dynamical processes, we mark the critical turning points ( $t_{i}$ for $i=1,\ldots ,9$ , see the caption of figure 14) on the energy curves. These turning points of bending and line energies largely coincide.

Figure 14.

(a) Time evolution of six vesicle dynamic regimes (I)–(VI) at $ \kappa _B = 0.65$ with varying ${\textit{Ca}}$ . (b) Corresponding evolution of bending and line energies. Key time points: $t_1 = 2.6$ , $t_2 = 9$ , $t_3 = 11.5$ , $t_4 = 15$ , $t_5 = 17.5$ , $t_6 = 18.7$ , $t_7 = 25$ , $t_8 = 26.7$ , $t_9 = 32$ .

Figure 15.

(a) Time evolution of subregimes (V $_a$ )–(V $_c$ ) within regime (V) at $ \kappa _B = 0.45$ with varying ${\textit{Ca}}$ . (b) Corresponding evolution of bending and line energies. Key time points: $t_1 = 8.2$ , $t_2 = 11.6$ , $t_3 = 14.3$ , $t_4 = 20.5$ , $t_5 = 25.7$ , $t_6 = 28.9$ , $t_7 = 43.6$ , $t_8 = 62$ , $t_9 = 66$ .

Comparing the $ l_d$ phase domain positions in figure 14(a) with the energy evolution in figure 14(b) provides valuable insights into vesicle dynamics. In regimes (I)–(III), $l_d$ domains (with lower rigidity $\kappa _{B_{l_d}}$ ) localize at high-curvature regions, minimizing bending energy and resulting in similar energy levels. However, due to the relatively low ${\textit{Ca}}$ , bending forces dominate over viscous forces, preventing further reduction in line energy. In regime (IV), the $l_d$ domains stretch from the high-curvature tips towards low-curvature regions away from the tips, leading to a temporary increase in bending energy from $t_1$ to $t_3$ . Unable to overcome the energy barrier required to form a vertical ring, both the bending and line energies stabilize thereafter at $t_7$ . In regime (V), the $l_d$ domains initially stretch towards low-curvature membrane regions away from the tips, increasing bending and line energies from $t_1$ to $t_2$ . At $t_4$ , the domains merge into a vertical ring band, minimizing both energies. After $t_9$ , the domains stabilize into a broken-ring configuration. Strictly, this case at $\kappa _B = 0.65$ corresponds to subregime (V $_{{a}}$ ), with detailed subregime dynamics presented in figure 15. In (VI), $l_d$ domains periodically migrate along the membrane. Bending energy peaks when domains align along flatter regions (e.g. $t_2$ , $t_7$ ) and decreases when they reach high-curvature tips (e.g. $t_4$ , $t_9$ ).

To further investigate the (V) vertical ring banding regime, figure 15 presents three subregimes: (V $_{{a}}$ ), (V $_{{b}}$ ) and (V $_{{c}}$ ). Figure 15(a) shows the temporal evolution of each dynamics subregime, while figure 15(b) shows the corresponding bending and line energies from $t=2$ to $t = 70$ , with key transition points ( $t_i$ , $i = 1, \ldots , 9$ , see the caption of figure 15) marked. Consistent with figure 14, the $l_{d}$ phase domain positions correlate with energy evolution. We find that the three subregimes exhibit similar energy evolution before $t_4$ , as all initially form a vertical ring. Subsequent divergence arises from differences in $l_d$ domain evolution. In (V $_{{a}}$ ), the ring ruptures and stabilizes as a broken-ring after $t_7$ , with both bending and line energies reaching stability. In (V $_{{b}}$ ), the ring undergoes periodic rupture–reconnection, leading to cyclic oscillations of the energy: bending energy increases during rupture (e.g. $t_3$ to $t_4$ ) as more low-rigidity $l_d$ domains migrate to low-curvature regions away from the tips, and decreases during reconnection (e.g. $t_5$ to $t_6$ ). In (V $_{{c}}$ ), after two rupture–reconnection cycles (from $t_3$ to $t_8$ ), the $l_d$ domains form a stable vertical ring at $t_8$ , with both bending and line energies stabilizing thereafter. Notably, the equilibrium energies in (V $_{{c}}$ ) are higher than in (V $_{{a}}$ ) due to (i) a larger fraction of low-rigidity $l_{d}$ domains remaining in low-curvature regions, increasing bending energy and (ii) a longer $l_d$ - $l_o$ interface, elevating line energy.

4.2. Surface tension heterogeneity-driven swinging and tumbling

We next investigate the role of membrane surface tension heterogeneity $\sigma _S = \sigma _{S_{l_d}} / \sigma _{S_{l_o}}$ in modulating vesicle dynamics under shear flow. Within the ‘low-rigidity full-dynamics belt’ at a fixed bending rigidity contrast $ \kappa _B = 0.65$ , we systematically vary $ \sigma _S \leq 1$ and the bending capillary number ${\textit{Ca}}$ .

The numerical set-up follows § 4.1. In all cases, the $l_d$ phase is assumed to possess lower surface tension than the $l_o$ phase, i.e. $\sigma _{S_{l_d}} \lt \sigma _{S_{l_o}}$ .

The phase diagram (figure 16) reveals eight dynamical regimes: the six previously identified (I–VI) and two additional regimes.

1. (vii) Swinging. Small-amplitude oscillations of the vesicle inclination angle, coupled with cyclic rupture and reconnection of the $l_d$ -formed vertical ring, similar to (V $_{{b}}$ ).

2. (viii) Tumbling. Rigid-body-like vesicle rotation, with $l_d$ domains stably anchored at high-curvature tips.

Swinging and tumbling dynamics have been previously reported in single-component vesicles with viscosity contrast, and in multicomponent vesicles with small-deformation bending heterogeneity (Gera et al. Reference Gera, Salac and Spagnolie2022; Venkatesh & Narsimhan Reference Venkatesh and Narsimhan2024) or with shape asymmetry arising from initial phase distributions and bending heterogeneity (Liu et al. Reference Liu, Marple, Allard, Li, Veerapaneni and Lowengrub2017). Here, we identify a novel mechanism: surface tension heterogeneity, without viscosity contrast or asymmetric phase distributions, can also induce swinging and tumbling. As shown in figure 16, regimes (VII) and (VIII) occur only under heterogeneous surface tension ( $\sigma _S \lt 1$ ).

Figure 16.

(a) Phase diagram of vesicle dynamics under shear flow as a function of surface tension contrast $\sigma _S$ and capillary number ${\textit{Ca}}$ at fixed bending rigidity contrast $\kappa _B = 0.65$ . Eight dynamic regimes (I)–(VIII) are identified and distinguished by colours and markers. (b) Inclination angle evolution in regime (VII) at $\sigma _S = 0.85$ , ${\textit{Ca}} = 0.02$ . (c) Inclination angle evolution in regime (VIII) at ${\textit{Ca}} = 1$ for varying $\sigma _S$ .

Under uniform surface tension ( $\sigma _S = 1$ ), moderate capillary numbers ( $0.0189 \lesssim {\textit{Ca}} \lesssim 0.037$ ) drive the vesicle into regime (IV), characterized by vertical broken-ring banding. As $\sigma _S$ decreases to 0.85, enhanced membrane fluidity enables $l_d$ domains to overcome the energy barrier, forming a complete ring. Periodic rupture and reconnection of this vertical ring dynamically alter the inclination angle, resulting in regime (VIII). At high capillary numbers ( $ 0.037 \lesssim Ca \leq 1$ ), uniform surface tension ( $\sigma _S = 1$ ) induces regimes (V), characterized by vertical ring banding, and (VI), defined by phase treading. As $\sigma _S$ decreases to 0.85, regime (VIII) tumbling emerges. This transition occurs because the lower-tension $l_d$ domains become strongly anchored at the high-curvature tips, suppressing their migration and introducing hydrodynamic asymmetry that promotes rigid-body-like rotation.

Overall, decreasing $\sigma _S$ lowers the energetic cost of deforming $l_d$ domains, enhancing their mobility and reducing the critical capillary number ${\textit{Ca}}$ required for successive dynamic regime transitions. As a result, at $0.85 \leq \sigma \lt 1$ only regimes (II)–(IV), (VII) and (VIII) persist, regime (I) becomes unsustainable under shear, while swinging and tumbling emerge as the dominant dynamic responses. In addition, for $0.35 \leq \sigma _S \lt 0.85$ , only regimes (II), (III) and (VIII) persist. In this range, strong surface tension contrast firmly anchors $l_d$ domains at the tips, promoting vesicle rotation and lowering the threshold for tumbling, while suppressing other regimes as domain evolution rapidly stabilizes. When $\sigma _S \lt 0.35$ , only (VIII) is observed across all tested ${\textit{Ca}}$ , indicating that extreme surface tension asymmetry robustly favours rigid-body-like rotation.

The distinct signatures of (VII) swinging and (VIII) tumbling are further quantified in figure 16(b), which show the time evolution of the vesicle inclination angle for regimes (VII) and (VIII). In (VII) Swinging, the angle exhibits step-like segments linked to the rupture and reconnection of the $l_d$ ring. As $l_d$ domains stretch towards low-curvature regions, surface tension gradients reduce the inclination angle; after rupture, domains shift back to the tips, increasing the angle. These non-smooth variations reflect the interplay between domain evolution and shear force. In (VIII) tumbling, further decreasing $\sigma _S$ at a fixed ${\textit{Ca}}$ accelerates the angular rotation, driven by stronger surface tension gradients.

In addition, figure 17 provides a comprehensive depiction of the (VII) swinging and (VIII) tumbling regimes. Figure 17(a) illustrates their temporal evolution (see supplementary Movie-7 and Movie-8, corresponding to regimes (VII)–(VIII), respectively), while figure 17(b) presents the corresponding bending, line and surface energies starting from $t = 2$ to $t = 40$ . The critical time points ( $t_i$ , for $i = 1, \ldots , 9$ ; see caption of figure 17) are marked on the energy curves. The evolution of bending and line energies, which reflects the migration of $l_d$ domains, is consistent with the discussion in § 4.1 and is not repeated here. The following analysis focuses on the variations in surface energy.

In (VII) swinging, as shown in figure 17(a), periodic rupture and reconnection of the ring cause surface energy to oscillate. When the ring reconnects, $l_d$ domains stretch into low-curvature regions, increasing both their area and the surface energy (e.g. around $t_4$ ). When the ring ruptures, $l_d$ domains cluster near the high-curvature tips, reducing their area and lowering surface energy (e.g. around $t_5$ ).

In (VIII) tumbling, the $l_d$ domains remain anchored at the tips throughout the vesicle’s rotation. The vesicle alternates between elongated shapes, tilted along the flow direction (with inclination angle $\theta \in (0, \pi /2)$ ), and flattened shapes, tilted against the flow ( $\theta \in (-\pi /2, 0)$ ). When elongated (e.g. at $t_3$ ), the vesicle adopts a thin, prolate shape, where the $l_o$ domains are located in relatively low-curvature regions. When flattened (e.g. at $t_5$ ), the vesicle deforms into a more oblate shape under shear flow, causing the curvature at the $l_o$ locations to decrease further. This curvature reduction lowers the surface energy. Such cyclic shape evolution, driven by the interplay between shear forces and surface tension heterogeneity, accounts for the periodic surface energy variations observed in the tumbling regime.

Figure 17.

(a) Time evolution of vesicle dynamic regimes (VII) and (VIII) at $ \sigma _S = 0.85$ and $ \kappa _B = 0.65$ . (b) Corresponding evolution of bending, line and surface energies. Key time points: $t_1 = 5.4$ , $t_2 = 10$ , $t_3 = 15.5$ , $t_4 = 17.5$ , $t_5 = 23.6$ , $t_6 = 27.4$ , $t_7 = 30$ , $t_8 = 37$ , $t_9 = 39.5$ .

We note that the phase diagram presented in this section is constructed using a symmetric initial distribution of the $l_d$ phase (see figure 8 b). Although the detailed dynamical trajectories may vary when different initial conditions are used, the qualitative dynamical regimes and their underlying physical mechanisms remain robust. Additional simulations starting from random initial conditions, provided in Appendix E, confirm that the vesicle consistently evolves towards dynamical states that are identifiable as the same regime classes reported here. These results demonstrate that the phase behaviour captured by our model is governed primarily by the interplay between membrane heterogeneity and shear flow, rather than strictly by the particular choice of initial phase distribution.