Derivation and simulation of a two-phase fluid deformable surface model

Abstract

To explore the impact of surface viscosity on coexisting fluid domains in biomembranes we consider two-phase fluid deformable surfaces as model systems for biomembranes. Such surfaces are modelled by incompressible surface Navier–Stokes–Cahn–Hilliard-like equations with bending forces. We derive this model using the Lagrange–d’Alembert principle considering various dissipation mechanisms. The highly nonlinear model is solved numerically to explore the tight interplay between surface evolution, surface phase composition, surface curvature and surface hydrodynamics. It is demonstrated that hydrodynamics can enhance bulging and furrow formation, which both can further develop to pinch-offs. The numerical approach builds on a Taylor–Hood element for the surface Navier–Stokes part, a semi-implicit approach for the Cahn–Hilliard part, higher-order surface parametrizations, appropriate approximations of the geometric quantities, and mesh redistribution. We demonstrate convergence properties that are known to be optimal for simplified subproblems.

1. Introduction

Coexisting fluid domains in biomembranes and in their model systems is an intensively studied field of research. With the possibility to visualize coexisting lipid phases in model membranes by high-resolution fluorescence imaging (Baumgart, Hess & Webb 2003) a strong correlation between domain composition and local membrane curvature can be established. These findings have supported earlier membrane models (Jülicher & Lipowsky 1996; Jiang, Lookman & Saxena 2000; Kumar, Gompper & Lipowsky 2001) and initiated a wealth of theoretical, numerical and experimental studies to understand the relation between composition and curvature (Baumgart et al. 2003, 2005; Veatch & Keller 2003; Veatch et al. 2007; Wang & Du 2008; Lowengrub, Rätz & Voigt 2009; Elliott & Stinner 2010b; Elson et al. 2010; Garcke et al. 2016; Gera & Salac 2018; Zimmermann et al. 2019). This relation has strong biological implications (McMahon & Gallop 2005). We refer to Veatch & Keller (2003), Deserno (2015) and Lipowsky & Dimova (2021) for reviews on the subject. While most studies address multicomponent giant unilaminar vesicles and consider phase separation and coarsening on spherical shapes, more recently multicomponent scaffolded lipid vesicles have been considered. In these systems, non-spherical shapes are stabilized and the effect of spatially varying curvature on phase separation and coarsening has been considered (Fonda et al. 2018, 2019; Rinaldin et al. 2020). These results demonstrate the strong influence of curvature on the spatial arrangement of the lipid phases. On the other hand also the shape evolution is considered. Coexisting fluid domains can lead to bulging and budding events (Baumgart et al. 2003; Lowengrub et al. 2009; Elliott & Stinner 2010b), indicating also the strong influence of composition on the evolving shape. Thus, the correlation between domain composition and membrane curvature acts in both directions. These effects can essentially be modelled by the Jülicher–Lipowsky model (Jülicher & Lipowsky 1996) or appropriate phase field approximations of it (Lowengrub et al. 2009; Haußer et al. 2013). These models essentially extend the classical Helfrich model to multiple components. Instead of a simple $L^2$-gradient flow of the Helfrich energy, with possible constraints, the model allows for dissipation through the simultaneous evolution of the shape and the lipid phases on the surface. Such models can be embedded in bulk flows and serve as interfacial conditions (Sohn et al. 2010; Zhang & Wolgemuth 2022). However, all these studies neglect the effect of surface viscosity. Surface viscosity has been shown to be a key property of biomembranes and their model systems controlling remodelling and with it, the coarsening process, see (Faizi, Dimova & Vlahovska 2022) for discussions and measurements. Already for flat membranes, it has been shown that experimental results on the coarsening rate of these fluid domains can only be quantitatively reproduced by simulations if the fluid properties of the membrane are taken into account (Fan, Han & Haataja 2010; Camlay & Brown 2011). Considering these effects on curved surfaces involves several modelling and numerical subtleties, due to the increased coupling between local curvature and surface fluid velocity. Although there are models in the literature dealing with the coupling between surface flow on curved surfaces and the surrounding bulk flow (Woodhouse & Goldstein 2012; Barrett, Garcke & Nürnberg 2015; Reuther & Voigt 2016), they do not consider coexisting fluid domains. Furthermore, it has been shown that the effect of the bulk fluid can be neglected if the Saffman–Delbrück number (Saffman & Delbrück 1975), which defines a hydrodynamic length relating the viscosity of the membrane and the surrounding bulk fluid, is large. In this case the hydrodynamics is effectively two-dimensional (2-D) on spatial scales smaller than the Saffman–Delbrück number (Henle & Levine 2010). We follow this simplification and consider only surface two-phase flow problems. On stationary surfaces these problems are addressed in Nitschke, Voigt & Wensch (2012), Ambrus et al. (2019), Olshanskii, Palzhanov & Quaini (2022) and Bachini, Krause & Voigt (2023b) and essentially reveal an enhanced coarsening rate due to hydrodynamic effects, similar to the situation in flat space (Fan et al. 2010; Camlay & Brown 2011). To consider evolving surfaces under the influence of surface viscosity requires models for so-called fluid deformable surfaces. Such models have been established as model systems for biomembranes (Arroyo & DeSimone 2009; Torres-Sánchez, Millán & Arroyo 2019). They exhibit a solid–fluid duality, while they store elastic energy when stretched or bent, as solid shells, they flow as viscous 2-D fluids under in-plane shear. This duality has several consequences: it establishes a tight interplay between tangential flow and surface deformation. In the presence of curvature, any shape change is accompanied by a tangential flow and, vice versa, the surface deforms due to tangential flow. Models describing this interplay between curvature and surface flow have been introduced and numerically solved in Torres-Sánchez et al. (2019), Reuther, Nitschke & Voigt (2020) and Krause & Voigt (2023). These numerical approaches provide the basis to explore the impact of surface viscosity on the dynamics of biomembranes. We here extend these models to two-phase flows by combining fluid deformable surfaces with a phase field approximation of the Jülicher–Lipowsky model. The model is systematically derived by a Lagrange–d’Alembert principle, accounting for dissipation by domain coarsening, shape evolution and surface viscosity. We relate the resulting model to known simplified models in the literature. We provide a detailed description of the numerical approach, which uses surface finite elements (SFEM) (Dziuk & Elliott 2013; Nestler, Nitschke & Voigt 2019) and builds on previous developments for one-component fluid deformable surfaces (Krause & Voigt 2023) and surface Navier–Stokes–Cahn–Hilliard-like models on stationary surfaces (Bachini et al. 2023b). The algorithm is used to demonstrate the strong interplay between composition, curvature and hydrodynamics and its implications for bulging and budding processes. Measuring the essential material properties, surface viscosity and bending rigidity for biomembranes is a challenging task, see (Faizi et al. 2022). Reported parameters vary by orders of magnitude. We therefore are not considering a specific system but the general properties by exploring the parameter space. Briefly, surface hydrodynamics enhances the onset of shape deformations and possible resulting topological changes, which has strong biological implications, e.g. in the case of endocytosis and exocytosis (Kaksonen & Roux 2018; Al-Izzi, Sens & Turner 2020).

The paper is structured as follows. In § 2 we introduce the used notation necessary to formulate the surface model, briefly describe the model derivation, formulate the two-phase fluid deformable surface model and relate the equations to known simplified models. Details are considered in the appendices. In § 3 we discuss the considered numerical approach. In § 4 we provide all used parameters, consider convergence studies and explore the parameter space and the implications of the coupling between composition, shape and surface flow. In § 5 we draw conclusions.

2. Continuous model

We derive the full model for two-phase fluid deformable surfaces by applying the Lagrange–d’Alembert principle. We first introduce the necessary notation, then motivate the use of the Lagrange–d’Alembert principle, introduce all ingredients and derive the model. Finally, we relate the model to known simplifications in the literature.

2.1. Notation

We consider a time dependent smooth and oriented surface $\mathcal {S} = \mathcal {S}(t)$ without boundary. Related to $\mathcal {S}$, we denote the outward pointing surface normal $\boldsymbol {\nu }$, the surface projection ${\boldsymbol {P}}=\boldsymbol {I}-\boldsymbol {\nu }\otimes \boldsymbol {\nu }$, the shape operator $\mathcal {B}= -\operatorname {\boldsymbol {\nabla }}_{\!{\boldsymbol {P}}}\!\boldsymbol {\nu }$ and the mean curvature ${\mathcal {H}}= \operatorname {tr}\mathcal {B}$. Note that, under these definitions the unit sphere has negative mean curvature ${\mathcal {H}}=-2$. Let $\phi$ be a continuously differentiable scalar field, ${\boldsymbol {\boldsymbol {u}}}$ a continuously differentiable ${\mathbb {R}}^3$-vector field, and $\boldsymbol {\sigma }$ a continuously differentiable ${\mathbb {R}}^{3\times 3}$-tensor field defined on $\mathcal {S}$. We define the surface tangential gradient $\operatorname {\boldsymbol {\nabla }}_{\!{\boldsymbol {P}}}$ as in Jankuhn, Olshanskii & Reusken (2018) and the componentwise surface gradient $\operatorname {\boldsymbol {\nabla }}_{\!C}$ as in Nitschke, Sadik & Voigt (2022), namely

(2.1) \begin{equation} \left.\begin{aligned} \operatorname{\boldsymbol{\nabla}}_{\!{\boldsymbol{P}}}\!\phi &= {\boldsymbol{P}}\boldsymbol{\nabla}\phi^e, \\ \operatorname{\boldsymbol{\nabla}}_{\!{\boldsymbol{P}}}\!{\boldsymbol{\boldsymbol{u}}} &= {\boldsymbol{P}}\boldsymbol{\nabla}{\boldsymbol{\boldsymbol{u}}}^e{\boldsymbol{P}}, \\ \operatorname{\boldsymbol{\nabla}}_{\!C}\!\boldsymbol{\sigma} &= \boldsymbol{\nabla}\boldsymbol{\sigma}^e {\boldsymbol{P}}, \end{aligned}\right\} \end{equation}

