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.


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 highresolution fluorescence imaging [11] a strong correlation between domain composition and local membrane curvature can be established.These findings have supported earlier membrane models [42,41,45] and initiated a wealth of theoretical, numerical, and experimental studies to understand the relation between composition and curvature [78,11,10,79,81,48,24,25,31,32,85].This relation has strong biological implications [52].We refer to [78,19,47] for reviews on the subject.While most studies address multicomponent giant unilaminar vesicles (GUV) and consider phase separation and coarsening on spherical shapes, more recently multicomponent scaffolded lipid vesicles (SLV) 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 [28,29,72].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 [11,48,24], 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 modeled by the Jülicher-Lipowsky model [42] or appropriate phase field approximations of it [48,36].These models essentially extend the classical Helfrich model to multiple components.Instead of a simple  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 [75,84].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 [26] 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 [27,17].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 [9,83,69], 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 [73], 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 2D on spatial scales smaller than the Saffmann-Delbrück number [37].We follow this simplification and consider only surface two-phase flow problems.On stationary surfaces these problems are addressed in [60,5,62,7] and essentially reveal an enhanced coarsening rate due to hydrodynamic effects, similar to the situation in flat space [27,17].To consider evolving surfaces under the influence of surface viscosity requires models for so-called fluid deformable surfaces.They exhibit a solid-fluid duality, while they store elastic energy when stretched or bent, as solid shells, they flow as viscous two-dimensional 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 [76,67,44].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) [21,55] and builds on previous developments for one-component fluid deformable surfaces [44] and surface Navier-Stokes-Cahn-Hilliard-like models on stationary surfaces [7].The algorithm is used to demonstrate the strong interplay between composition, curvature and hydrodynamics and its implications for bulging and budding processes.Briefly, surface hydrodynamics enhances the onset of such shape deformations and possible resulting topological changes, which has strong biological implications, e.g., in the case of endocytosis and exocytosis [43,3].
The paper is structured as follows.In Section 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 Section 3 we discuss the considered numerical approach.In Section 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 Section 5 we draw conclusions.

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.

Notation
We consider a time dependent smooth and oriented surface S = S() without boundary.Related to S, we denote the outward pointing surface normal , the surface projection  =  −  ⊗ , the shape operator B = − ∇  , and the mean curvature H = tr B. Note that, under these definitions the unit sphere has negative mean curvature H = −2.Let  be a continuously differentiable scalar field,  a continuously differentiable R 3 -vector field, and  a continuously differentiable R 3×3 -tensor field defined on S. We define the surface tangential gradient ∇  as in [39] and the componentwise surface gradient ∇  as in [58], namely: where   ,   and   are arbitrary smooth extensions of ,  and  in the normal direction and ∇ is the gradient of the embedding space R 3 .The fields ∇  , ∇   are purely tangential vector and tensor fields, respectively.We define the corresponding divergence operators for a vector field  and a tensor field  by: where tr is the trace operator.The divergence of the tensor , div  , leads to a non-tangential vector field even if  is a tangential tensor field.Let ∇ S be the gradient with respect to the covariant derivative on S, as used in [67].This operator is defined for scalar fields and tangential vector fields and relates to the tangential operators by ∇   = ∇ S  and div   = div S () − ( • )H , respectively.We further clarify the relation between the above operators in Appendix A.
The surface S is given by a parametrization .The material on the surface is described by a material parametrization   , as in [59].Both parametrizations relate to each other by     •  =    • , 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      and the relative material velocity by  :=  −   .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.

