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 (Reference Jankuhn, Olshanskii and Reusken2018) and the componentwise surface gradient $\operatorname {\boldsymbol {\nabla }}_{\!C}$ as in Nitschke, Sadik & Voigt (Reference Nitschke, Sadik and Voigt2022), 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. (Reference Reuther, Nitschke and Voigt2020). 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 (Reference Nitschke and Voigt2022). 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 (Reference Marsden and Ratiu2013). 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 (Reference Marsden and West2001), Marsden & Ratiu (Reference Marsden and Ratiu2013) and Hairer et al. (Reference Hairer, Hochbruck, Iserles and Lubich2006) 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. (Reference Torres-Sánchez, Millán and Arroyo2019) 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. Reference Garcke, Kampmann, Rätz and Röger2016; Elliott, Hatcher & Stinner Reference Elliott, Hatcher and Stinner2022; Benes et al. Reference Benes, Kolar, Sischka and Voigt2023), $\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. (Reference Fonda, Rinaldin, Kraft and Giomi2018) and Bachini et al. (Reference Bachini, Krause and Voigt2023b), we consider a diffuse interface approximation of the Jülicher–Lipowski model (Jülicher & Lipowsky Reference Jülicher and Lipowsky1996),

(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 (Reference Elliott and Stinner2010b) and Bachini et al. (Reference Bachini, Krause and Voigt2023b),

(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. (Reference Reuther, Nitschke and Voigt2020) 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. (Reference Torres-Sánchez, Millán and Arroyo2019) 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. (Reference Jankuhn, Olshanskii and Reusken2018). 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. (Reference Magaletti, Picano, Chinappi, Marino and Casciola2013) and consider this parameter to be constant. For alternative approaches, we refer to Abels (Reference Abels2009). 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. (Reference Torres-Sánchez, Millán and Arroyo2019) and Reuther et al. (Reference Reuther, Nitschke and Voigt2020). 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 (Reference Marsden and West2001), Udwadia & Kalaba (Reference Udwadia and Kalaba2002) and Izadi & Sanyal (Reference Izadi and Sanyal2014). 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 (Reference Lowengrub and Truskinovsky1998), Abels, Garcke & Grün (Reference Abels, Garcke and Grün2012), Aland & Voigt (Reference Aland and Voigt2012) and Ten Eikelder et al. (Reference Ten Eikelder, Van Der Zee, Akkerman and Schillinger Chillinger2023). 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. (Reference Reuther, Nitschke and Voigt2020), Krause & Voigt (Reference Krause and Voigt2023) and in the Stokes limit in Torres-Sánchez et al. (Reference Torres-Sánchez, Millán and Arroyo2019). For further numerical and analytical approaches under additional symmetry assumptions we refer to Al-Izzi et al. (Reference Al-Izzi, Sens and Turner2020), where a linear stability analysis of a tube geometry is considered, and to Olshanskii (Reference Olshanskii2023), where potential rotational symmetric equilibrium configurations are addressed. For a comparison of different derivations of this model we refer to Reuther & Voigt (Reference Reuther and Voigt2015, Reference Reuther and Voigt2018a) and Brandner, Reusken & Schwering (Reference Brandner, Reusken and Schwering2022b).

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. (Reference Jankuhn, Olshanskii and Reusken2018) and used in Reuther et al. (Reference Reuther, Nitschke and Voigt2020). 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. (Reference Bachini, Krause and Voigt2023b), where it is used to study the influence of surface hydrodynamics on coarsening in multicomponent scaffolded lipid vesicles, see Fonda et al. (Reference Fonda, Rinaldin, Kraft and Giomi2018), Fonda et al. (Reference Fonda, Rinaldin, Kraft and Giomi2019) and Rinaldin et al. (Reference Rinaldin, Fonda, Giomi and Kraft2020). Without the bending terms, the system has already been addressed in Nitschke et al. (Reference Nitschke, Voigt and Wensch2012), Ambrus et al. (Reference Ambrus, Busuioc, Wagner, Paillusson and Kusumaatmaja2019) and Olshanskii et al. (Reference Olshanskii, Palzhanov and Quaini2022).

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. (Reference Nitschke, Voigt and Wensch2012) and Reuther & Voigt (Reference Reuther and Voigt2015, Reference Reuther and Voigt2018a) (and in the Stokes limit in Olshanskii et al. (Reference Olshanskii, Quaini, Reusken and Yushutin2018), Bonito, Demlow & Licht (Reference Bonito, Demlow and Licht2020) and Brandner & Reusken (Reference Brandner and Reusken2020)), and for general surfaces in Nitschke, Reuther & Voigt (Reference Nitschke, Reuther and Voigt2017), Reusken (Reference Reusken2018), Fries (Reference Fries2018), Reuther et al. (Reference Reuther, Nitschke and Voigt2020), Reuther & Voigt (Reference Reuther and Voigt2018b) and Lederer, Lehrenfeld & Schoeberl (Reference Lederer, Lehrenfeld and Schoeberl2020) (and in the Stokes limit in Olshanskii & Yushutin (Reference Olshanskii and Yushutin2019)).

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. (Reference Haußer, Marth, Li, Lowengrub, Rätz and Voigt2013), 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 (Reference Wang and Du2008) and Elliott & Stinner (Reference Elliott and Stinner2010a), 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 (Reference Elliott and Stinner2010a) 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.1. Mesh movement

By following the work in Krause & Voigt (Reference Krause and Voigt2023), we combine the system in Problem 2.2 with a mesh redistribution approach introduced in detail in Barrett, Garcke & Nürnberg (Reference Barrett, Garcke and Nürnberg2008). 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 (Reference Praetorius and Stenger2020) and Brandner et al. (Reference Brandner, Jankuhn, Praetorius, Reusken and Voigt2022a). 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. (Reference Bachini, Krause and Voigt2023b), 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. (Reference Demont, van Zwieten, Diddens and van Brummelen2022) for a detailed analysis. Similar properties have been found in Bachini et al. (Reference Bachini, Krause and Voigt2023b) on surfaces. In contrast to Bachini et al. (Reference Bachini, Krause and Voigt2023b), where the surface was stationary, here the surface evolves and we adapt the numerical schemes used in Krause & Voigt (Reference Krause and Voigt2023) 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 (Reference Shen and Yang2010a,Reference Shen and Yangb). These methods and more an advanced scalar auxiliary variable (SAV) approaches, for example, Zhu et al. (Reference Zhu, Chen, Yao and Sun2019), 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. (Reference Barrett, Garcke and Nürnberg2008) and Krause & Voigt (Reference Krause and Voigt2023). 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 (Reference Krause and Voigt2023), 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. (Reference Khanwale, Lofquist, Sundar, Rossmanith and Ganapathysubramanian2020, Reference Khanwale, Saurabh, Ishii, Sundar, Rossmanith and Ganapathysubramanian2023), 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 Reference Vey and Voigt2007; Witkowski et al. Reference Witkowski, Ling, Praetorius and Voigt2015) based on DUNE (Alkämper et al. Reference Alkämper, Dedner, Klöfkorn and Nolte2014; Sander Reference Sander2020). The surface approximation is done by the Dune–CurvedGrid library developed in Praetorius & Stenger (Reference Praetorius and Stenger2020). 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 Reference Jankuhn and Reusken2020), cut finite element methods, level-set methods or diffuse interface methods (Nestler et al. Reference Nestler, Nitschke, Praetorius and Voigt2018). 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. Reference Brandner, Jankuhn, Praetorius, Reusken and Voigt2022a). 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 Reference Nestler and Voigt2023a). These results are supported in a benchmark problem for vector-valued diffusion on surfaces with large curvature (Bachini et al. Reference Bachini, Brandner, Jankuhn, Nestler, Praetorius, Reusken and Voigt2023a). This motivates the use of SFEM and the needs to consider simulations only before a potential topological change.

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 Reference Demlow2009). 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 (Reference Praetorius and Stenger2020) and Demlow (Reference Demlow2009) 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. (Reference Bachini, Krause and Voigt2023b). 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 (Reference Krause and Voigt2023) for one-component fluid deformable surfaces. In Brandner et al. (Reference Brandner, Jankuhn, Praetorius, Reusken and Voigt2022a), 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. Reference Hansbo, Larson and Larsson2020; Hardering & Praetorius Reference Hardering and Praetorius2022), 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. (Reference Bachini, Krause and Voigt2023b) 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 Reference Dziuk2008). 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. Reference Bachini, Krause and Voigt2023b), 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. Reference Reuther, Nitschke and Voigt2020; Krause & Voigt Reference Krause and Voigt2023; Nestler & Voigt Reference Nestler and Voigt2023b; Olshanskii Reference Olshanskii2023). 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 (Reference Nestler and Voigt2023b) 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 (Reference Elliott and Stinner2010a). 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 (Reference Elliott and Stinner2010a) with additional local inextensibility constraint. These differences are further enhanced by hydrodynamics.