where $\phi ^e$, ${\boldsymbol {\boldsymbol {u}}}^e$ and $\boldsymbol {\sigma }^e$ are arbitrary smooth extensions of $\phi$, ${\boldsymbol {\boldsymbol {u}}}$ and $\boldsymbol {\sigma }$ in the normal direction and $\boldsymbol {\nabla }$ is the gradient of the embedding space ${\mathbb {R}}^3$. The fields $\operatorname {\boldsymbol {\nabla }}_{\!{\boldsymbol {P}}}\!\phi, \operatorname {\boldsymbol {\nabla }}_{\!{\boldsymbol {P}}}\!{\boldsymbol {\boldsymbol {u}}}$ are purely tangential vector and tensor fields, respectively. We define the corresponding divergence operators for a vector field ${\boldsymbol {\boldsymbol {u}}}$ and a tensor field $\boldsymbol {\sigma }$ by

(2.2) \begin{equation} \left.\begin{aligned} \operatorname{div}_{\!{\boldsymbol{P}}}\!{\boldsymbol{\boldsymbol{u}}} &= \operatorname{tr}(\operatorname{\boldsymbol{\nabla}}_{\!{\boldsymbol{P}}}\!{\boldsymbol{\boldsymbol{u}}}), \\ \operatorname{div}_{\!C}(\boldsymbol{\sigma}{\boldsymbol{P}}) &= \operatorname{tr}\operatorname{\boldsymbol{\nabla}}_{\!C}(\boldsymbol{\sigma}{\boldsymbol{P}}), \end{aligned}\right\} \end{equation}

where $\operatorname {tr}$ is the trace operator. The divergence of the tensor $\boldsymbol {\sigma }$, $\operatorname {div}_{\!C}\!\boldsymbol {\sigma }$, leads to a non-tangential vector field even if $\boldsymbol {\sigma }$ is a tangential tensor field. Let $\operatorname {\boldsymbol {\nabla }}_{\!{\mathcal {S}}}$ be the gradient with respect to the covariant derivative on $\mathcal {S}$, as used in Reuther et al. (2020). This operator is defined for scalar fields and tangential vector fields and relates to the tangential operators by $\operatorname {\boldsymbol {\nabla }}_{\!{\boldsymbol {P}}}\!\phi =\operatorname {\boldsymbol {\nabla }}_{\!{\mathcal {S}}}\!\phi$ and ${\operatorname {div}_{\!{\boldsymbol {P}}}\!{\boldsymbol {\boldsymbol {u}}} = \operatorname {div}_{\! {\mathcal {S}}}({\boldsymbol {P}}{\boldsymbol {\boldsymbol {u}}})-({\boldsymbol {\boldsymbol {u}}}\boldsymbol {\cdot }\boldsymbol {\nu }){\mathcal {H}}}$, respectively. We further clarify the relation between the above operators in Appendix A.

The surface $\mathcal {S}$ is given by a parametrization $\boldsymbol {X}$. The material on the surface is described by a material parametrization $\boldsymbol {X}_\mathfrak {m}$, as in Nitschke & Voigt (2022). Both parametrizations relate to each other by $\partial _t\boldsymbol {X}_\mathfrak {m}\boldsymbol {\cdot }\boldsymbol {\nu } =\partial _t\boldsymbol {X} \boldsymbol {\cdot } \boldsymbol {\nu }$, which leads to a Lagrangian perspective in normal direction. In tangential direction the surface and the material can move independently. We define the material velocity by ${\boldsymbol {\boldsymbol {u}}}\unicode{x2254} \partial _t\boldsymbol {X}_\mathfrak {m}$ and the relative material velocity by ${\boldsymbol {\boldsymbol {w}}} := {\boldsymbol {\boldsymbol {u}}} - \partial _t\boldsymbol {X}$. The relative material velocity is a pure tangential vector field. Additional information and the relation to the time derivative used can be found in Appendix B.6.

2.2. Model derivation by the Lagrange–d’Alembert principle

The Lagrange–d’Alembert principle has been explained in Marsden & Ratiu (2013). The approach is a combination of the Lagrange and the Onsager variational principles. The concept of the Lagrange principle has been introduced in Marsden & West (2001), Marsden & Ratiu (2013) and Hairer et al. (2006) and it models the interplay of kinetic and potential energies under total energy conservation. On the other side, the concept of the Onsager variational principle models dissipative systems, where potential energy decreases under a dissipation potential, e.g. $L^2$-gradient flows. Within our context, the Onsager principle is used in Torres-Sánchez et al. (2019) to derive fluid deformable surfaces with the surface Stokes model.

We consider a density function $\phi$ that describes the two-phases on the surface $\mathcal {S}$, e.g. liquid-ordered and disordered phases, where one phase is represented by $\phi =1$, the other one by $\phi =-1$, and we assume a mixture of both if $\phi \in (-1,1)$. We consider a conserved evolution for $\phi$ with respect to the following energies.

The first is a Ginzburg–Landau energy modelling phase-separation on $\mathcal {S}$,

(2.3)\begin{equation} {{\mathcal{F}}_{GL}} = \int_{\mathcal{S}}\tilde{\sigma} \left(\frac{\epsilon}{2}\,{\left\| \operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!\phi \right\|}^2 +\frac{1}{\epsilon}W(\phi) \right) {\operatorname{d}\!{\mathcal{S}}}, \end{equation}

where $\tilde {\sigma }>0$ is a rescaled line tension with $\tilde {\sigma } = 3 / (2 \sqrt {2}) \sigma$, $\sigma$ the line tension in the corresponding Jülicher–Lipowski model (Garcke et al. 2016; Elliott, Hatcher & Stinner 2022; Benes et al. 2023), $\epsilon >0$ is related to the diffuse interface width and $W(\phi ) = \frac {1}{4}(\phi ^2 - 1)^2$ is a double-well potential.

The second energy we consider accounts for bending properties. As in Fonda et al. (2018) and Bachini et al. (2023b), we consider a diffuse interface approximation of the Jülicher–Lipowski model (Jülicher & Lipowsky 1996),

(2.4)\begin{equation} {\mathcal{F}}_{H} = \int_{\mathcal{S}} \frac{1}{2} \kappa(\phi)({\mathcal{H}}-{\mathcal{H}}_0(\phi))^2{\operatorname{d}\!{\mathcal{S}}}, \end{equation}

where $\kappa (\phi )$ is the bending stiffness and ${\mathcal {H}}_0(\phi )$ the spontaneous curvature. We neglect additional contributions due to the Gaussian curvature of $\mathcal {S}$. This is justified as long as the Gaussian bending stiffness is independent on the phase $\phi$ and no topological changes occur. To get a continuously differentiable dependency on $\phi$ we consider an interpolation function as in Elliott & Stinner (2010b) and Bachini et al. (2023b),

(2.5)\begin{equation} f(\phi)=\begin{cases} f_1 & \mbox{if } \phi = 1,\\ \displaystyle \frac{f_1+f_2}{2}+ \frac{f_1-f_2}{4}\phi(3-\phi^2) & \mbox{if } -1< \phi< 1,\\ f_2 & \mbox{if } \phi ={-}1, \end{cases} \end{equation}

for $f \in \{ \kappa, {\mathcal {H}}_0 \}$ and $f_1,f_2$ the material parameter of $\kappa$ or ${\mathcal {H}}_0$ in the separated phases. Together, the energies ${{\mathcal {F}}_{GL}}$ and ${\mathcal {F}}_{H}$ define the potential energy ${\mathcal {F}}= {\mathcal {F}}_{H}+{{\mathcal {F}}_{GL}}$ of the system.

We define the surface material velocity by ${\boldsymbol {\boldsymbol {u}}}$ and the kinetic energy as in Reuther et al. (2020) by

(2.6)\begin{equation} {{\mathcal{F}}_{K}} = \int_{\mathcal{S}}\frac{\rho}{2}{\left\| {\boldsymbol{\boldsymbol{u}}} \right\|}^2{\operatorname{d}\!{\mathcal{S}}}, \end{equation}

with $\rho$ the surface density. For simplicity, we assume $\rho$ to be constant. The Lagrangian $\mathcal {L}={{\mathcal {F}}_{K}}-{\mathcal {F}}$ is defined as the difference between the kinetic and the potential energy.

We next consider the various sources of dissipation. As in Torres-Sánchez et al. (2019) we define the dissipation potential of the viscous stress,

(2.7)\begin{equation} \mathcal{D}_V = \int_{\mathcal{S}} \eta{\left\| \boldsymbol{\sigma} \right\|}^2, \end{equation}

where $\eta$ denotes the viscosity and $\boldsymbol {\sigma }({\boldsymbol {\boldsymbol {u}}}) = \frac {1}{2} (\operatorname {\boldsymbol {\nabla }}_{\!{\boldsymbol {P}}}\! {\boldsymbol {\boldsymbol {u}}} + (\operatorname {\boldsymbol {\nabla }}_{\!{\boldsymbol {P}}}\! {\boldsymbol {\boldsymbol {u}}})^\textrm {T})$ is the rate of deformation tensor as considered in Jankuhn et al. (2018). For simplicity, we assume $\eta$ to be constant. In addition, we consider the friction with the surrounded material, which is modelled by

(2.8)\begin{equation} \mathcal{D}_R = \int_{\mathcal{S}} \frac{\gamma}{2}{\left\| {\boldsymbol{\boldsymbol{u}}} \right\|}^2, \end{equation}

where $\gamma \ge 0$ is a friction coefficient, again assumed to be constant. The third component of the dissipation potential is associated with the dissipation due to phase separation. We assume the immobility potential of the phase field given by

(2.9)\begin{equation} \mathcal{D}_{\phi} = \int_{\mathcal{S}} \frac{1}{2m} \left\Vert \dot{\phi} \right\Vert_{H^{{-}1}}^2, \end{equation}