Model derivation by Lagrange-D'Alembert principle
The Lagrange-D'Alembert principle has been explained in [50].The approach is a combination of the Lagrange and the Onsager variational principles.The concept of the Lagrange principle has been introduced in [51,50,33] and it models the interplay of kinetic and potential energies under total energy conservation.On the opposite, the concept of the Onsager variational principle models dissipative systems, where potential energy decreases under a dissipation potential, e.g. 2gradient flows.Within our context, the Onsager principle is used in [76] to derive fluid deformable surfaces with the surface Stokes model.
We consider a density function  that describes the two-phases on the surface S, e.g., liquidordered and disordered phases, where one phase is represented by  = 1, the other one by  = −1, and we assume a mixture of both if  ∈ (−1, 1).We consider a conserved evolution for  with respect to the following energies.
The second energy we consider accounts for bending properties.As in [28,7], we consider a diffuse interface approximation of the Jülicher-Lipowski model [42]: where () is the bending stiffness and H 0 () the spontaneous curvature.We neglect additional contributions due to the Gaussian curvature of S. This is justified as long as the Gaussian bending stiffness is independent on the phase  and no topological changes occur.To get a continuously differentiable dependency on  we consider an interpolation function as in [24] and [7]: for  ∈ {, H 0 } and  1 ,  2 the material parameter of  or H 0 in the separated phases.Together, the energies F  and F  define the potential energy F = F  + F  of the system.We define the surface material velocity by  and the kinetic energy as in [67] by: with  the surface density.For simplicity, we assume  to be constant.The Lagrangian L = F  −F is defined as the difference between the kinetic and the potential energy.
We next consider the various sources of dissipation.As in [76] we define the dissipation potential of the viscous stress: where  denotes the viscosity and () = 1 2 (∇   + (∇  )  ) is the rate of deformation tensor as considered in [39].For simplicity, we assume  to be constant.In addition, we consider the friction with the surrounded material, which is modeled by: where  ≥ 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 give by: with  > 0 a resistance associated with  =    + ∇   with respect to  −1 norm.Thereby,  denotes the material time derivative of  and ∇   = (∇ S , ), see Appendix B.6.The parameter  plays the role of a mobility.Here, we follow the approach of [49] and consider this parameter to be constant.For alternative approaches, we refer to [1].We now define the dissipation potential by D = D  + D  + D  .Moreover, we add to L and D the following constraints.We here only assume local inextensibility of the material as in [76,67].This is incorporated by the Lagrange-function , which serves as the surface pressure and is related to the surface tension.It enters in the constraint by: This constraint induces a conservation of surface area |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 [51,77,38].Here we consider this principle for space depended functions, for which it provides an elegant way to derive thermodynamically consistent models.It reads: for all test functions  ∈ R 3 | S and ,  ∈  0 S. Here, the symbol D • • denotes the  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 1 Find (, , , ) such that: , with tangential and normal bending forces defined by: and  the relative material velocity as explained above.
This problem provides a tight coupling between the phase field , the surface velocity , the surface pressure , and the geometry S. Exploring the main implications of these couplings is one of the goals of this paper.

Remark 1 (Thermodynamic consistency)
We refer to Appendix B.5 to show that the relation holds for the model in Problem 1.This demonstrates thermodynamic consistency.
Considering a characteristic length and a characteristic velocity, Problem 1 can be formulated in non-dimensional form, see Appendix C for details.Keeping the same notation also in nondimensional form, we obtain: Problem 2 Find (, , , ) such that: where Re denotes the Reynolds number and   and   are defined as in Problem 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.

Model simplifications
The model in Problem 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 , the viscosity , and the friction coefficient , are constant.It is possible to extend the model and consider phase-dependent parameters by following the approach described in [2].In the following, we link our model to simplified models already discussed in the literature.

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  ≡ const, the system simplifies to: where  =  −    as before, and the normal bending forces reduce to: Note that, bending stiffness and spontaneous curvature are now constant parameters.In the case of H 0 = 0 and  = 0, this model has been introduced and simulated in [67,44] and in the Stokes-limit in [76].For further numerical and analytical approaches under additional symmetry assumptions we refer to [3], where a linear stability analysis of a tube-geometry is considered, and to [61], where potential rotational symmetric equilibrium configurations are addressed.For a comparison of different derivations of this model we refer to [68,70] and [16].

Two component fluid on a stationary surface
If we assume a stationary surface S, i.e.  •  = 0, the model in Problem 2 restricts to a pure tangential problem.We follow the approach of the directional splitting as shown in detail in [39] and used in [67].We note that   =  and   =  •  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: where the relative material velocity simplifies to the Eulerian case  =   .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  = 0, the model has been introduced and discussed in [7], where it is used to study the influence of surface hydrodynamics on coarsening in multicomponent scaffolded lipid vesicles (SLV), see [28,29,72].Without the bending terms, the system has already been addressed in [60,5,62].

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: where again  =   .These equations have been solved numerically for simply connected surfaces in [60,68,70] (and in the Stokes limit in [64,15,13]), and for general surfaces in [57,66,30,67,71,46] (and in the Stokes limit in [63]).

Two component surface without fluid behavior
A larger literature exists for the evolution of two-component surfaces if surface hydrodynamics is neglected.These models follow as the overdamped limit of Problem 2. An explicit derivation is shown in Appendix D. The resulting model reads: where b = b •.In this model, p serves as a Lagrange-function to ensure the local inextensibility constraint.The model is similar to a model discussed in [36], 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 [81] and [23].

Numerical discretization
We consider a surface finite element method (SFEM) [21,55] to solve the highly nonlinear set of geometric and surface partial differential equations (PDEs) in Problem 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 vectorvalued surface PDEs ensuring convergence of surface vector-Laplace and surface Stokes problems on stationary surfaces [34,35] and reproduces the same expected optimal order of convergences for the full two component fluid deformable surface model.

Mesh movement
By following the work in [44], we combine the system in Problem 2 with a mesh redistribution approach introduced in detail in [8] .We consider the initial surface given  (0) =  0 and an additional equation for the time evolution of the parametrization: This approach generates a tangential mesh movement that maintains the shape regularity and additionally provides an implicit representation of the mean curvature H .

Surface Approximation
We assume that the smooth surface S is approximated by a discrete -th order approximation S ℎ .Let S  ℎ be a piecewise linear reference surface given by shape regular triangulation We define  as bijective map  : S  ℎ → S such that S = ∪   =1  (  ).The construction of such maps is discussed in [65] and [14].We get a -th order approximation of  by the -th order interpolation S  ℎ =   ℎ ( ), which defines a higher order triangulation such that S  ℎ = ∪   =1   ℎ (  ).We use each geometrical quantity like the normal vector  ℎ , the shape operator B ℎ , and the inner products (•, •) ℎ with respect the S  ℎ where we will drop the index  in the following.We denote the size of the grid by ℎ, i.e. the longest edge of the mesh.

Discrete function spaces
We define the discrete function spaces for scalar function by: with P   the space polynomials, where we set the element order as   = .We define   (S ℎ ) = [  (S ℎ )] 3 as space of discrete vector fields.We discretize and  ℎ ∈  1 (S ℎ ).This corresponds to a Taylor-Hood element for the pair of velocity and pressure.

Discrete model
Following the strategy in [7], we separate the surface phase field and surface Navier-Stokes equations by an operator splitting approach.In contrast with [7], where the surface was stationary, here the surface evolves and we adapt the numerical schemes used in [44] for one-component systems.
Let {  }   be time interval discretization, with   =  and  > 0 be the time steps.We first solve the surface phase field problem, eqs.( 10)- (11), and then solve the surface Navier-Stokes, eqs.( 12) and ( 13), together with eqs.( 14)-( 15) for mesh regularization.The inner product (• , •) ℎ 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 S −1 ℎ = S ℎ ( −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 (Discrete surface Cahn-Hilliard problem)
).The surface Navier-Stokes equations and the mesh redistribution are considered together.Moreover, we define a discrete surface update variable   ℎ =   ℎ −  −1 ℎ , which is considered as unknown instead of the surface parametrization   .For time discretization, we again consider a semi-implicit Euler scheme following the strategy presented in [8,44].We obtain the following discrete problem.