with $m>0$ a resistance associated with $\dot {\phi }=\partial _t\phi +\boldsymbol {\nabla }_{{\boldsymbol {\boldsymbol {w}}}}\phi$ with respect to $H^{-1}$ norm. Thereby, $\dot {\phi }$ denotes the material time derivative of $\phi$ and $\boldsymbol {\nabla }_{{\boldsymbol {\boldsymbol {w}}}}\phi =(\operatorname {\boldsymbol {\nabla }}_{\!{\mathcal {S}}}\!\phi,{\boldsymbol {\boldsymbol {w}}})$, see Appendix B.6. The parameter $m$ plays the role of a mobility. Here, we follow the approach of Magaletti et al. (2013) and consider this parameter to be constant. For alternative approaches, we refer to Abels (2009). We now define the dissipation potential by $\mathcal {D}= \mathcal {D}_V + \mathcal {D}_R + \mathcal {D}_{\phi }$.

Moreover, we add to $\mathcal {L}$ and $\mathcal {D}$ the following constraints. We here only assume local inextensibility of the material as in Torres-Sánchez et al. (2019) and Reuther et al. (2020). This is incorporated by the Lagrange-function $p$, which serves as the surface pressure and is related to the surface tension. It enters in the constraint by

(2.10)\begin{equation} \mathcal{C}_{IE}={-}\int_{\mathcal{S}}p\operatorname{div}_{\!{\boldsymbol{P}}}\!{\boldsymbol{\boldsymbol{u}}}. \end{equation}

This constraint induces a conservation of surface area $|\mathcal {S}|$. We neglect possible additional constraints, e.g. on the enclosed volume.

With these ingredients, the Lagrange–d’Alembert principle can be applied. The concept is introduced and used for rigid body models in Marsden & West (2001), Udwadia & Kalaba (2002) and Izadi & Sanyal (2014). Here we consider this principle for space-dependent functions, for which it provides an elegant way to derive thermodynamically consistent models. It reads

(2.11)\begin{gather} 0 = (\rho (\partial_t{\boldsymbol{\boldsymbol{u}}}+\boldsymbol{\nabla}_{{\boldsymbol{\boldsymbol{w}}}}{\boldsymbol{\boldsymbol{u}}}) + \boldsymbol{D}_{\boldsymbol{X}}{\mathcal{F}} + \boldsymbol{D}_{{\boldsymbol{\boldsymbol{u}}}}\mathcal{D} + \boldsymbol{D}_{{\boldsymbol{\boldsymbol{u}}}} \mathcal{C}_{IE}, \boldsymbol{Y}) + (\boldsymbol{D}_{\phi}{\mathcal{F}} + \boldsymbol{D}_{\dot{\phi}}\mathcal{D}, \psi) + (\boldsymbol{D}_{p} \mathcal{C}_{IE}, q), \end{gather}

for all test functions $\boldsymbol {Y}\in T{\mathbb {R}}^3\vert _{\mathcal {S}}$ and $\psi,q\in T^0\mathcal {S}$. Here, the symbol $\boldsymbol{D}_{\cdot}$ denotes the $L^2$-gradient of the functional. See Appendix B for the explicit definitions and calculations. The approach is broadly applicable and provides an alternative to other frameworks, such as the Onsager principle. We obtain the following problem.

Problem 2.1 Find $({\boldsymbol {\boldsymbol {u}}},p,\phi,\mu )$ such that

(2.12) \begin{align} \left.\begin{aligned} \partial_t{\phi}+ \boldsymbol{\nabla}_{{\boldsymbol{\boldsymbol{w}}}}\phi &= m{\rm \Delta}_{{\mathcal{S}}}\mu,\\ \mu &=\tilde{\sigma} \left(-\epsilon{\rm \Delta}_{{\mathcal{S}}}\phi+\frac{1}{\epsilon}W^{\prime}(\phi)\right) +\frac{1}{2}\kappa^{\prime}(\phi)({\mathcal{H}}-{\mathcal{H}}_0(\phi))^2 \\ &\quad -\,\kappa(\phi){\mathcal{H}}^{\prime}_0(\phi)({\mathcal{H}}-{\mathcal{H}}_0(\phi)),\\ \rho(\partial_t {\boldsymbol{\boldsymbol{u}}} + \boldsymbol{\nabla}_{{\boldsymbol{\boldsymbol{w}}}}{\boldsymbol{\boldsymbol{u}}}) &={-}\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\! p - p{\mathcal{H}}\boldsymbol{\nu} + 2\eta\operatorname{div}_{\!C}\!\boldsymbol{\sigma} - \gamma {\boldsymbol{\boldsymbol{u}}} + \boldsymbol{b}_{T}+ \boldsymbol{b}_N,\\ \operatorname{div}_{\!{\boldsymbol{P}}}\!{\boldsymbol{\boldsymbol{u}}} &=0, \end{aligned}\right\} \end{align}

with tangential and normal bending forces defined by

(2.13) \begin{align} \left.\begin{aligned} \boldsymbol{b}_T &= \mu \operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!\phi,\\ \boldsymbol{b}_N &={-}\left({\rm \Delta}_{{\mathcal{S}}}(\kappa(\phi)({\mathcal{H}}-{\mathcal{H}}_0(\phi))) +\kappa(\phi)({\mathcal{H}}-{\mathcal{H}}_0(\phi)) \vphantom{\frac{1}{2}}\right.\\ &\left.\quad \times\left(\Vert\mathcal{B}\Vert^2-\frac{1}{2}{\mathcal{H}}({\mathcal{H}}-{\mathcal{H}}_0(\phi))\right)\right)\boldsymbol{\nu} \\ &\quad +~\tilde{\sigma}\left(\frac{\epsilon}{2}\Vert\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!\phi\Vert^2+\frac{1}{\epsilon} W(\phi)\right){\mathcal{H}}\boldsymbol{\nu} -\tilde{\sigma}\epsilon\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!\phi^{\rm T}\mathcal{B}\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!\phi\boldsymbol{\nu}, \end{aligned}\right\} \end{align}

and ${\boldsymbol {\boldsymbol {w}}}$ the relative material velocity as explained above.

This problem provides a tight coupling between the phase field $\phi$, the surface velocity ${\boldsymbol {\boldsymbol {u}}}$, the surface pressure $p$ and the geometry $\mathcal {S}$. Exploring the main implications of these couplings is one of the goals of this paper.

Remark 2.1 (Thermodynamic consistency)

We refer to Appendix B.5 to show that the relation

(2.14)\begin{equation} \frac{{\rm d}}{{\rm d} t} ( {{\mathcal{F}}_{K}} + {\mathcal{F}} ) ={-}2 \mathcal{D} \leq 0, \end{equation}

holds for the model in Problem 2.1. This demonstrates thermodynamic consistency.

Considering a characteristic length and a characteristic velocity, Problem 2.1 can be formulated in non-dimensional form, see Appendix C for details. Keeping the same notation also in non-dimensional form, we obtain the following.

Problem 2.2 Find $({\boldsymbol {\boldsymbol {u}}},p,\phi,\mu )$ such that

(2.15)\begin{equation} \partial_t{\phi}+ \boldsymbol{\nabla}_{{\boldsymbol{\boldsymbol{w}}}}\phi = m {\rm \Delta}_{{\mathcal{S}}}\mu,\end{equation}

(2.16)\begin{align} \mu &=\tilde{\sigma} \left(-\epsilon{\rm \Delta}_{{\mathcal{S}}}\phi+ \frac{1}{\epsilon}W^{\prime}(\phi) \right) +\frac{1}{2}\kappa^{\prime}(\phi)\left({\mathcal{H}}-{\mathcal{H}}_0(\phi)\right)^2 \nonumber\\ &\quad -\kappa(\phi){\mathcal{H}}^{\prime}_0(\phi)\left({\mathcal{H}}-{\mathcal{H}}_0(\phi)\right)\!, \end{align}

(2.17)\begin{gather} \partial_t {\boldsymbol{\boldsymbol{u}}} + \boldsymbol{\nabla}_{{\boldsymbol{\boldsymbol{w}}}}{\boldsymbol{\boldsymbol{u}}} ={-}\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\! p - p{\mathcal{H}}\boldsymbol{\nu} + \frac{2}{Re}\operatorname{div}_{\!C}\!\boldsymbol{\sigma} -\gamma{\boldsymbol{\boldsymbol{u}}}+ \boldsymbol{b}_{T}+ \boldsymbol{b}_N, \end{gather}

(2.18)\begin{gather}\operatorname{div}_{\!{\boldsymbol{P}}}\!{\boldsymbol{\boldsymbol{u}}} = 0, \end{gather}

where $Re$ denotes the Reynolds number and $\boldsymbol {b}_T$ and $\boldsymbol {b}_N$ are defined as in Problem 2.1.

In this form, the problem is suitable for numerical approximation. However, before addressing the model in its discrete form, we consider several model simplifications.

2.3. Model simplifications

The model in Problem 2.2 is a two-component fluid deformable surface model with phase-dependent elasticity. For simplicity, we assume that the material parameters, such as the density $\rho$, the viscosity $\eta$ and the friction coefficient $\gamma$, are constant. It is possible to extend the model and consider phase-dependent parameters by following the approach described, for example, in Lowengrub & Truskinovsky (1998), Abels, Garcke & Grün (2012), Aland & Voigt (2012) and Ten Eikelder et al. (2023). In the following, we link our model to simplified models already discussed in the literature.

2.3.1. One component fluid deformable surface

In the case of a single phase, the simplified model considers only surface hydrodynamics and its interaction with the geometry due to bending. Explicitly, by considering $\phi \equiv \mbox {const}$, the system simplifies to

(2.19) \begin{equation} \left.\begin{aligned} \partial_t {\boldsymbol{\boldsymbol{u}}} + \operatorname{\boldsymbol{\nabla}}_{{\boldsymbol{\boldsymbol{w}}}}{\boldsymbol{\boldsymbol{u}}} &={-}\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\! p - p{\mathcal{H}}\boldsymbol{\nu} + \frac{2}{Re}\operatorname{div}_{\!C}\!\boldsymbol{\sigma} - \gamma{\boldsymbol{\boldsymbol{u}}}+\boldsymbol{b}_N,\\ \operatorname{div}_{\!{\boldsymbol{P}}}\!{\boldsymbol{\boldsymbol{u}}} &= 0, \end{aligned}\right\} \end{equation}

where ${\boldsymbol {\boldsymbol {w}}}={\boldsymbol {\boldsymbol {u}}}-\partial _t\boldsymbol {X}$ as before, and the normal bending forces reduce to

(2.20)\begin{equation} \boldsymbol{b}_N ={-}\kappa\big({\rm \Delta}_{{\mathcal{S}}}{\mathcal{H}} +({\mathcal{H}}-{\mathcal{H}}_0)\big(\Vert\mathcal{B}\Vert^2-\tfrac{1}{2} {\mathcal{H}}({\mathcal{H}}-{\mathcal{H}}_0) \big) \big)\boldsymbol{\nu}. \end{equation}

Note that, bending stiffness and spontaneous curvature are now constant parameters. In the case of ${\mathcal {H}}_0=0$ and $\gamma =0$, this model has been introduced and simulated in Reuther et al. (2020), Krause & Voigt (2023) and in the Stokes limit in Torres-Sánchez et al. (2019). For further numerical and analytical approaches under additional symmetry assumptions we refer to Al-Izzi et al. (2020), where a linear stability analysis of a tube geometry is considered, and to Olshanskii (2023), where potential rotational symmetric equilibrium configurations are addressed. For a comparison of different derivations of this model we refer to Reuther & Voigt (2015, 2018a) and Brandner, Reusken & Schwering (2022b).

2.3.2. Two-component fluid on a stationary surface

If we assume a stationary surface $\mathcal {S}$, i.e. ${\boldsymbol {\boldsymbol {u}}}\boldsymbol {\cdot }\boldsymbol {\nu }=0$, the model in Problem 2.2 restricts to a pure tangential problem. We follow the approach of the directional splitting as shown in detail in Jankuhn et al. (2018) and used in Reuther et al. (2020). We note that ${\boldsymbol {\boldsymbol {u}}}_T={\boldsymbol {P}}{\boldsymbol {\boldsymbol {u}}}$ and $u_N={\boldsymbol {\boldsymbol {u}}}\boldsymbol {\cdot }\boldsymbol {\nu }$ for the tangential and normal surface velocity, respectively. We project the equations in the surface tangent space and this results in a surface two-phase flow problem,

(2.21) \begin{equation} \left.\begin{aligned} \partial_t{\phi}+ \boldsymbol{\nabla}_{{\boldsymbol{\boldsymbol{w}}}}\phi &= m{\rm \Delta}_{{\mathcal{S}}}\mu,\\ \mu&=\tilde{\sigma} \left(-\epsilon{\rm \Delta}_{{\mathcal{S}}}\phi+\frac{1}{\epsilon}W^{\prime}(\phi)\right) +\frac{1}{2}\kappa^{\prime}(\phi)\left({\mathcal{H}}-{\mathcal{H}}_0(\phi)\right)^2 \\ &\quad -\,\kappa(\phi){\mathcal{H}}^{\prime}_0(\phi)\left({\mathcal{H}}-{\mathcal{H}}_0(\phi)\right)\!,\\ {\boldsymbol{P}}\partial_t {\boldsymbol{\boldsymbol{u}}}_T + \boldsymbol{\nabla}_{{\boldsymbol{\boldsymbol{w}}}}{\boldsymbol{\boldsymbol{u}}}_T &={-}\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\! p + \frac{2}{Re}\operatorname{div}_{\! {\mathcal{S}}}\!\boldsymbol{\sigma}({\boldsymbol{\boldsymbol{u}}}_T) -\gamma{\boldsymbol{\boldsymbol{u}}}_T + \mu \operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!\phi,\\ \operatorname{div}_{\! {\mathcal{S}}}\!{\boldsymbol{\boldsymbol{u}}}_T &= 0, \end{aligned}\right\} \end{equation}

where the relative material velocity simplifies to the Eulerian case ${\boldsymbol {\boldsymbol {w}}}={\boldsymbol {\boldsymbol {u}}}_T$. Because of the tangentiality of the vector and tensor fields, all differential operators fall back to the operators with respect to the covariant derivative. In the case of $\gamma =0$, the model has been introduced and discussed in Bachini et al. (2023b), where it is used to study the influence of surface hydrodynamics on coarsening in multicomponent scaffolded lipid vesicles, see Fonda et al. (2018), Fonda et al. (2019) and Rinaldin et al. (2020). Without the bending terms, the system has already been addressed in Nitschke et al. (2012), Ambrus et al. (2019) and Olshanskii et al. (2022).

2.3.3. One component fluid on a stationary surface

By considering a single phase on a stationary surface, the system simplifies even further and leads to the inextensible surface Navier–Stokes equations on a stationary surface,

(2.22) \begin{equation} \left.\begin{aligned} {\boldsymbol{P}}\partial_t {\boldsymbol{\boldsymbol{u}}}_T + \boldsymbol{\nabla}_{{\boldsymbol{\boldsymbol{w}}}}{\boldsymbol{\boldsymbol{u}}}_T &={-}\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\! p + \frac{2}{Re}\operatorname{div}_{\! {\mathcal{S}}}\!\boldsymbol{\sigma}({\boldsymbol{\boldsymbol{u}}}_T) - \gamma{\boldsymbol{\boldsymbol{u}}},\\ \operatorname{div}_{\! {\mathcal{S}}}\!{\boldsymbol{\boldsymbol{u}}} &= 0, \end{aligned}\right\} \end{equation}

where again ${\boldsymbol {\boldsymbol {w}}}={\boldsymbol {\boldsymbol {u}}}_T$. These equations have been solved numerically for simply connected surfaces in Nitschke et al. (2012) and Reuther & Voigt (2015, 2018a) (and in the Stokes limit in Olshanskii et al. (2018), Bonito, Demlow & Licht (2020) and Brandner & Reusken (2020)), and for general surfaces in Nitschke, Reuther & Voigt (2017), Reusken (2018), Fries (2018), Reuther et al. (2020), Reuther & Voigt (2018b) and Lederer, Lehrenfeld & Schoeberl (2020) (and in the Stokes limit in Olshanskii & Yushutin (2019)).

2.3.4. Two-component surface without fluid behaviour

A larger literature exists for the evolution of two-component surfaces if surface hydrodynamics is neglected. These models are related to the overdamped limit of Problem 2.2. An explicit derivation is shown in Appendix D. The resulting model reads

(2.23) \begin{equation} \left.\begin{aligned} \partial_t{\phi}+ \boldsymbol{\nabla}_{{\boldsymbol{\boldsymbol{w}}}}\phi &= \tilde m\,{\rm \Delta}_{{\mathcal{S}}}\tilde\mu, \\ \tilde\mu &={-}\tilde{\tilde{\sigma}}\epsilon{\rm \Delta}_{{\mathcal{S}}}\phi+ \frac{\tilde{\tilde{\sigma}}}{\epsilon}W^{\prime}(\phi), \\ &\quad +\,\frac{1}{2}\tilde\kappa^{\prime}(\phi)\left({\mathcal{H}}-{\mathcal{H}}_0(\phi)\right)^2 -\tilde\kappa(\phi){\mathcal{H}}^{\prime}_0(\phi)\left({\mathcal{H}}-{\mathcal{H}}_0(\phi)\right)\!, \\ u_N &={-}\tilde{p}{\mathcal{H}} +\tilde{b}_N,\\ {\boldsymbol{\boldsymbol{u}}}_T &={-}\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\! \tilde{p} + \tilde\mu\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!\phi,\\ -{\rm \Delta}_{{\mathcal{S}}}\tilde{p} + \tilde{p}{\mathcal{H}}^2 &={-}\operatorname{div}_{\! {\mathcal{S}}}\! \left(\tilde{\mu}\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!\phi\right) + b_N{\mathcal{H}}, \end{aligned}\right\} \end{equation}

where $\tilde {b}_N = \tilde {\boldsymbol {b}}_N\boldsymbol {\cdot }\boldsymbol {\nu }$. In this model, $\tilde {p}$ serves as a Lagrange function to ensure the local inextensibility constraint. The model is similar to a model discussed in Haußer et al. (2013), see Appendix D for a detailed comparison. If the constraint on local inextensibility is dropped and only a global area constraint is considered, the model further simplifies. See Appendix D for the relations with the models considered in Wang & Du (2008) and Elliott & Stinner (2010a), which are derived by only considering variations in the normal direction. One essential difference is the tangential velocity component $\tilde \mu \operatorname {\boldsymbol {\nabla }}_{\!{\mathcal {S}}}\!\phi$, which is not present in these models. We consider the model in Elliott & Stinner (2010a) but with an additional local inextensibility constraint for a qualitative comparison to highlight the differences of the model derivation and the importance of surface hydrodynamics on the evolution of two-component membranes.

3. Numerical discretization

We consider a SFEM (Dziuk & Elliott 2013; Nestler et al. 2019) to solve the highly nonlinear set of geometric and surface partial differential equations (PDEs) in Problem 2.2. The approach uses higher-order surface discretizations, mesh regularization, a Taylor–Hood element for the surface Navier–Stokes equations and established approaches in flat space to split the Navier–Stokes–Cahn–Hilliard-like problem. The discretization on the surface addresses numerical analysis results for vector-valued surface PDEs ensuring convergence of surface vector-Laplace and surface Stokes problems on stationary surfaces (Hansbo, Larson & Larsson 2020; Hardering & Praetorius 2022) and reproduces the same expected optimal order of convergences for the full two-component fluid deformable surface model.

3.1. Mesh movement

By following the work in Krause & Voigt (2023), we combine the system in Problem 2.2 with a mesh redistribution approach introduced in detail in Barrett, Garcke & Nürnberg (2008). We consider the initial surface given $\boldsymbol {X}(0)=\boldsymbol {X}_0$ and an additional equation for the time evolution of the parametrization

(3.1)\begin{gather} \partial_{t}\, \boldsymbol{X} \boldsymbol{\cdot} \boldsymbol{\nu} = {\boldsymbol{\boldsymbol{u}}}\boldsymbol{\cdot}\boldsymbol{\nu}, \end{gather}