Problem 4 (Discrete surface Navier-Stokes and surface update problem)
The solutions of the discrete problems above are intermediate solutions with respect to the old surface S −1 ℎ .For each variable ψ ∈  ℎ (S −1 ℎ ), we define the lift  ∈  ℎ (S  ℎ ) by the evaluation with respect to the linear reference geometry where (   ℎ ) = ψ(  −1 ℎ ).This corresponds to a nodal interpolation with respect to the degrees of freedom.We consider the following steps in each time step:

Compute the new surface S 𝑛
ℎ by updating its parametrization

Implementational aspects
The discrete systems are implemented within the finite element toolbox AMDiS [80,82] based on DUNE [74,4].The surface approximation is done by the Dune-CurvedGrid library developed in [65].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  ℎ ∈ R and extend all material parameters constantly for  ℎ < −1 and  ℎ > 1.
The SFEM approach discussed in Section 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 K.We use this property as indicator to define our simulated time interval [0, ], where the final simulation time  is defined by  = inf{ > 0|K () < K 0 } where K 0 < 0 is a chosen lower bound.If this criteria is not met, we consider  = ∞ 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 (TraceFEM) [40], cut finite element methods (CutFEM), level-set methods or diffuse interface methods [54].However, these methods are computationally more expensive and currently not applicable for the considered problem.A comparison between SFEM and TraceFEM 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 [14].Applications of CutFEM 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 [56].These results are supported in a benchmark problem for vector-valued diffusion on surfaces with large curvature [6].This motivates the use of SFEM and the needs to consider simulations only before a potential topological change.

Considered parameters and initial conditions
Let S be the unit sphere with radius  = 1.The velocity field is initialized by  0 = 0 and we define two different initial configurations for the phase field.The first one is  0 (x) = tanh( 0 ) that defines a symmetric and equal distributed phase field where the phases are separated.The second one defines a random initial condition:  1 = min{max{ φ1 , −1}, 1}, where φ is given by Gaussian random field: where {x  }  0 ⊂ S and  =  = 100.This allows 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  but set H 0 =  = 0. Other parameters are set to  = 0.02,  = 0.001 and , which provides a good compromise between computational effort and physical accuracy.Numerical parameters are such that ℎ 3 ∼ .

Convergence study
Due to the tight coupling between the phase field , the surface velocity , the surface pressure , and the geometry 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  2 norm in space and the  ∞ norm in time.The norms are calculated with respect to the linear reference geometry S  ℎ , where the evaluation of the norms on different surfaces S, S ℎ , S  ℎ leads to equivalent norms with the same order of convergence [18].We define the following errors: with discrete solutions  ℎ ,  ℎ ,  ℎ with varying ℎ and , ,  the discrete reference solutions with the finest mesh size.Additionally, we measure the area conservation error   , the local inextensibility error  div  , and the error of the total energy  F , where again F  and F denote the discrete reference solution on the finest mesh size.The simulations are done for the random initial condition  1 with Re = 1.0 and  1 =  2 =  = 0.02.
The results are shown in Figure 1.They indicate third order convergence for the surface error   and the error of the phase field   , which corresponds to the results in [65,18] for simple mean curvature flow and scalar-valued surface PDEs.The order of convergence for the phase field   also agrees with the one obtained for the corresponding problem on stationary surfaces considered in [7].For the error of the velocity field   the results indicate second order convergence in accordance with the results in [44] for one-component fluid deformable surfaces.In [14], third order convergence of   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 [34,35], 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   and the inextensibility error  div  we get second order convergence with respect to ℎ.Both are related to each other.While   = 0, the same order for  div  has been obtained in [7] 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 [20].However, overall we see experimentally the expected optimal order of convergence for the full problem.