(3.2)\begin{gather}{\mathcal{H}} \boldsymbol{\nu} = {\rm \Delta}_C \boldsymbol{X}. \end{gather}

This approach generates a tangential mesh movement that maintains the shape regularity and additionally provides an implicit representation of the mean curvature ${\mathcal {H}}$.

3.2. Surface approximation

We assume that the smooth surface $\mathcal {S}$ is approximated by a discrete $k$th-order approximation $\mathcal{S}_h$. Let $\mathcal {S}_h^{lin}$ be a piecewise linear reference surface given by shape regular triangulation ${\mathcal {T}_{h}}^{lin}=\{T_i\}_{i=1}^{{N_{T}}}$. We define $\boldsymbol {X}$ as bijective map $\boldsymbol {X} \colon \mathcal {S}_h^{lin} \rightarrow \mathcal {S}$ such that $\mathcal {S} = \cup _{i=1}^{{N_{T}}} \boldsymbol {X}(T_i)$. The construction of such maps is discussed in Praetorius & Stenger (2020) and Brandner et al. (2022a). We get a $k$th-order approximation of $\boldsymbol {X}$ by the $k$th-order interpolation $\mathcal {S}_h^k=I_{h}^k(\boldsymbol {X})$, which defines a higher-order triangulation such that $\mathcal {S}_h^k = \cup _{i=1}^{{N_{T}}} \boldsymbol {X}^k_h(T_i)$. We use each geometrical quantity like the normal vector $\boldsymbol {\nu }_h$, the shape operator $\mathcal {B}_h$ and the inner products $({\cdot }, {\cdot })_h$ with respect the $\mathcal {S}_h^k$ where we will drop the index $k$ in the following. We denote the size of the grid by $h$, i.e. the longest edge of the mesh.

3.3. Discrete function spaces and weak formulation

We define the discrete function spaces for scalar function by

(3.3)\begin{equation} V_{k_e}(\mathcal{S}_h)=\{ \psi \in C^0(\mathcal{S}_h) \vert \psi\vert_{T}\in\mathcal{P}_{k_e}(T)\}, \end{equation}

with $\mathcal {P}_{k_e}$ the space polynomials, where we set the element order as $k_e=k$. We define $\boldsymbol {V}_{k}(\mathcal {S}_h)=[V_{k}(\mathcal{S}_h)]^3$ as space of discrete vector fields. We discretize ${\boldsymbol {\boldsymbol {u}}}_h,\boldsymbol {X}_h\in \boldsymbol {V}_2(\mathcal {S}_h)$, ${\mathcal {H}}_h,\phi _h,\mu _h\in V_2(\mathcal {S}_h)$ and $p_h\in V_1(\mathcal {S}_h)$. This corresponds to a Taylor–Hood element for the pair of velocity and pressure. The resulting weak formulation reads

(3.4)\begin{gather} (\partial_t{\phi_h}+ \boldsymbol{\nabla}_{{\boldsymbol{\boldsymbol{w}}}_h}\phi_h,\psi_h)_h ={-}m(\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!\mu, \operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!\psi_h)_h ,\end{gather}

(3.5)\begin{gather} (\mu_h,\xi_h)=\tilde{\sigma} \epsilon(\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!\phi_h,\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!\xi_h)_h -\frac{\tilde{\sigma}}{\epsilon}(W^{\prime}(\phi_h),\xi_h)_h \nonumber\\ \qquad \quad+\,\frac{1}{2}(\kappa^{\prime}(\phi_h)\left({\mathcal{H}}_h-{\mathcal{H}}_0(\phi_h)\right)^2,\xi_h)_h \nonumber\\ \qquad\qquad\qquad\quad -\left(\kappa(\phi_h){\mathcal{H}}^{\prime}_0(\phi_h) \left({\mathcal{H}}_h-{\mathcal{H}}_0(\phi_h)\right), \xi_h\right)_h, \end{gather}

(3.6)\begin{gather} (\partial_t {\boldsymbol{\boldsymbol{u}}}_h + \boldsymbol{\nabla}_{{\boldsymbol{\boldsymbol{w}}}_h}{\boldsymbol{\boldsymbol{u}}}_h,{\boldsymbol{v}_{\scriptscriptstyle{h}}})_h = (p_h,\operatorname{div}_{\!{\boldsymbol{P}}}\!{\boldsymbol{v}_{\scriptscriptstyle{h}}})_h - \frac{2}{Re}(\boldsymbol{\sigma}_h,\operatorname{\boldsymbol{\nabla}}_{\!{\boldsymbol{P}}}\!{\boldsymbol{v}_{\scriptscriptstyle{h}}})_h + (\boldsymbol{b}_{T}+\boldsymbol{b}_N-\gamma{\boldsymbol{\boldsymbol{u}}}_h,{\boldsymbol{v}_{\scriptscriptstyle{h}}})_h, \end{gather}

(3.7)\begin{gather}(\operatorname{div}_{\!{\boldsymbol{P}}}\!{\boldsymbol{\boldsymbol{u}}}_h,q_h)_h = 0, \end{gather}

(3.8)\begin{gather}(\partial_t \boldsymbol{X}_h \boldsymbol{\cdot} \boldsymbol{\nu}_h,h_h)_h = ({\boldsymbol{\boldsymbol{u}}}_h\boldsymbol{\cdot}\boldsymbol{\nu}_h,h_h)_h, \end{gather}

(3.9)\begin{gather}({\mathcal{H}}_h \boldsymbol{\nu}_h, {\boldsymbol{Z}_{\scriptscriptstyle{h}}})_h ={-}(\operatorname{\boldsymbol{\nabla}}_{\!C}\!\boldsymbol{X}_h,\operatorname{\boldsymbol{\nabla}}_{\!C}\!{\boldsymbol{Z}_{\scriptscriptstyle{h}}})_h, \end{gather}

for all $(\psi _h,\xi _h,{\boldsymbol {v}_{\scriptscriptstyle {h}}},q_h,h_h,{\boldsymbol {Z}_{\scriptscriptstyle {h}}})\in [V_k\times V_k \times \boldsymbol {V}_{k} \times V_{k-1}\times V_k \times \boldsymbol {V}_{k}](\mathcal {S}_h)$.

3.4. Discrete model

Following the strategy in Bachini et al. (2023b), we separate the surface phase field equations (3.4)–(3.5) and surface Navier–Stokes equations (3.6)–(3.7) by an operator-splitting approach. In 2-D and three-dimensional (3-D) settings such splitting approaches are common strategies for Navier–Stokes–Cahn–Hilliard systems, and have been shown to be robust and fast converging, see, for example, Demont et al. (2022) for a detailed analysis. Similar properties have been found in Bachini et al. (2023b) on surfaces. In contrast to Bachini et al. (2023b), where the surface was stationary, here the surface evolves and we adapt the numerical schemes used in Krause & Voigt (2023) for one-component systems.

Let $\{ t^n \}_i^N$ be time interval discretization, with $t^n=n\tau$ and $\tau >0$ be the time steps. We first solve the surface phase field problem (3.4)–(3.5), and then solve the surface Navier–Stokes (3.6) and (3.7), together with (3.8)–(3.9) for mesh regularization. The inner product $\left ({\cdot }\,,\,{\cdot }\right )_{h}$ is approximated by a quadrature rule with an order that is chosen high enough such that test and trial functions and area elements are well integrated. We denote the time discrete surface by $\mathcal {S}_h^{n-1}=\mathcal {S}_h(t^{n-1})$.

The equations are discretized in time by a semi-implicit Euler scheme where the nonlinear terms are chosen explicitly apart from the double well potential. As is common for Cahn–Hilliard equations, we linearize the derivative of the double-well potential by a Taylor expansion of order one.

Problem 3.1 (Discrete surface Cahn–Hilliard problem)

Find $(\phi _{h},\mu _{h})\in [{{V}_{h}}\times {{V}_{h}}](\mathcal {S}_h^{n-1})$ such that

(3.10) \begin{equation} \left.\begin{gathered} \frac{1}{\tau}\left({\phi_{h}^n},\,{{\psi_{\scriptscriptstyle{h}}}}\right)_{h} + \left({\operatorname{\boldsymbol{\nabla}}_{{\boldsymbol{\boldsymbol{w}}}^{n-1}}\phi_{h}^n},\,{{\psi_{\scriptscriptstyle{h}}}}\right)_{h} = \frac{1}{\tau}\big({\phi_{h}^{n-1}},\,{{\psi_{\scriptscriptstyle{h}}}}\big)_{h} - m \left({\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!\mu_{h}^{n}},\,{\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!{\psi_{\scriptscriptstyle{h}}}}\right)_{h}\!,\\ \left({\mu_{h}^n},\,{{\xi_{\scriptscriptstyle{h}}}}\right)_{h} = \epsilon \left({\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!\phi_{h}^n},\,{\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!{\xi_{\scriptscriptstyle{h}}}}\right)_{h} + \frac{1}{\epsilon} \big({-2(\phi_{h}^{n-1})^3 + \big(3(\phi_{h}^{n-1})^2-1\big)\phi_{h}^n },\,{{\xi_{\scriptscriptstyle{h}}}}\big)_{h}\\ \hskip-.8pc +\, \frac{1}{2}\big({ \kappa^{\prime}(\phi_{h}^{n-1}) \big({\mathcal{H}}^{n-1}_{h}-{\mathcal{H}}_0(\phi)(\phi_{h}^{n-1})\big)^2},\,{{\xi_{\scriptscriptstyle{h}}}}\big)_{h} \\ \quad\quad\quad\,\,\, - \big({\kappa(\phi_{h}^{n-1}){\mathcal{H}}_0(\phi)^{\prime}(\phi_{h}^{n-1}) \big({\mathcal{H}}^{n-1}_{h}-{\mathcal{H}}_0(\phi)(\phi_{h}^{n-1})\big)},\,{{\xi_{\scriptscriptstyle{h}}}}\big)_{h}\!, \end{gathered}\right\} \end{equation}

for all $({\psi _{\scriptscriptstyle {h}}},{\xi _{\scriptscriptstyle {h}}})\in [{{V}_{h}}\times {{V}_{h}}](\mathcal {S}_h^{n-1})$.

This approach can be related to analysed schemes for the Cahn–Hilliard and Navier–Stokes–Cahn–Hilliard equations with explicit treatment of the double-well potential and additional stabilization terms, for example, Shen & Yang (2010a,b). These methods and more an advanced scalar auxiliary variable (SAV) approaches, for example, Zhu et al. (2019), have been shown to be unconditionally energy-stable in 2-D and 3-D settings.

The surface Navier–Stokes equations and the mesh redistribution are considered together. Moreover, we define a discrete surface update variable $\boldsymbol {Y}_h^{n}=\boldsymbol {X}_h^{n}-\boldsymbol {X}_h^{n-1}$, which is considered as unknown instead of the surface parametrization $\boldsymbol {X}^{n}$. For time discretization, we again consider a semi-implicit Euler scheme following the strategy presented in Barrett et al. (2008) and Krause & Voigt (2023). We obtain the following discrete problem.

Problem 3.2 (Discrete surface Navier–Stokes and surface update problem)

Find $({\boldsymbol {\boldsymbol {u}}_{h}},p_{h},{\mathcal {H}}_h,\boldsymbol {Y}_{h})\in [{\boldsymbol {V}_{h}}\times {{V}_{h}}\times {{V}_{h}}\times {\boldsymbol {V}_{h}}](\mathcal {S}_h^{n-1})$ such that

(3.11) \begin{gather} \left.\begin{aligned} \frac{1}{\tau} \left({{\boldsymbol{\boldsymbol{u}}}^{n}_{h}},\,{{\boldsymbol{v}_{\scriptscriptstyle{h}}}}\right)_{h} + \big({\operatorname{\boldsymbol{\nabla}}_{{\boldsymbol{\boldsymbol{w}}}_{h}^{n-1}}{\boldsymbol{\boldsymbol{u}}}^n_{h}},\,{{\boldsymbol{v}_{\scriptscriptstyle{h}}}}\big)_{h} &= \left({p_{h}^n},\,{\operatorname{div}_{\!{\boldsymbol{P}}}\!{\boldsymbol{v}_{\scriptscriptstyle{h}}}}\right)_{h} - \frac{2}{Re} \left({(\boldsymbol{\sigma}({\boldsymbol{\boldsymbol{u}}}^{n}_{h})},\,{\operatorname{\boldsymbol{\nabla}}_{\!{\boldsymbol{P}}}\!{\boldsymbol{v}_{\scriptscriptstyle{h}}}}\right)_{h} \\ &\quad + \left({\mu_{h}^n\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!\phi_{h}^n},\,{{\boldsymbol{v}_{\scriptscriptstyle{h}}}}\right)_{h} \\ &\quad + \left({ \operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}\!\left(\kappa(\phi^{n}) ({\mathcal{H}}_h-{\mathcal{H}}_0(\phi^{n}))\right) },\,{\operatorname{\boldsymbol{\nabla}}_{\!{\mathcal{S}}}({\boldsymbol{v}_{\scriptscriptstyle{h}}}\boldsymbol{\cdot}\boldsymbol{\nu}_h)}\right)_{h} \\ &\quad -\big({\kappa(\phi^{n})({\mathcal{H}}_h-{\mathcal{H}}_0(\phi^{n}))B^{n-1}},\,{{\boldsymbol{v}_{\scriptscriptstyle{h}}}\boldsymbol{\cdot}\boldsymbol{\nu}_h}\big)_{h} \\ &\quad +\,\frac{1}{\tau} \big({{\boldsymbol{\boldsymbol{u}}}^{n-1}_{h}},\,{{\boldsymbol{v}_{\scriptscriptstyle{h}}}}\big)_{h}-\left({\gamma{\boldsymbol{\boldsymbol{u}}}_h^n},\,{{\boldsymbol{v}_{\scriptscriptstyle{h}}}}\right)_{h}\!,\\ \left({\operatorname{div}_{\!{\boldsymbol{P}}}\!{\boldsymbol{\boldsymbol{u}}}^{n}_{h}},\,{{q_{\scriptscriptstyle{h}}}}\right)_{h} &= 0,\\ \frac{1}{\tau}\left({\boldsymbol{Y}_h^n \boldsymbol{\cdot} \nu_{h}},\,{{{h}_{\scriptscriptstyle{h}}}}\right)_{h} &= \left({{\boldsymbol{\boldsymbol{u}}}_h^n\boldsymbol{\cdot}\nu_{h}},\,{{{h}_{\scriptscriptstyle{h}}}}\right)_{h}\!,\\ \left({{\mathcal{H}}_h^n\boldsymbol{\nu}_h},\,{{\boldsymbol{Z}_{\scriptscriptstyle{h}}}}\right)_{h} + \left({\kappa(\phi_h^n)\operatorname{\boldsymbol{\nabla}}_{\!C}\! \boldsymbol{Y}_h^n},\,{\operatorname{\boldsymbol{\nabla}}_{\!C}\! {\boldsymbol{Z}_{\scriptscriptstyle{h}}}}\right)_{h} &={-}\big({\kappa(\phi_h^n)\operatorname{\boldsymbol{\nabla}}_{\!C}\! \boldsymbol{X}_h^{n-1}},\,{\operatorname{\boldsymbol{\nabla}}_{\!C}\! {\boldsymbol{Z}_{\scriptscriptstyle{h}}}}\big)_{h}\!, \end{aligned}\right\} \end{gather}

for all $({\boldsymbol {v}_{\scriptscriptstyle {h}}},{q_{\scriptscriptstyle {h}}},{{h}_{\scriptscriptstyle {h}}},Z_{h})\in [{\boldsymbol {V}_{h}}\times {{V}_{h}}\times {{V}_{h}}\times {\boldsymbol {V}_{h}}](\mathcal {S}_h^{n-1})$, where $B^{n-1}=(\Vert \mathcal {B}_h\Vert ^2- \frac {1}{2}\operatorname {tr} \mathcal {B}_h (\operatorname {tr}\mathcal {B}_h -{\mathcal {H}}_0(\phi ^n)) )$. (We note a typing error in Krause & Voigt (2023), which has been confirmed by the authors.)

The solutions of the discrete problems above are intermediate solutions with respect to the old surface $\mathcal {S}^{n-1}_h$. For each variable $\hat \psi \in V_h(\mathcal {S}_h^{n-1})$, we define the lift $\psi \in V_h(\mathcal {S}_h^{n})$ by the evaluation with respect to the linear reference geometry where $\psi (\boldsymbol {X}_h^n) = \hat \psi (\boldsymbol {X}_h^{n-1})$. This corresponds to a nodal interpolation with respect to the degrees of freedom. We consider the following steps in each time step:

1. (i) compute $(\hat \phi _h,\hat \mu _h)\in [{{V}_{h}}\times {{V}_{h}}](\mathcal {S}_h^{n-1})$ as intermediate solution of Problem 3.1;

2. (ii) compute $(\hat {\boldsymbol {\boldsymbol {u}}}_h,\hat p,\hat {\mathcal {H}}_h,\hat {\boldsymbol {Y}}_h)\in [{\boldsymbol {V}_{h}}\times {{V}_{h}}\times {{V}_{h}}\times {\boldsymbol {V}_{h}}](\mathcal {S}_h^{n-1})$ as intermediate solution of Problem 3.2;

3. (iii) compute the new surface $\mathcal {S}^{n}_h$ by updating its parametrization $\boldsymbol {X}_h^n=\boldsymbol {X}_h^{n-1}+\boldsymbol {Y}_h^{n}$;

4. (iv) lift the intermediate solutions on the new surface to get $(\phi _h,\mu _h)\in [{{V}_{h}}\times {{V}_{h}}](\mathcal {S}_h^{n})$ and $({\boldsymbol {\boldsymbol {u}}}_h,p,{\mathcal {H}}_h,\boldsymbol {Y}_h)\in [{\boldsymbol {V}_{h}}\times {{V}_{h}}\times {{V}_{h}}\times {\boldsymbol {V}_{h}}](\mathcal {S}_h^{n})$.

As numerical analysis results for the overall scheme are missing we justify the considered approach, essentially to perform just one iteration within each time step and the linearization of the nonlinear terms by comparing with an iterative Euler scheme. In this scheme we iterate Problems 3.1 and 3.2 within each time step until convergence and update all unknowns in each iteration. This can be considered as a fully implicit scheme. Similar fully implicit approaches for Navier–Stokes–Cahn–Hilliard equations have been considered in 2-D and 3-D settings, for example, Khanwale et al. (2020, 2023), and as the approaches mentioned above can be shown to be unconditionally energy stable. We compare the fully implicit iterative Euler scheme with our semi-implicit scheme for selected simulations, see Appendix E. The results are almost identical and indicate a reduction of computing time, which makes the following numerical studies to explore the complex interplay between surface phase composition, surface curvature and surface hydrodynamics feasible.

3.5. Implementational aspects

The discrete systems are implemented within the finite element toolbox AMDiS (Vey & Voigt 2007; Witkowski et al. 2015) based on DUNE (Alkämper et al. 2014; Sander 2020). The surface approximation is done by the Dune–CurvedGrid library developed in Praetorius & Stenger (2020). For a straightforward mesh parallelization and multiple processor computation we used the PETSc library and solved the linear system by using a direct solver. Due to the nature of the equations, which lack a maximum principle, we allow $\phi _h \in {\mathbb {R}}$ and extend all material parameters constantly for $\phi _h < -1$ and $\phi _h>1$.