Phase separation, bulging and induced flow
To demonstrate the strong interplay between composition, curvature, and hydrodynamics we consider  1 as initial condition.The other parameters are Re = 1.0 and either  1 =  2 =  = 0.02 or  1 = 0.02 and  2 = 0.5.In the last case, the red colored phase has a lower bending stiffness and is therefore expected to be guided to or initiate regions with higher mean curvature, while the blue color 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) and 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 shape changes, 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 a) and b) are also visible.While in a) the red and blue phases evolve similarly, in 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 behavior of the simulation shown in Figure 2 are quantified in Figure 3 by different energies.In the current setting, the evolution is dominated by F  , which is reduced by coarsening but also by bulging, which is associated with an increase in F  .Hydrodynamics seems to play a minor role in the evolution.F  is 1-2 orders of magnitude lower than F  .The strong increase of F  and F  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., H and H , is shown in Figure 3 c).The deviation between H 1 and 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.

Variation of parameters
We now explore the influence of the parameters in more details.To reduce the number of parameters we only consider the situation of a constant bending stiffness  1 =  2 = , set  = {0.02,0.1, 0.5}, and vary the Reynolds number with Re = {0.1, 1, 3}.We consider both initial conditions  0 and  1 .The behavior is demonstrated for the symmetric initial phase field  0 in Figure 4 and for the random initial phase field  1 in Figure 5.We only show the final configuration, this is either the equilibrium configuration or the configuration at the critical time  before a potential pinch-off.A  The corresponding critical times  are shown in Table 1.
table with critical final times is presented in Table 1.In the first case ( 0 ) the final configuration is mainly determined by the interplay of the Helfrich energy F  and the Ginsburg-Landau energy F  .For  = 0.5, F  dominates.Any deformation of the initial sphere requires to increase the Helfrich energy.Only small deformations are possible.The interface length, and thus F  , 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 F  .This is only realized for settings with lower  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  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,  = {0.1,0.02}, the critical time is drastically reduced if the Reynolds number Re is increased.
For the random initial condition  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  = 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 .For  = 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 less islands present.While it is known in flat space that hydrodynamics can enhance coarsening and this is also demonstrated on stationary surfaces [7], 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 .For  = 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 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.

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 [44] and two-phase flow models on stationary surfaces [7], 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 , the Reynolds number Re, and the bending stiffness  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 [53]; in the context of bulges with pinch-offs, this can be associated with endocytosis and exocytosis [43,3].While quantitative comparisons with experimental data in these context require further studies, the mathematical model and the numerical algorithm to solve it with the required accuracy are provided.
Funding This work was supported by the German Research Foundation (DFG) within the Research Unit "Vector-and Tensor-Valued Surface PDEs" (FOR3013).

Declaration of interests
The authors report no conflict of interest.
where Γ • •• are the Christoffel symbols determined by the metric tensor.The relation to the tangential and componentwise derivatives is given by: By definition, the covariant and tangential derivatives of vector fields are purely tangential operators but the componentwise derivative contains normal components too.We define the corresponding divergence operators by: div

Tensor fields
Let  ∈  2 R 3 | S be a tensor field with arbitrary smooth extension   , i. e.   | S = .The definition of the componentwise, tangential, or covariant derivative is possible for arbitrary tensor fields but we will show only the definitions of the operators used in this paper.To do so, we define the componentwise derivative by: The covariant divergence for the pure tangential tensor field  =  ∈  2 S is defined by: We define a componentwise divergence for right side tangential tensor field  by: = div S () − B + ((, B) + div S ())  .
The reason for defining only this restricted version of the componentwise divergence is that tr • ∇  would not be the adjoint operator of ∇  and thus the integration by parts formula (31) would no longer hold in general.

Laplace operators
For scalar fields, we define a Laplace operator by the Laplace Beltrami operator Δ S : and the componentwise Laplace operator for vector fields  =     by : where   are the Cartesian components with respect to the Cartesian basis vectors   .

B Lagrange-D'Alembert principle for the full model B.1 General proceeding
For state variables  (parametrization of S) and  ∈  0 S (phase field), and process variables  ∈ R 3 | S (material velocity) and  ∈  0 S (phase field rate), the Lagrange-D'Alembert principle generally reads (see [50]): where D is the dissipation potential and In a first step, we localize (35) in time.Assuming that the potential energy F depends only on state variables, we only have to consider the kinetic energy F  for this purpose.Note that time integration commutes with spatial variations in a non-relativistic setting.For instance, the following holds ,  .In contrast, spatial integration does clearly not commutes with spatial variation.From [58], we adopt the identity: where ð  :  0 S →  0 S is the deformation derivative in direction of , i. e. , We assume a variational mass conservation, i. e. , valid for mass density  ∈  0 S. Since the subset  ⊆ S is arbitrary, ð   = − div   holds locally.In addition, we assume temporal mass conservation,  = − div  , and vanishing variation directions at times  0 and  1 , i. e.  |  0 =  |  1 = 0.The fundamental theorem of calculus, temporal integration by parts, transport formula, and ð   =  Cartesian-componentwise together yields: Finally, the temporal localized version of the Lagrange-D'Alembert principle (35) reads: Note that, in our derivations, we assume a priori independence between the phase field  and the surface given by the parametrization .Therefore, ð   = 0 is valid for all  ∈ R 3 | S , see [58] for more details, where this a priori condition is referred to the scalar gauge of surface independence.In contrast to  2 -gradient flow techniques, here this assumption is not necessary to determine (37), since all terms containing ð   would vanish anyway.However, taken the generalized applied forces individually would comprises terms of the undetermined quantity ð   in a weak sense.Consequently, we would have to choose ð   = 0 to give this forces a determined meaning.

B.2 Lagrangian forces
In this section we consider the potential energy F = F  + F  , which is composed by the Ginzburg-Landau energy and the Helfrich energy 2).Without much mathematical effort, computing the negative generalized applied forces D  F  , D  F  ∈  0 S gives: To derive the Ginzburg-Landau force −D  F  , we use that ð     = −( [∇  ]   + [∇  ]  ) is valid for    = (  ,   ) [58].Moreover, the relation ð     = 0 holds, according to the assumption ð   = 0.This results in the following: by symmetry of the outer product.With (36) and (25), we obtain: Eventually, integration by parts (30) yields: In terms of stress, this means that the Ginzburg-Landau energy induces a trace-free and symmetric tangential stress for phase separations and a volumetric stress for the double-well potential.To separate the induced tangential and normal forces, we use metric compatibility ∇ S  = 0 and surface divergence (29).This results in: To calculate the Helfrich force −D  F  we need the deformation derivative of the mean curvature.From [58], we already know the tangential part of the deformation derivative of the shape operator.Explicitly, (ð  B) = ∇ S ( ∇  ) − B ∇  , where ð  B is to be read Cartesiancomponentwise.Since tr and ð  commute in this way, we obtain ð  H = div S ( ∇  ) − (B, ∇  ).With (36) and integration by parts for the covariant divergence, we obtain: and then: by integration by parts (31).In terms of stress, we get an orthogonal decomposition of a trace-free and symmetric tangential stress tensor, a volumetric stress tensor scaled by H 0 (), and an additional stress tensor operating in the normal-tangential space  ⊗ S.Using (29), metric compatibility ∇ S  = 0, the fact that B is curl-free and thus div S B = ∇ S H , and ∇ S  () =  ′ () ∇ S  for  ∈ {, H 0 }, the results in the tangential and normal Helfrich forces are:

B.3 Dissipative forces
In Variation with respect to  gives  D   ,  =  ∇   + (∇  )  , ∇   by symmetry.Therefore, by applying integration by parts (30) we obtain the negative viscous force: containing twice the tangential strain rate tensor.The friction potential does not depend on .Therefore, we have: We approach the immobility potential by defining a scalar field [  ] ∈  0 S implicitly so that it solve the equation Δ S [  ] =  ∈  0 S. Due to this, we can write the immobility potential in a  2 -manner by since this is the only mechanism that allows energy exchange with the surrounding of S, and the total energy F  + F must not increase with time.Using ( 52) -( 54), including  = 0, and applying the chain rule yields: Integration by parts ( 32) and ( 31), D  D  (45), and D  D  (47) results in: As a consequence, we finally obtain: which satisfies our requirements above.