The SFEM approach discussed in § 3.2 does not support topological changes of the domain. Such changes can occur due to pinch offs associated with budding events or the formation of furrows. The pinch offs are characterized by a negative Gaussian curvature $\mathcal {K}$. We use this property as an indicator to define our simulated time interval $[0,T]$, where the final simulation time $T$ is defined by $T=\inf \{ t>0\vert \mathcal {K}(t) < \mathcal {K}_0\}$ where $\mathcal {K}_0<0$ is a chosen lower bound. If this criterion is not met, we consider $T=\infty$ and solve until an equilibrium configuration is reached. Numerical methods which are able to handle topological changes, or at least able to recover after such an event in a meaningful physical state, require an implicit description of the surface. Candidates are trace finite element methods (Jankuhn & Reusken 2020), cut finite element methods, level-set methods or diffuse interface methods (Nestler et al. 2018). However, these methods are computationally more expensive and currently not applicable for the considered problem. A comparison between SFEM and trace finite element methods for Stokes flow on stationary surfaces shows a factor of $10^2$–$10^3$ lower errors for SFEM for comparable mesh sizes, depending on the considered error measure (Brandner et al. 2022a). Applications of cut finite element methods and level set methods for this problem class are not known and the first detailed numerical investigations for diffuse interface methods for vector-valued surface PDEs show strong limitations for the required accuracy for the reconstructed geometric terms (Nestler & Voigt 2023a). These results are supported in a benchmark problem for vector-valued diffusion on surfaces with large curvature (Bachini et al. 2023a). This motivates the use of SFEM and the needs to consider simulations only before a potential topological change.

4. Numerical results

4.1. Considered parameters and initial conditions

Let $\mathcal {S}$ be the unit sphere with radius $R = 1$. The velocity field is initialized by ${\boldsymbol {\boldsymbol {u}}}_0=0$ and we define two different initial configurations for the phase field. The first one is $\phi _0(\boldsymbol {x})=\tanh (x_0/\sqrt {2})$ that defines a symmetric and equal distributed phase field where the phases are separated. The second one defines a random initial condition: $\phi _1 = \min \{\max \{\hat {\phi }_1,-1\},1\}$, where $\hat {\phi }$ is given by Gaussian random field,

(4.1)\begin{equation} \hat{\phi}_1(\boldsymbol{x})={-}1.0 + \sum_{i=0}^N \exp\big({-\frac{\beta}{2}\Vert \boldsymbol{x} - \boldsymbol{x}_i\Vert^2}\big)\quad \forall \boldsymbol{x}\in\mathcal{S}, \end{equation}

where $\{\boldsymbol {x}_i\}_0^N \subset \mathcal {S}$ and $\beta =100$. This allows us to create a reproducible random initial condition. Both configurations consider an equal distributed phase field. We vary the Reynolds number $Re$ and the bending stiffness $\kappa$ but set ${\mathcal {H}}_0=\gamma =0$. Other parameters are set to $\epsilon = 0.02$, $m = 0.001$ and $\sigma =1.0$ ($\tilde {\sigma } = \frac {3}{2} \sqrt {2}$), which provides a good compromise between computational effort and physical accuracy. Numerical parameters are such that $h^3 \sim \tau$. This relation results from the properties of the discretization scheme, which can be expected to be first order in time and third order in space, which will be computationally confirmed below.

4.2. Convergence study

Due to the tight coupling between the phase field $\phi$, the surface velocity ${\boldsymbol {\boldsymbol {u}}}$, the surface pressure $p$ and the geometry $\mathcal {S}$ we expect the numerical solution to sensitively depend on the approximations. Therefore, we start by considering a convergence study. Due to a lack of analytical solutions for the full problem, we compute a numerical solution with a fine discretization and study convergence with respect to this solution. In addition, we address general properties of the solution.

We compute the errors with respect to the $L^2$ norm in space and the $L^{\infty }$ norm in time. The norms are calculated with respect to the linear reference geometry $\mathcal {S}_h^{lin}$, where the evaluation of the norms on different surfaces $\mathcal {S},\mathcal {S}_h,\mathcal {S}_h^{lin}$ leads to equivalent norms with the same order of convergence (Demlow 2009). We define the following errors:

(4.2) \begin{equation} \left.\begin{aligned} e_{\boldsymbol{X}}&=\Vert \boldsymbol{X}_h - \boldsymbol{X} \Vert_{L^{\infty}(L^2(\mathcal{S}_h^{lin}))}\!\kern 0.09em , \\ e_{{\boldsymbol{\boldsymbol{u}}}}&=\Vert {\boldsymbol{\boldsymbol{u}}}_h(\boldsymbol{X}_h) - {\boldsymbol{\boldsymbol{u}}}(\boldsymbol{X}) \Vert_{L^{\infty}(L^2(\mathcal{S}_h^{lin}))}\! \kern 0.09em , \\ e_{\phi}&=\Vert \phi_h(\boldsymbol{X}_h) - \phi(\boldsymbol{X}) \Vert_{L^{\infty}(L^2(\mathcal{S}_h^{lin}))}\! \kern 0.09em , \\ e_A &= \vert \Vert\mathcal{S}_h(t)\vert - \vert \mathcal{S}_h(0) \vert \vert_{L^{\infty}}\!\kern 0.09em , \\ e_{\operatorname{div}_{\!{\boldsymbol{P}}}\!}&=\Vert \operatorname{div}_{\!{\boldsymbol{P}}}\!{\boldsymbol{\boldsymbol{u}}}_h \Vert_{L^{\infty}(L^2(\mathcal{S}_h))} \!\kern 0.09em , \\ e_{{\mathcal{F}}}&=\vert ({{\mathcal{F}}_{K}}_h+{\mathcal{F}}_h) - ({{\mathcal{F}}_{K}}+{\mathcal{F}}) \vert_{L^{\infty}} \!\kern 0.09em , \end{aligned}\right\} \end{equation}

with discrete solutions $\boldsymbol {X}_h,{\boldsymbol {\boldsymbol {u}}}_h,\phi _h$ with varying $h$ and $\boldsymbol {X},{\boldsymbol {\boldsymbol {u}}},\phi$ the discrete reference solutions with the finest mesh size. Additionally, we measure the area conservation error $e_A$, the local inextensibility error $e_{\operatorname {div}_{\!{\boldsymbol {P}}}\!}$ and the error of the total energy $e_{{\mathcal {F}}}$, where again ${{\mathcal {F}}_{K}}$ and ${\mathcal {F}}$ denote the discrete reference solution on the finest mesh size. The simulations are done for the random initial condition $\phi _1$ with $Re = 1.0$ and $\kappa _1 = \kappa _2 = \kappa = 0.02$.

The results are shown in figure 1. They indicate third-order convergence for the surface error $e_{\boldsymbol {X}}$ and the error of the phase field $e_{\phi }$, which corresponds to the results in Praetorius & Stenger (2020) and Demlow (2009) for simple mean curvature flow and scalar-valued surface PDEs. The order of convergence for the phase field $e_{\phi }$ also agrees with the one obtained for the corresponding problem on stationary surfaces considered in Bachini et al. (2023b). For the error of the velocity field $e_{{\boldsymbol {\boldsymbol {u}}}}$ the results indicate second-order convergence in accordance with the results in Krause & Voigt (2023) for one-component fluid deformable surfaces. In Brandner et al. (2022a), third-order convergence of $e_{{\boldsymbol {\boldsymbol {u}}}}$ is shown for the tangential flow of the surface Stokes equations on stationary surfaces. However, this requires a higher-order approximation of the normal vector (Hansbo et al. 2020; Hardering & Praetorius 2022), which is not fulfilled in our approach. Furthermore, it remains open if this increased order also emerges for the surface Navier–Stokes equations and on evolving surfaces. For the area conservation error $e_A$ and the inextensibility error $e_{\operatorname {div}_{\!{\boldsymbol {P}}}\!}$ we get second-order convergence with respect to $h$. Both are related to each other. While $e_A = 0$, the same order for $e_{\operatorname {div}_{\!{\boldsymbol {P}}}\!}$ has been obtained in Bachini et al. (2023b) for a stationary surface. The convergence error of the total energy indicates third order. This might be due to the dominating effect of Ginzburg–Landau energy as the expected order for the underlying Helfrich model would be lower (Dziuk 2008). However, overall we see experimentally the expected optimal order of convergence for the full problem.

Figure 1. Convergence study for two-phase fluid deformable surfaces with respect to the mesh size $h^3 \sim \tau$ for $\boldsymbol {X}_h,{\boldsymbol {\boldsymbol {u}}}_h,\phi _h$ (a–c); and error for area $A$, divergence ${\rm div}_{\boldsymbol {P}} \boldsymbol {u}_h$ and total energy ${\mathcal {F}}_K + {\mathcal {F}}$ (d–f).

4.3. Phase separation, bulging and induced flow

To demonstrate the strong interplay between composition, curvature and hydrodynamics we consider $\phi _1$ as an initial condition. The other parameters are $Re = 1.0$ and either $\kappa _1 = \kappa _2 = \kappa = 0.02$ or $\kappa _1 = 0.02$ and $\kappa _2 = 0.5$. In the last case, the red coloured phase has a lower bending stiffness and is therefore expected to be guided to or initiate regions with higher mean curvature, while the blue coloured phase has a larger bending stiffness and can be expected to prefer regions of lower mean curvature. The evolution is shown in figure 2(a,b), respectively. Both cases show the composition and the tangential velocity for selected time instances. In both cases red and blue islands form and a wavy interface between larger red and blue regions is established. Some coarsening events can be spotted and the wavy interface flattens over time. However, the main evolution is in the shape, which strongly deforms and forms bulges. Especially for circular islands the interface length is reduced by bulging leading to strong deformation from the sphere. These shapes change, but also the coarsening events induce flow. The simulations are only shown for a relatively short period of time and are terminated before a potential topological change would happen. Differences between figure 2(a) and figure 2(b) are also visible. While in figure 2(a) the red and blue phases evolve similarly, in figure 2(b) the blue phase, the one associated with a larger bending stiffness, forms less curved regions. Especially, islands of this phase do not bulge out. Instead, they become relatively flat. The behaviour of the simulation shown in figure 2 is quantified in figure 3 by different energies. In the current setting, the evolution is dominated by ${{\mathcal {F}}_{GL}}$, which is reduced by coarsening but also by bulging, which is associated with an increase in ${\mathcal {F}}_{H}$. Hydrodynamics seems to play a minor role in the evolution. Here ${{\mathcal {F}}_{K}}$ is one to two orders of magnitude lower than ${{\mathcal {F}}_{GL}}$. The strong increase of ${{\mathcal {F}}_{K}}$ and ${\mathcal {F}}_{H}$ at the beginning is associated with the increased driving force resulting from the phase-dependent bending stiffness. The spatial average over the mean curvature for the red and the blue phase, i.e.