B.6 Material time derivative
In this section we explain the material time derivative in details.Following the description in [59] for scalar fields  ∈  0 S, the material derivative can be written as: where  ∈  1 S is the so-called relative velocity.To be able to evaluate both summands on the right-hand-side, we use the two different parametrizations  : (,  1 , the path of a single material particle in time.In contrast,  is arbitrary as long as it depicts the same moving surface.This means that  acts as an observer of S. With this clarified, the scalar field  can be represented by evaluating (,  1 ,  2 ) as well as φ(, ()), where (()), respectively.To make it clear,  and φ are the very same scalar field from a physical point of view, only the mathematical representation differs.Partial derivatives depend on the mathematical description of their argument, contrarily to total derivatives, for example.Hence, we have to be careful when evaluating the first summand in (56).This partial time derivative reads: or in terms of a differential quotient: with respect to global coordinates () ∈ S(), where    is equal to    due to the relation () =  (,  1 , 2 ).As a consequence,    is well defined and can be discretize in time with respect to global coordinates.Since    represents only the observer rate of , the second summand in (56) represents the correction to obtain the material rate.The relative velocity is given by  =  − , where  is the material velocity and  the observer velocity.In local observer coordinates, they are given by: As a consequence, we have: , where the relation between global and material local coordinates is given by , respectively.Note that for geometrical reasons, (, ) = (, ) holds.Due to this, the relative velocity is a tangential vector field, i. e.  =  −  ∈  1 S, and therefore ∇   = (∇ S , ) is valid.Putting all together, we can write the material derivative ( 56 in terms of global coordinates and differential quotients.This representation is very suitable to implement flow equations on a surface grid that does not necessarily have to follow the material flow.Here, () is given by the surface grid coordinates and the material velocity ũ is part of the solution defined on these grid coordinates.To simplify the notation in all the other sections, we do not make use of the tilde, since all the field quantities with or without the tilde are exactly the same fields.They differ only within the argument formulation, which is not used in any case outside of this section.Note that we can also apply (56), or (57) respectively, for the material acceleration  with respect to a Cartesian frame   .Since   = 0, the relation  =     holds, for  =     .
There are other representations of the material derivative (56).For instance [21] when re-written in our notation2, where φ(,   ()) = (,  1 ,  2 ) for ( 1 ,  2 ) =  | −1  (  ()) and û(,   ()) =     ().As already mentioned in [21], both summands on the right-hand-side are not defined on a geometrically nonstationary surface, if considered individually.For the partial time derivative we have to evaluate φ( + ,   ()) for a small time step .However, φ| + is only defined at future locations   ( + ) ∈ S( + ) and not at the current S() if no further assumptions are considered.The same applies to the second summand, which depends on the normal derivative (∇ φ, ) if the material velocity û comprises a normal velocity.Such a normal derivative cannot be obtained without considering its extensions outside of the surface.One might be tempted to interpret (58) as a fully Eulerian perspective, since (  φ)|   ( )= ( , 1 , 2 ) =    is valid for    = 0.However, a fully Eulerian surface observer does not exists if the surface comprises a flow in normal direction.Any normal part of the surface motion has to be treated in a Lagrangian perspective at least.In general, ( 58) is equal to (56), since (  φ)|   ( )= ( , 1 , 2 ) =    − (∇, ) holds formally in arbitrary observer coordinates.

C Non-dimensionalization
In order to non-dimensionalize the model in Problem 1, we write the rescaled variables by using the symbol •.We define the rescaled parametrization  = X , the velocity  = û, the mean curvature H = Ĥ /, the phase field  = φ, and the chemical potential  =  μ 2 for the characteristic length  and a characteristic velocity , which defines a rescaled time  = / t.According to the parametrization, we rescale the gradient of the embedded space by ∇ = ∇/.This rescales the surface operators accordingly, e.g., the material time derivative reads (  + ∇  ) = / ( t + ∇ ŵ ).We introduce the rescaled material functions: κ() = ()/( 2  2 ), Ĥ0 () = H 0 (), and Ŵ () = 1/ 2  (); and the rescaled material parameters: Re = /, denoting the Reynolds number, m = /, σ = σ/( 2  2 ), and γ = /().Considering the rescaled variables and replacing the material functions and parameters into the model, we obtain the following non-dimensional system: where b and b correspond to the rescaled tangential and normal bending terms, respectively.
In the main sections of the paper, we drop the symbol • for a better readability.

D Derivation of the overdamped limit model
We consider the overdamped limit of the model in Problem 2. We define the rescaled parameters  = 1/, σ =  σ, κ = , and m = /, and the rescaled variables p =   and μ = .
In [36], the material derivative  is reduced to the partial time derivative   , and this leads to the differences noted between the models.More commonly used models reduce the local inextensibility constraint and only consider a global area constraint: with  =   +    −   .This system has also been considered in [36] and shows the same conceptional differences discussed above.Other approaches, with H 0 = 0, κ() = κ, and that use a penalization approach to enforce the global area constraint instead of the Lagrange multiplier, have been considered in [23].

1 Figure 1 :
Figure1: Convergence study for two-phase fluid deformable surfaces with respect to the mesh size ℎ 3 ∼  for  ℎ ,  ℎ ,  ℎ (top row) and error for area , divergence  P u ℎ and total energy F  + F (bottom row).

Figure 2 :
Figure 2: (a),(b) Snapshots of the relaxation of the two-component fluid deformable surface with random initial condition  1 and Re = 1.0 for  = 0, 0.3, 0.8, 1.1 (from left to right), with constant bending stiffness  1 =  2 =  = 0.02 in (a) and phase-depended bending stiffness  1 = 0.02 and  2 = 0.5 in (b).In (b) the red colored phase is less stiff than the blue colored phase.Top row: phase field .Bottom row: tangential velocity .The flow is visualized by a LIC filter and color coding represents the magnitude.Corresponding movies are provided in the supplementary data.

1 Figure 3 :
Figure 3: (a),(b) The energies F  , F  , F  and F  + F over time where (a) corresponds to fig 2 (a) and (b) corresponds to fig 2 (b).The time instances are highlighted in the plots.(c) Averaged mean curvature H 1,2 for the different phases with respect to the simulations done in fig 2 (a) and (b), the color corresponds to the colored phases.

Figure 4 : 8 Table 1 :
Figure 4: Final configuration obtained for different parameters and initial condition  0 .Either the simulation reaches the equilibrium configuration or the criteria for a potential pinch-off is reached.The corresponding critical times  are shown in Table1.

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

𝜌 2 ∥𝒖∥ 2 ,
in which the Lagrangian is defined by L = F  − F , with kinetic energy F  = ∫ S potential energy F , and time interval [ 0 ,  1 ].
where the negative generalized applied forces D  F , D  D ∈ R 3 | S and D  F , D  D ∈  0 S are given as  2 -gradients, i. e. (D  F , ) = (  F   , ) and (D  D, ) = (  D  , ) for all virtual displacements  ∈ R 3 | S , as well as (D  F , ) = (  F   , ) and (D  D, ) = (  D   , ) for all virtual displacements  ∈  0 S. Next, in Section B.2, we calculate the negative Lagrangian forces D  F and D  F , and the negative dissipative forces D  D and D  D in Section B.3, according to the potential energy F and the dissipation potential D given in this paper.Moreover, we implement local inextensibility of the material by the Lagrange-multiplier technique in Section B.4.Eventually, the strong formulation of (37) and a tangential-normal splitting of the Lagrangian force to   := −D  F and   := −(, D  F ) leads to the full model given in Problem 1.