(4.3a,b)\begin{equation} \bar{{\mathcal{H}}}_1 = \int_{\{\phi<0\}} {\mathcal{H}}\quad \text{and} \quad \bar{{\mathcal{H}}}_2 = \int_{\{\phi>0\}} {\mathcal{H}}, \end{equation}

is shown in figure 3(c). The deviation between $\bar {{\mathcal {H}}}_1$ and $\bar {{\mathcal {H}}}_2$ for the setting of figure 2(a) results from asymmetries in the initial condition. For the setting of figure 2(b), these quantities show strong differences between the different phases. The blue phase on average forms flat regions, while the red phase forms strongly curved bulges.

Figure 2. (a,b) Snapshots of the relaxation of the two-component fluid deformable surface with random initial condition $\phi _1$ and $Re=1.0$ for $t=0,0.3,0.8,1.1$ (from left to right), with constant bending stiffness $\kappa _1 = \kappa _2 = \kappa =0.02$ in (a) and phase-depended bending stiffness $\kappa _1 = 0.02$ and $\kappa _2 = 0.5$ in (b). In (b) the red coloured phase is less stiff than the blue coloured phase. Here, in each of (a,b), is shown the (top) phase field $\phi$; (bottom) tangential velocity ${\boldsymbol {P}}{\boldsymbol {\boldsymbol {u}}}$. The flow is visualized by a LIC (Line Integral Convolution) filter and colour coding represents the magnitude. Corresponding movies are provided in the Supplementary data available at https://doi.org/10.1017/jfm.2023.943.

Figure 3. (a,b) The energies ${{\mathcal {F}}_{K}}$, ${{\mathcal {F}}_{GL}}$, ${\mathcal {F}}_{H}$ and ${\mathcal {F}}_K+{\mathcal {F}}$ over time where (a) corresponds to figure 2(a) and (b) corresponds to figure 2(b). The time instances are highlighted in the plots. (c) Averaged mean curvature $\bar {{\mathcal {H}}}_{1,2}$ for the different phases with respect to the simulations done in figure 2(a,b), the colour corresponds to the coloured phases.

4.4. Variation of parameters

We now explore the influence of the parameters in more detail. To reduce the number of parameters we only consider the situation of a constant bending stiffness $\kappa _1 = \kappa _2 = \kappa$, set $\kappa = \{0.02, 0.1, 0.5\}$, and vary the Reynolds number with $Re = \{0.1, 1, 3\}$. We consider both initial conditions $\phi _0$ and $\phi _1$. The behaviour is demonstrated for the symmetric initial phase field $\phi _0$ in figure 4 and for the random initial phase field $\phi _1$ in figure 5. We only show the final configuration, this is either the equilibrium configuration or the configuration at the critical time $T$ before a potential pinch-off. A table with critical final times is presented in table 1. In the first case ($\phi _0$) the final configuration is mainly determined by the interplay of the Helfrich energy ${\mathcal {F}}_{H}$ and the Ginsburg–Landau energy ${{\mathcal {F}}_{GL}}$. For $\kappa = 0.5$, ${\mathcal {F}}_{H}$ dominates. Any deformation of the initial sphere requires us to increase the Helfrich energy. Only small deformations are possible. The interface length, and thus ${{\mathcal {F}}_{GL}}$, is reduced by deforming the sphere to an ellipsoid-like shape with the interface placed along the short axis. Further reduction of the interface length either requires to increase the mean curvature at the tips of the long axis of this shape or to form a furrow alone the interface. Both effects strongly increase ${\mathcal {F}}_{H}$. This is only realized for settings with lower $\kappa$ and eventually will lead to a pinch-off. Hydrodynamics seems to have no significant effect in these evolutions. However, a closer look at the consider critical times $T$ in table 1 shows the opposite. While the final configurations look similar, the time to reach these configurations is strongly influenced by hydrodynamics. For low bending stiffness, $\kappa = \{0.1, 0.02\}$, the critical time is drastically reduced if the Reynolds number $Re$ is increased.

Figure 4. Final configuration obtained for different parameters and initial condition $\phi _0$. Either the simulation reaches the equilibrium configuration or the criteria for a potential pinch-off is reached. The corresponding critical times $T$ are shown in table 1.

Figure 5. Final configuration obtained for different parameters and initial condition $\phi _1$. Either the simulation reaches the equilibrium configuration or the criteria for a potential pinch-off is reached. The corresponding critical times $T$ are shown in table 1.

Table 1. Critical times $T$ for the different Reynolds numbers $Re$ and bending stiffness $\kappa$ for the symmetric initial value $\phi _0$ in (a) and the random initial value $\phi _1$ in (b).

For the random initial condition $\phi _1$ the results are shown in figure 5. Here, we have an interplay between coarsening, shape deformation, bulging and furrow formation. For large bending stiffness $\kappa = 0.5$, where strong shape deformations are suppressed, the same final configuration as in figure 4 is reached, meaning an ellipsoidal-like shape with the interface placed along the short axis. This changes for lower $\kappa$. For $\kappa =0.1$, we observe partial coarsening, bulging of islands and the formation of a furrow, which eventually leads to a pinch-off. The configurations also change with hydrodynamics. For $Re = 3$ the potential pinch-off happens earlier. Several islands are still present in the red as well as the blue phase. Decreasing $Re$ also leads to possible pinch-offs but at a later coarsening state. There are fewer islands present. While it is known in flat space that hydrodynamics can enhance coarsening and this is also demonstrated on stationary surfaces (Bachini et al. 2023b), these results indicate that this also holds on evolving surfaces and that hydrodynamics enhances furrow formation and potential pinch-off. This is confirmed in table 1, which again shows a drastic reduction of the critical time $T$. For $\kappa = 0.02$ the evolution is dominated by bulging. As discussed in the previous section bulging is associated with large local absolute mean curvature values but allows us to reduce the interface length. The critical time results from potential pinch-offs of the bulges. Coarsening is only a minor effect. Again hydrodynamic drastically enhances the evolution. Increasing $Re$ again drastically reduces the time for potential pinch-off, see table 1.

We also would like to comment on the axisymmetric ellipsoidal-like configurations reached for $\kappa = 0.5$. With friction $\gamma > 0$ these configurations would be characterized by zero tangential velocity $\boldsymbol {P} \boldsymbol {u} = 0$ and they would be independent on the initial condition and the Reynolds number $Re$. This equilibrium state coincides with the corresponding result for the reduced model without hydrodynamics, see Appendix D. However, without friction $\gamma = 0$ the dissipation potential in the reached axisymmetric state is invariant under surface rigid body motion. Such configurations have been explored for one-component fluid deformable surfaces (Reuther et al. 2020; Krause & Voigt 2023; Nestler & Voigt 2023b; Olshanskii 2023). Indeed the shown configurations in figures 4 and 5 for $\kappa = 0.5$ undergo slight rigid body motions and as the resulting forces interact with the shape and the phase composition also slightly differ. This result can be viewed as an extension of the phenomenon of rotating equilibrium states from one-component to two-component fluid deformable surfaces. However, as discussed in Nestler & Voigt (2023b) any perturbation in shape or velocity that destroys the axisymmetric configurations leads to dissipation and we therefore can expect for $t \to \infty$ to reach the mentioned equilibrium state with $\boldsymbol {P} \boldsymbol {u} = 0$.

However, the main focus of this paper is on the dynamics. One could ask about the behaviour if the Reynolds number $Re$ is further reduced. This can be considered in the Stokes limit. This limit cannot be directly derived by the Lagrange–d’Alembert principle and requires some further thoughts on the fluid–solid duality of the fluid deformable surface model, which will be explored in future research. However, we can already remark that the dynamics for two-phase fluid deformable surfaces, also for low $Re$, strongly differ from previous models neglecting surface viscosity. The derivation by the Lagrange–d’Alembert principle accounts for normal and tangential variations, which differs from various previous approaches, see, for example, Elliott & Stinner (2010a). Qualitative differences already result from the tangential flow induced by the phase boundary, see Appendix D for corresponding results of the model proposed in Elliott & Stinner (2010a) with additional local inextensibility constraint. These differences are further enhanced by hydrodynamics.

5. Conclusions

While the literature on coexisting fluid domains in model systems for biomembranes and their dynamics is rich, the influence of surface viscosity has not been discussed in this context. We have filled this gap and derived a thermodynamically consistent model that accounts for surface hydrodynamics. The derivation of the model by a Lagrange–d’Alembert principle is applicable in general. Simplifications of the derived two-phase fluid deformable surface model lead to known models in the literature, and this further confirms the validity of the approach.

By combining numerical approaches for fluid deformable surfaces (Krause & Voigt 2023) and two-phase flow models on stationary surfaces (Bachini et al. 2023b), we obtain a numerical approach for the full model, which is demonstrated to converge with expected optimal order. This provides the basis for a detailed investigation of the strong interplay of surface phase composition, surface curvature and surface hydrodynamics. Depending on the material parameters, the line tension $\sigma$, the Reynolds number $Re$ and the bending stiffness $\kappa$, the evolution is dominated by coarsening or shape evolution. Both phenomena are shown to be strongly enhanced by hydrodynamics. In situations where the line tension and the bending stiffness are compatible, the interface length is reduced by the formation of a furrow or the formation of bulges. Both potentially lead to pinch-offs. The enhanced evolution towards such topological changes driven by hydrodynamics has various biological implications. In the context of a furrow formation and its subsequent shrinkage, this can be associated with cell division (Mietke et al. 2019); in the context of bulges with pinch-offs, this can be associated with endocytosis and exocytosis (Kaksonen & Roux 2018; Al-Izzi et al. 2020). While quantitative comparisons with experimental data in these contexts require further studies, the mathematical model and the numerical algorithm to solve it with the required accuracy are provided.