Numerical modelling of phase separation on dynamic surfaces

TL;DR

This paper develops a fully Eulerian numerical method for simulating phase separation on dynamic, deforming surfaces using a surface-independent mesh, enabling analysis of complex surface phenomena like splitting and pattern formation.

Contribution

It introduces a hybrid finite difference and trace finite element method for Cahn-Hilliard equations on evolving surfaces, avoiding surface triangulation and handling topological changes.

Findings

Method accurately models phase separation on moving surfaces.

Capable of simulating topological transitions like surface splitting.

Successfully applied to pattern formation and droplet splitting scenarios.

Abstract

The paper presents a model of lateral phase separation in a two component material surface. The resulting fourth order nonlinear PDE can be seen as a Cahn-Hilliard equation posed on a time-dependent surface. Only elementary tangential calculus and the embedding of the surface in $\mathbb{R}^3$ are used to formulate the model, thereby facilitating the development of a fully Eulerian discretization method to solve the problem numerically. A hybrid method, finite difference in time and trace finite element in space, is introduced and stability of its semi-discrete version is proved. The method avoids any triangulation of the surface and uses a surface-independent background mesh to discretize the equation. Thus, the method is capable of solving the Cahn-Hilliard equation numerically on implicitly defined surfaces and surfaces undergoing strong deformations and topological transitions. We assess the approach on a set of test problems and apply it to model spinodal decomposition and pattern formation on colliding surfaces. Finally, we consider the phase separation on a sphere splitting into two droplets.

Equations

$$\mathbf{u}={\mathbf{u}}_T+u_N\mathbf{n},{\mathbf{u}}_T\cdot \mathbf{n}=0,$$

$$\Psi (0,\mathbf{y})=\mathbf{y},\frac{\partial \Psi (t,\mathbf{y})}{\partial t}=\mathbf{u}(t,\Psi (t,\mathbf{y})),t\in [0,T].$$

$${\mathrm{div}}_{\Gamma }\mathbf{f}:=\mathrm{tr}({\nabla }_{\Gamma }\mathbf{f}),{\Delta }_{\Gamma }f:={\mathrm{div}}_{\Gamma }({\nabla }_{\Gamma }f),$$

$${\int }_{S(t)}f\text{div}{}_{\Gamma }\mathbf{v}ds={\int }_{\partial S(t)}f\mathbf{v}\cdot \mathbf{m}d\gamma -{\int }_{S(t)}\mathbf{v}\cdot {\nabla }_{\Gamma }fds+{\int }_{S(t)}\kappa f\mathbf{v}\cdot \mathbf{n}ds,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{S(t)}fds={\int }_{S(t)}\left(\dot{f}+f{\mathrm{div}}_{\Gamma }\mathbf{u}\right)ds,$$

$$\dot{f}=\frac{\partial f}{\partial t}+\mathbf{u}\cdot \nabla f.$$

$$\mathcal{S}=\bigcup _{t\in [0,T]}\Gamma (t)\times \{t\}.$$

$$0=\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{S(t)}\rho ds={\int }_{S(t)}\left(\dot{\rho }+\rho {\mathrm{div}}_{\Gamma }\mathbf{u}\right)ds.$$

$$\dot{\rho }+\rho {\mathrm{div}}_{\Gamma }\mathbf{u}=0\text{on}\Gamma (t).$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{S(t)}\rho cds=-{\int }_{\partial S(t)}\mathbf{j}\cdot \mathbf{m}d\gamma .$$

$${\int }_{S(t)}(\dot{\rho }c+\dot{c}\rho +c\rho {\mathrm{div}}_{\Gamma }\mathbf{u})ds=-{\int }_{S(t)}{\mathrm{div}}_{\Gamma }\mathbf{j}ds$$

$${\int }_{S(t)}\dot{c}\rho ds=-{\int }_{S(t)}{\mathrm{div}}_{\Gamma }\mathbf{j}ds$$

$$\dot{c}+\frac{1}{\rho }{\mathrm{div}}_{\Gamma }\mathbf{j}=0\text{on}\Gamma (t),$$

$$\mathbf{j}=-M{\nabla }_{\Gamma }\mu \text{on}\Gamma ,\mu =\frac{\delta f}{\delta c},$$

$$f(c)=\frac{1}{ϵ}{f}_0(c)+\frac{ϵ}{2}|{\nabla }_{\Gamma }c{|}^2.$$

$$\dot{\rho }+\rho {\mathrm{div}}_{\Gamma }\mathbf{u}$$

$$\dot{c}-{\rho }^{-1}{\mathrm{div}}_{\Gamma }\left(M{\nabla }_{\Gamma }\left(\frac{1}{ϵ}{f}_0^{\prime }-ϵ{\Delta }_{\Gamma }c\right)\right)$$

$$\rho (\cdot ,t)={\rho }_0,\text{and}c(\cdot ,t)=c_0\text{on}\Gamma (0).$$

$$\dot{\rho }=-\kappa u_N\rho .$$

$$\dot{c}-{\rho }^{-1}{\mathrm{div}}_{\Gamma }\left(M{\nabla }_{\Gamma }\mu \right)=0\text{on}\Gamma (t),$$

$$\mu =\frac{1}{ϵ}{f}_0^{\prime }-ϵ{\Delta }_{\Gamma }c\text{on}\Gamma (t).$$

$$M=M(c)=\sigma c(1-c),$$

$$f_0(c)=\frac{1}{4}{c}^2(1-c{)}^2.$$

$${\int }_{\Gamma (t)}\rho \dot{c}vds+{\int }_{\Gamma (t)}M{\nabla }_{\Gamma }\mu {\nabla }_{\Gamma }vds=0,$$

$${\int }_{\Gamma (t)}\mu qds-\frac{1}{ϵ}{\int }_{\Gamma (t)}{f}_0^{\prime }(c)qds-ϵ{\int }_{\Gamma (t)}{\nabla }_{\Gamma }c{\nabla }_{\Gamma }qds=0.$$

$$H^{1,\Gamma }$$

$$L_{1,\Gamma }^{\infty }$$

$$\dot{c}=\frac{\partial c}{\partial t}+\mathbf{u}\cdot \nabla c,$$

$${\mathcal{O}}_{\delta }(\Gamma (t))=\{\mathbf{x}\in {\mathbb{R}}^3:\text{mbox}dist(\mathbf{x},\Gamma (t))<\delta \}.$$

$$\frac{\partial c^e}{\partial {\mathbf{n}}^e}=\frac{\partial {\mu }^e}{\partial {\mathbf{n}}^e}=0\text{in}\mathcal{O}(\mathcal{S}).$$

$${\int }_{\Gamma (t)}\rho (\frac{\partial c}{\partial t}+\mathbf{u}\cdot \nabla c)vds+{\int }_{\Gamma (t)}M{\nabla }_{\Gamma }\mu {\nabla }_{\Gamma }vds$$

$${\int }_{\Gamma (t)}\mu qds-\frac{1}{ϵ}{\int }_{\Gamma (t)}{f}_0^{\prime }(c)qds-ϵ{\int }_{\Gamma (t)}{\nabla }_{\Gamma }c{\nabla }_{\Gamma }qds$$

$$\frac{\partial \mu }{\partial \mathbf{n}}=\frac{\partial c}{\partial \mathbf{n}}$$

Full text

1 Introduction

The Cahn–Hilliard (CH) equation on a stationary, flat domain has been introduced in the late 50s to model phase separation in binary alloy systems [8, 9]. As a prototype model for segregation of two components in a mixture, over the years it has been used in many areas beyond materials science. Applications in the biology are particularly numerous. For example, CH type equations have been used to model and simulate tumor growth [68, 29], dynamics of plasma membranes and multicomponent vesicles [5, 40, 59, 36, 43, 3, 4], and lipid rafts formation [22]. In some applications, such as sorting in biological membranes, phase separation and coarsening happen in a thin, evolving layer of self-organizing molecules, which in continuum-based approach can be modeled as a material surface. This motivates the recent interest in the CH equation posed on time-dependent surfaces.

The Cahn–Hilliard equation is challenging to solve numerically due to non-linearity, stiffness, and the presence of a fourth order derivative in space. For some recent publications on the CH equation in planar and volumetric domains, we refer to [28, 64, 38, 10] and references therein. The numerical solution of the CH equation posed on surfaces is further complicated by the need to discretize tangential differential operators and to approximately recover complex shapes. Several authors have opted for a finite difference method. For example, the closest point finite difference method was applied to solve the CH equation on a stationary torus in [23] and on more general stationary domains in [31]. A finite difference method for a diffuse volumetric representation of the surface CH equation was introduced in [27]. However, a finite element method (FEM) is often considered to be the most flexible numerical approach to handle complex geometries. Concerning the CH equation posed on a stationary surface, the convergence of a FEM was studied in [14], where numerical examples for a sphere and saddle surface are provided; results obtained by FEM on more general surfaces are presented in [22, 37]. In [46], solutions to the surface Cahn–Hilliard–Navier–Stokes equation were computed on a sphere and torus. All of the above references use a sharp surface representation and a discretization mesh fitted to the surface. In [69], we studied for the first time a geometrically unfitted finite element method for the surface Allen–Cahn and Chan–Hilliard equations. In the unfitted FEM, the surface is not triangulated in the common sense and may overlap a background computational (bulk) mesh in an arbitrary way. Moreover, stability and accuracy of the discretization do not depend on the position of the surface relative to the bulk mesh. Thus, the next natural development is to allow the surface to evolve through the time-independent bulk mesh that is used to discretize the equation posed on the dynamic surface itself. Such development is presented in this paper, which extends and studies the unfitted FEM for CH equations on time-dependent surfaces.

Very few works deal with a CH equation posed on an evolving surface. In [19], the authors show a rigorous well-posedness result for a continuous CH-type equation, which is a simplification of the model for surface dissolution set out in [17]. The FEM used for the space discretization in [19], known as evolving surface finite element method, relies on evolving an initial surface triangulation by moving the nodes according to a prescribed velocity. The asymptotic limit (as the interfacial width parameter tends to zero) of the CH equation on a surface evolving with prescribed velocity is studied theoretically and numerically in [47]. The authors of [4] studied a finite element method for the bulk Navier–Stokes equations coupled to the surface CH model, where the surface evolution is driven by the bulk fluid dynamics and a curvature energy. Again, in [17, 19, 47, 4], the discretization mesh is fitted to the computational surface and evolves with it. Although this approach offers many advantages, problems arise when strong deformations and topological changes of the surface occur. The numerical approximation of the surface CH equation using an isogeometric approach was studied in [70]. The point of using isogeometric analysis is that spline bases with high-order and high-continuity allow to treat the fourth order problem, without resorting to a mixed formulation. In [70], an isogeometric finite element formulation was introduced for the CH equation written in intrinsic surface variables on a surface evolved by the PDE derived according to Kirchhoff–Love shell theory. Finally, we mention a different approach to study surface phase distribution [67, 40, 65]: a pair of phase-field variables is introduced such that one variable characterizes the surface (in a diffuse manner) while the other describes the distribution of the surface phases.

In this paper, we use elementary tangential calculus to derive a CH equation posed on an evolving material surface. There is a (non-essential) difference between the equation we get compared to the one in [19, 47], which is explained in Remark 2.1. The more substantial difference with the previous works is that we study a geometrically unfitted finite element method. This method builds upon earlier work on an unfitted FEM, known as trace finite element method (TraceFEM), for elliptic PDEs on stationary surfaces in [50] and evolving surfaces in [35]. Unlike some other geometrically unfitted methods for surface PDEs, TraceFEM employs a sharp surface representation. The surface can be defined implicitly, e.g. as the zero of a level-set function, and no knowledge of the surface parametrization is required. We shall see that the developed method is very well suited for CH equation posed on evolving surfaces.

The work presented in here is a first step in the direction of understanding the evolution of multicomponent lipid vesicles used as drug carriers and their fusion with the target cell. Membrane fusion is recognized as a potentially efficient mechanism for the delivery of macromolecular therapeutics to the cellular cytoplasm. The process of lipid membrane phase separation to concentrate binding lipids within distinct regions of the membrane surface has been shown to enhance membrane fusion [30]. In addition, recent developments on targeted lipid vesicles suggest that the formation of reversible phase-separated patterns on the vesicle surface increase target selectivity, cell uptake and overall efficacy [1, 32]. The numerical experiments presented in Sec. 4.2.2 mimic phase separation occurring on a two component lipid vesicle, leading to fusion with another vesicle. These and other numerical results in Sec. 4 showcase the ease with which the numerical method handles changes in the surface topology, making it a perfect computational tool to support and complement experimental practice in the design of drug carriers.

The outline of the paper is as follows. In Section 2, we will derive from conservation laws a Cahn–Hilliard equation on an evolving surface. The hybrid finite difference in time – finite element in space variant of the TraceFEM for the CH equation is introduced in Section 3. After the assessment of the numerical method for a set of model problems, in Section 4 we will study phase separation modeled by our surface Cahn–Hilliard equation on evolving surfaces with topological transitions: two colliding spheres and a sphere splitting into two droplets. Section 5 closes the paper with a few concluding remarks.

2 Mathematical model

2.1 Preliminaries

Let $\Gamma(t)$ be a closed, smooth, evolving surface in $\mathbb{R}^{3}$ for $t\in[0,T]$, where $T>0$ is a final time. Let the 3D bulk domain $\Omega$ be such that $\Gamma(t)\subset\Omega$, for all $t\in[0,T]$. We assume surface $\Gamma(t)$ evolves according to a given smooth velocity field $\mathbf{u}$, i.e., $\mathbf{u}(\mathbf{x},t)$ is the velocity of point $\mathbf{x}\in\Gamma(t)$. Consider the decomposition of $\mathbf{u}$ into normal $(\mathbf{u}_{N}=u_{N}\mathbf{n})$ and tangential $(\mathbf{u}_{T})$ components:

[TABLE]

where $\mathbf{n}$ is the outward normal vector on $\Gamma$. While the normal velocity $u_{N}$ completely defines the geometric evolution of the surface, we are interested in the motion of $\Gamma(t)$ as a material surface, which also includes tangential deformations. The material deformation can be defined through the Lagrangian mapping $\Psi(t,\cdot)$ from $\Gamma(0)$ to $\Gamma(t)$, i.e. for $\mathbf{y}\in\Gamma(0)$, $\Psi(t,\mathbf{y})$ solves the ODE system

[TABLE]

For any sufficiently smooth function $f$ in a neighborhood of $\Gamma(t)$ its tangential gradient is defined as $\nabla_{\Gamma}f=\nabla f-(\mathbf{n}\cdot\nabla f)\mathbf{n}$. The tangential (surface) gradient $\nabla_{\Gamma}f$ then depends only on values of $f$ restricted to $\Gamma(t)$ and $\mathbf{n}\cdot\nabla_{\Gamma}f=0$ holds. For a vector field $\mathbf{f}$ on $\Gamma(t)$ we define $\nabla_{\Gamma}\mathbf{f}$ componentwise. The surface divergence operator for $\mathbf{f}$ and the Laplace–Beltrami operator for $f$ are given by:

[TABLE]

where ${\rm tr}(\cdot)$ is the trace of a matrix. For a fixed $t\in[0,T]$, $L^{2}(\Gamma(t))$ is the Lebesgue space of square-integrable functions on $\Gamma(t)$ and $H^{1}(\Gamma(t))$ is the Sobolev space of all functions $f\in L^{2}(\Gamma(t))$ such that $\nabla_{\Gamma}f\in L^{2}(\Gamma(t))^{3}$. For a subdomain $S(t)\subset\Gamma(t)$, we recall the integration by parts identity:

[TABLE]

where $\kappa$ is the sum of principle curvatures and $\mathbf{m}$ is the normal vector on $\partial S(t)$ that is tangential to $\Gamma(t)$ and outward for $S(t)$; vector field $\mathbf{v}$ and scalar function $f$ are such that all the quantities in (2) exist.

Let $S(t)\subset\Gamma(t)$ be a material area evolving according to (1). Then, the following surface analogue of the Reynolds transport theorem holds (see, e.g., Lemma 2.1 in [16]):

[TABLE]

where $\dot{f}$ is the material derivative of $f$. One can write the material derivative in terms of partial derivatives

[TABLE]

The terms on the right-hand side of (4) are well defined if one identifies $f$ with its arbitrary smooth extension from the surface to a neighborhood of the space-time manifold $\mathcal{S}\subset\mathbb{R}^{4}$,

[TABLE]

Of course, $\dot{f}$ is an intrinsic surface quantity and the quantity on right-hand side of (4) is independent of the choice of smooth extension.

Now, we are prepared to set up the mathematical model to describe the separation of two conserved phases on $\Gamma$.

2.2 Model setup

On $\Gamma(t)$ we consider a heterogeneous mixture of two species with densities $\rho_{i}$, $i=1,2$. Let $\rho=\rho_{1}+\rho_{2}$ be the total density. The conservation of mass written for an arbitrary material area $S(t)\subset\Gamma(t)$ and (3) yield

[TABLE]

Since the above identity holds for arbitrary $S(t)\subset\Gamma(t)$, it implies

[TABLE]

Following [8, 9, 41], to describe the dynamics of phases we introduce specific mass concentrations $c_{i}=m_{i}/m$, $i=1,2$, where $m_{i}$ are the masses of the components and $m$ is the total mass. Since $m=m_{1}+m_{2}$, it holds $c_{1}+c_{2}=1$. Let $c_{1}$ be the representative concentration $c$ (order parameter), that is $c=c_{1}$ and $c\in[0,1]$. Mass conservation for one component on $S(t)\subset\Gamma(t)$ takes the form

[TABLE]

where $\mathbf{j}$ is a mass flux.

By applying the transport formula (3) to the left-hand side of (6) and the integration by parts formula (2) to the right-hand side of (6), we obtain

[TABLE]

Thanks to (5), the above equality simplifies to

[TABLE]

Since $S(t)$ can be taken arbitrary, we get

[TABLE]

The classical assumption for the flux $\mathbf{j}$ is Fick’s law:

[TABLE]

where $M$ is the so-called mobility coefficient (see [33]) and $\mu$ is the chemical potential, which is defined as the functional derivative of the total specific free energy $f$ with respect to the concentration $c$. The choice of $f$ defines a specific model of mixing. Following [8], we set

[TABLE]

In eq. (9), $f_{0}(c)$ is the free energy per unit surface, a non-convex function of $c$, while the second term represents the interfacial free energy based on the concentration gradient. In model (7)–(9), the interface between the two components is a layer of size ${\epsilon}$.

By combining eqs. (5), (7)–(9), we obtain the Cahn–Hilliard equation posed on evolving material surface $\Gamma(t)$:

[TABLE]

for all $t\in(0,T)$. We close system (10)–(11) with initial conditions

[TABLE]

There are no boundary conditions since the boundary of $\Gamma(t)$ is empty.

Note that eq. (10) is decoupled from eq. (11), in the sense that the evolution of the total density distribution depends only on the given surface motion, but not on the order parameter $c$. In the particular case of only normal velocities, i.e. $\mathbf{u}_{T}=0$, eq. (10) simplifies to

[TABLE]

For inextensible membranes, it holds ${\mathop{\,\rm div}}_{\Gamma}\mathbf{u}=0$ and eq. (10) further simplifies to $\dot{\rho}=0$.

Remark 2.1**.**

The system (10)–(11) is similar to the Cahn-Hilliard equation on evolving surface found in [19, 47], but different. In [19, 47], the equation is written for the conserved variable $u=\rho\,c$ and under the assumption that the free energy functional depends on $u$ rather than on the order parameter $c$, i.e. $f=f(u)$ and $\mu=\frac{\delta f}{\delta u}$. Formally, system (10)–(11) and the surface Cahn–Hilliard equation in [19, 47] are the same problem for inextensible membranes and if one assumes $\rho_{0}=const$. Here the situation partially resembles the coupling of a two-component compressible fluid flow with dissipative Ginzburg–Landau interface dynamics, where, depending on the choice of the variables to define the energy functional, the effects of compressibility has to be considered in the definition of the reactive stress tensor and the chemical potential; see, e.g., [41] for further discussion and references.

2.3 Weak form of the equations

Eq. (10) is a transport problem for the total density. It can be integrated independently of eq. (11) as long as the evolution of the surface is given. Since $\rho$ is completely determined by the initial data and $\mathbf{u}$, but not by the phase composition of $\Gamma(t)$, for the purposes of this paper we shall assume that $\rho$ is given. Eq. (11) is a fourth-order nonlinear evolutionary equation, which is more challenging numerically. In particular, casting it in a weak form would lead to second-order spatial derivatives. From the numerical point of view, it is beneficial to avoid higher order spatial derivatives. Hence, following the common practice we rewrite eq. (11) in mixed form, i.e. as two coupled second-order equations:

[TABLE]

System (13)–(14) needs to be supplemented with the definitions of mobility $M$ and free energy per unit surface $f_{0}$. A possible choice for $M$ is given by

[TABLE]

with $\sigma>0$. This mobility is referred to as a degenerate mobility, since it is not strictly positive. As for the free energy per unit surface, a common choice is given by the double-well potential

[TABLE]

We are interested in a finite element numerical method for the evolving surface Cahn–Hilliard problem (13)–(16). We first need a weak formulation of system (13)–(14). To devise it, we start with multiplying (13) by $\rho$ and a smooth $v:\mathcal{S}\to\mathbb{R}$, while we multiply (14) by a smooth $q:\mathcal{S}\to\mathbb{R}$. Next, we integrate over $\Gamma(t)$ and employ the integration by parts identity (2) with $S(t)=\Gamma(t)$ (this implies $\partial S(t)=\emptyset$). The integration by parts is applied to the diffusion terms in (13) and (14). The curvature terms vanish since the fluxes $M\nabla_{\Gamma}\mu$ and $\nabla_{\Gamma}c$ are tangential. We obtain the following equalities:

[TABLE]

A rigorous week formulation requires the definition of the test and trial functional spaces. Suitable spaces are formulated for functions defined on the space-time manifold $\mathcal{S}$: $L^{q}(\mathcal{S})$, $1\leq q\leq\infty$ and $H^{1}(\mathcal{S})$ are standard Lebesgue and Sobolev spaces on $\mathcal{S}$. We also need

[TABLE]

Note that $\nabla_{\mathcal{S}}$ and $\int_{\mathcal{S}}$ can be written in terms of tangential calculus on $\Gamma$ and the geometric motion of $\Gamma$, i.e. $u_{N}$ and $\mathbf{n}$. See [51, 19] for details and properties of the above spaces.

We introduce the week formulation of (13)–(14):

Find $c\in L^{\infty}_{1,\Gamma}\cap H^{1}(\mathcal{S})$ and $\mu\in H^{1,\Gamma}$ satisfying (17)–(18) for almost every $t\in(0,T)$ and for all $(v,q)\in H^{1,\Gamma}\times H^{1,\Gamma}$.*

For a closely related problem (which coincides with (17)–(18) if $\textrm{div}\ \!_{\Gamma}\mathbf{u}=0$, $\rho=1$, $M=\mbox{const}$), the above weak formulation was shown in [19] to have the unique solution. Extending the well-posedness analysis for (17)–(18) is out of the scope for this paper.

3 A hybrid finite difference / finite element numerical method

In this section, we present a completely Eulerian numerical method for solving problem (10)–(11). In our approach, the mesh does not follow the evolution of the surface, as is typical for Lagrangian methods. Examples of finite element methods for surface PDEs based on the Lagrangian description can be found, e.g., in [15, 2, 20, 19, 60, 42]. We rather allow the surface to travel through a background mesh without restrictions. Furthermore, here we look for a numerical scheme that extends, in some sense, the method of lines, which is (arguably) the most popular computational approach for parabolic problems in stationary domains. We recal that the method of lines provides the convenience of a separated numerical treatment of spatial and temporal variables. This is by no means straightforward for problem (10)–(11), since the equations are defined on a time-dependent surface and there is no evident way of separating variables in the Eulerian setting. To circumvent these difficulties, we build on ideas from [54, 35] and suggest a hybrid, finite-difference in time and finite element in space discretization of problem (10)–(11).

We start with several important observations about the differential problem (10)–(11) and its alternative integral formulation.

3.1 Extending $c$ and $\mu$ off the surface

We are able to decompose the material derivative into the sum of Eulerian terms,

[TABLE]

if we assume an arbitrary smooth extension of $c$ to a neighborhood of $\mathcal{S}$, denoted further by $\mathcal{O}(\mathcal{S})$. Let us define $\mathcal{O}(\mathcal{S})=\bigcup\limits_{t\in(0,T)}\mathcal{O}_{\delta}(\Gamma(t))\times\{t\}$, where

[TABLE]

We let $\delta$ be sufficiently small such that $\mathcal{O}_{\delta}(\Gamma(t))\subset\Omega$ for all $t\in(0,T)$.

Among infinitely many possible smooth extensions of $c$ to $\mathcal{O}(\mathcal{S})$, it will be convenient to assume the one given by the close point projection on $\Gamma(t)$: Fix $t\in(0,T)$, then for $\mathbf{x}\in\mathcal{O}_{\delta}(\Gamma(t))$ we denote its close point projection on $\Gamma(t)$ by $\mathbf{p}(\mathbf{x})$. For smoothly evolving $\Gamma(t)$, the projection $\mathbf{p}(\mathbf{x})$ is well defined in $\mathcal{O}(\mathcal{S})$, and we consider the extensions given by $c^{e}(\mathbf{x},t)=c(\mathbf{p}(\mathbf{x}),t)$ and $\mu^{e}(\mathbf{x},t)=\mu(\mathbf{p}(\mathbf{x}),t)$ for $(\mathbf{x},t)\in\mathcal{O}(\mathcal{S})$. We also introduce the extension of the spatial normal vector field, $\mathbf{n}^{e}(\mathbf{x},t)=\mathbf{n}(\mathbf{p}(\mathbf{x}),t)$. Recall that $\mathbf{n}(\cdot,t)$ is the outward normal for $\Gamma(t)$. Now $c^{e}$ and $\mu^{e}$ can be characterized as extensions along the directions given by the normal field, i.e. the so called normal extensions:

[TABLE]

If $\Gamma(t)$ is a $C^{2}$ surface, then its normal field is $C^{1}$-smooth and $c\in H^{1}(\Gamma(t))$ implies $c^{e}\in H^{1}(\mathcal{O}_{\delta}(\Gamma(t)))$. In turn, the regularity of $\Gamma(t)$ together with assumptions on its evolution (e.g., in terms of smoothness of the mapping $\Psi$) yield the space-time regularity of $c^{e}$ (if $c$ is smooth); see more details in, e.g., [48, Sec. 6.1].

We shall identify functions defined on $\mathcal{S}$ with their normal extensions and skip the upper index $e$. Once we do this, the weak formulation (17)–(18) yields the following identities for $c$ and $\mu$:

[TABLE]

for all smooth $v$ and $q$ defined in $\Gamma(t)$. We shall base our discretization method on equalities (21)–(23). The crucial idea is that our numerical solution should approximate the concentration and chemical potential on $\Gamma(t)$ together with their normal extensions in a suitable neighborhood of the discrete surface. This will allow for a completely Eulerian method with a separate treatment of spatial and temporal derivatives.

Next, we turn to describing the discrete formulation.

3.2 Time discretization

In the classical method of lines, one first discretizes spatial variables and leaves time continuous; next the resulting Cauchy problem is integrated numerically. Here, we follow the reverse order by first treating time derivatives, which is a popular alternative [63] and often leads at the end to the same fully discrete scheme. To this end, consider the uniform time step $\Delta t=T/N$, and let $t_{n}=n\Delta t$ and $I_{n}=[t_{n-1},t_{n})$. Denote by $c^{n}$, $\mu^{n}$ an approximation of $c(t_{n})$ and $\mu(t_{n})$, define $\Gamma_{n}:=\Gamma(t_{n})$, $n=0,\dots,N$. We assume that $\mathcal{O}(\mathcal{S})$ is a sufficiently large neighborhood of $\mathcal{S}$ such that

[TABLE]

see Fig. 1. In this case, $c^{n-1}$ is well-defined on $\Gamma_{n}$ and we can approximate the time derivative by simple finite difference

[TABLE]

This brings us to the implicit Euler method for problem (21)–(23):

[TABLE]

It is important that all terms in (25)–(26) are defined on the ‘current’ surface $\Gamma_{n}$ and its neighborhood. In particular, $c^{n-1}$ is well-defined on $\Gamma_{n}$ thanks to condition (24), which gives the meaning to $\left[{c}\right]_{t}^{n}$ in (25)–(26),

In practice, one can use a semi-implicit Euler method, which compromises the stability of the implicit method in favor of computational time savings:

[TABLE]

Following the idea of first order stabilization in [58], we included term $-\beta_{s}\Delta t\left[{c}\right]_{t}^{n}$ in eq. (26). On a stationary domain, for constant mobility and a slightly modified double-well potential $f_{0}$ than (16), the scheme above was shown in [58] to be stable for large enough $\beta_{s}=O(\epsilon^{-1})$. In the semi-implicit Euler method both nonlinear terms, i.e. the mobility coefficient in eq. (28) and the potential in eq. (29), are extrapolated from the previous time step. So, technically at each time step the method requires solving a linear system (28)–(29) and further extending the solution along normals, i.e. solving (30) with given data on $\Gamma_{n}$. We will show in Sec. 3.4 that both steps (solving for $\mu^{n}$ and $c^{n}$ on $\Gamma_{n}$ and extension) can be naturally combined in one linear solve on the finite element level.

3.3 Stability of the semi-discrete scheme

It is well known that in a steady domain the Cahn-Hilliard problem defines the $H^{-1}$-gradient flow of an energy functional. However, we are not aware of a minimization property for the Cahn-Hilliard problem in time-dependent domains. In Section 4, we shall illustrate by numerical examples that the evolution of $\Gamma(t)$ can actually produce a (local) increase of the system energy.

It is natural to say that (25)–(27) or (28)–(30) are numerically stable if the solution satisfy an energy bound uniform in the discretization parameter $\Delta t$. In this Section, we demonstrate such a bound for (25)–(27) subject to the following simplifications. First, we assume constant density and mobility: $\rho^{n}=1$, $M^{n}=1$. We also assume a tangential incompressibility condition:

[TABLE]

which corresponds to the motion of an inextensible membrane.

We start by taking $v=\Delta t\,\mu^{n}$ in (28) and $q=-\Delta t(\left[{c}\right]_{t}^{n}+\mathbf{u}^{n}\cdot\nabla c^{n})$ in (29) and sum the resulting equalities. After few cancellations, we obtain:

[TABLE]

where $\|\cdot\|_{\Gamma_{n}}$ stands for the $L^{2}(\Gamma_{n})$ norm. First, we treat the terms with ${f}_{0}^{\prime}$ on the right-hand side. Here and throughout the rest of this section, we use the fact that (27) yields $\nabla c^{n}=\nabla_{\Gamma}c^{n}$, $\nabla{f}_{0}(c^{n})=\nabla_{\Gamma}{f}_{0}(c^{n})$, etc. Using integration by parts, formula (2), and condition (31), we have:

[TABLE]

With the help of trunceted Taylor expansion and ${f}_{0}^{\prime\prime}\geq-\frac{1}{4}$, we get

[TABLE]

Next, we consider the terms on the left-hand side of (3.3). Using Cauchy–Schwarz inequality and estimate $ab\leq\frac{1}{2\delta}a^{2}+\frac{\delta}{2}b^{2}$, $\forall~{}a,b\in\mathbb{R}$ with $\delta\in(0,+\infty)$, we obtain:

[TABLE]

where constant $C$ depends on $\mathbf{u}(\mathbf{x},t)$. The last term on the left-hand side of (3.3) is handled by the integration by parts

[TABLE]

Combining (3.3)–(3.3), for $\beta_{s}>\frac{1}{2\epsilon}$ we get

[TABLE]

where the constant $C$ depends on $\mathbf{u}(\mathbf{x},t)$.

Before applying the discrete Gronwall inequality to pass from (37) to a priori stability estimate, we need to relate the quantities $\|\nabla_{\Gamma}c^{n-1}\|^{2}_{\Gamma_{n}}$ and $\int_{\Gamma_{n}}{f}_{0}(c^{n-1})\,ds$ to $\|\nabla_{\Gamma}c^{n-1}\|^{2}_{\Gamma_{n-1}}$ and $\int_{\Gamma_{n-1}}{f}_{0}(c^{n-1})\,ds$, respectively. To this end, consider the closest point projection $\mathbf{p}^{n-1}:\mathcal{O}_{\delta}(\Gamma_{n-1})\to\Gamma_{n-1}$. Since $c^{n-1}$ satisfies (30), we can write $c^{n-1}(\mathbf{x})=c^{n-1}(\mathbf{p}^{n-1}(\mathbf{x}))$ in $\mathcal{O}_{\delta}(\Gamma_{n-1})$. In particular, due to (24) this representation holds for $\mathbf{x}\in\Gamma_{n}$. For the surface measures on $\Gamma_{n-1}$ and $\Gamma_{n}$, we have [35, Lemma 1]:

[TABLE]

with some $C>0$ depending only on surface velocity $\mathbf{u}$. Therefore, for $\Delta t$ small enough we obtain

[TABLE]

In turn, the surface gradients on $\Gamma_{n}$ and $\Gamma_{n-1}$ are related through (cf., e.g., [12, section 2.3])

[TABLE]

where $\mathbf{P}_{\Gamma_{n}}$ is the orthogonal projector on the tangential space on $\Gamma_{n}$, $d(\mathbf{x})$ a signed distance function for $\Gamma_{n-1}$, and $\mathbf{H}(\mathbf{x})=\nabla^{2}d(\mathbf{x})$ is the shape operator. Our assumptions on $\Gamma(t)$ and its smooth evolution imply the uniform in $n$ bounds $\|\mathbf{H}\|_{L^{\infty}(\mathcal{O}_{\delta}(\Gamma_{n-1}))}\leq C$ and $\|d\|_{L^{\infty}(\Gamma_{n})}\leq C\Delta t$. Together with (38), this leads to estimate

[TABLE]

By employing (39) and (40) in (37), we are led to the following bound

[TABLE]

We sum the inequalities (3.3) over $n=1,\dots,N$ to get

[TABLE]

Now assume that $\Delta t$ is small enough such that $1-\Delta t\widetilde{C}(1+\frac{|\Delta t|}{\epsilon^{2}})=\alpha>0$ and apply the discrete Gronwall inequality to obtain

[TABLE]

where the constant on the right-hand side depends on given problem data, i.e., on $\mathbf{u}$, $\Gamma(0)$, $c^{0}$, but does not depend on $n$ and $\Delta t$. The constant $C({\rm data})$ is also uniformly bounded with respect to $\epsilon$ if $|\Delta t|\lesssim\epsilon$.

3.4 Fully discrete method

The (quasi)-stationary system (28)–(30) involves only integrals over $\Gamma^{n}$ and $\mathcal{O}_{\delta}(\Gamma^{n})$. This enables us to apply a surface finite element method developed for a steady surface. Below we consider an unfitted finite element method, known as the TraceFEM [49], to discretize problem (28)–(30) in space.

We start by assuming a shape regular triangulation $\mathcal{T}_{h}$ of the bulk computational domain $\Omega\subset\mathbb{R}^{3}$. At this point, we consider the same surface-independent mesh for all $t\in(0,T)$. However, in practice dynamic refinement is possible and the numerical examples in Section 4 demonstrate its feasibility and use. We will address this option later. The domain formed by all the elements in $\mathcal{T}_{h}^{\Gamma}$ is denoted by $\Omega^{\Gamma}_{h}$. On $\mathcal{T}_{h}$ we consider a standard finite element space of continuous functions that are piecewise-polynomials of degree $1$. This bulk (volumetric) finite element space is denoted by $V_{h}^{\rm bulk}$:

[TABLE]

Let the surface be given implicitly as the zero level of the level set function $\phi:\Omega\times[0,T]\to\mathbb{R}$:

[TABLE]

with $|\nabla\phi|>0$ in $\mathcal{O}(\mathcal{S})$. Note that $\phi$ can be defined as the solution to the level-set equation [56]:

[TABLE]

where $\mathbf{u}^{\rm bulk}$ is a given velocity field in $\Omega$ such that $\mathbf{u}^{\rm bulk}\cdot\mathbf{n}=\mathbf{u}\cdot\mathbf{n}$ on $\Gamma(t)$. For the purpose of numerical integration, we approximate $\Gamma$ with a “discrete” surface $\Gamma_{h}$, which is defined as the zero level set of a $P_{1}$ Lagrangian interpolant for $\phi$ on the one time refined mesh:

[TABLE]

where $I^{2}_{h/2}(\phi)$ is a $P_{2}$ nodal interpolant of the level set function. See Fig. 2. Note that $\mathbf{n}_{h}$ is defined in a neighborhood of $\Gamma_{h}$. The discrete surface $\Gamma_{h}$ satisfies the following geometric approximation properties:

[TABLE]

The analysis of the TraceFEM for scalar elliptic problem on a steady surface [49] and parabolic problems on evolving surfaces [35] shows that the geometric consistency error (45) allows to prove optimal order convergence for $P^{1}$ elements (both in energy and $L^{2}$ norms). The full convergence analysis of TraceFEM for fourth order problems as the surface Cahn-Hilliard equation (11) is an open problem that we shall address elsewhere. When only the initial position of the surface and the bulk advection field $\mathbf{u}^{\rm bulk}$ are given, the evolution of $\Gamma$ can be recovered solving eq. (44) numerically. Then, the calculated solution $\phi_{h}$ defines $\Gamma_{h}$.

The numerical method approximates the solution to the surface Cahn-Hilliard problem and its extension to a neighborhood of the surface. In practice, at time $t^{n}$ we extend the current solution to a narrow band of $\Gamma^{n}_{h}$ such that $\Gamma^{n+1}_{h}$ stays inside this narrow band and so all discrete quantities at $t^{n+1}$ are computable. This narrow band $\mathcal{O}(\Gamma^{n}_{h})$ is defined as the union of tetrahedra within distance $\delta_{n}$ from $\Gamma^{n}_{h}$:

[TABLE]

To ensure

[TABLE]

we set the minimum thickness of the extension layer to be

[TABLE]

where $c_{\delta}\geq 1$ is an $O(1)$ mesh-independent constant. If $\phi_{h}$ is an approximate distance function, then the condition $\mbox{dist}(S,\Gamma_{h})\leq\delta_{n}$ in (46) can be replaced by $|\phi_{h}^{n}(\mathbf{x})|\leq\delta_{n}$ for any vertex $\mathbf{x}$ of $T$.

We define finite element spaces

[TABLE]

These spaces are the restrictions of the time-independent bulk space $V_{h}^{\rm bulk}$ on all tetrahedra from $\mathcal{O}(\Gamma^{n}_{h})$.

Our hybrid finite difference in time / finite element in space method is based on the semi-discrete formulation (21)–(23). It reads: Given $c_{h}^{0}\in V_{h}^{0}$, for $n=1,\dots,N$ find $c_{h}^{n}\in V_{h}^{n}$ and $\mu_{h}^{n}\in V_{h}^{n}$ satisfying

[TABLE]

for all $v_{h},\,q_{h}\in V_{h}^{n}$. Here $\mathbf{u}^{e}(\mathbf{x})=\mathbf{u}(\mathbf{p}^{n}(\mathbf{x}))$, i.e. the lifted data on $\Gamma^{n}_{h}$ from $\Gamma^{n}$. The first term in (50) is well-defined thanks to condition (47) with the index shifted $n\to n-1$. In accordance to the analysis of the scalar advection-diffusion problems [35], parameters $\rho_{\mu}$ and $\rho_{c}$ are set to be

[TABLE]

The role of the term $\rho_{c}\int_{\mathcal{O}(\Gamma^{n}_{h})}\frac{\partial c^{n}_{h}}{\partial\mathbf{n}_{h}}\frac{\partial q_{h}}{\partial\mathbf{n}_{h}}\,dx$ in eq. (51) is twofold. First, form

[TABLE]

is positive definite on the narrow-band finite element space $V_{h}^{n}$, rather than only on the space of traces. Therefore, at each time step we obtain a finite element solution $c^{n}_{h}$ defined in $\mathcal{O}(\Gamma^{n}_{h})$. This can be seen as an implicit extension procedure of finite element solution from $\Gamma^{n}_{h}$ to the narrow band $\mathcal{O}(\Gamma^{n}_{h})$. The same observation holds for $\mu^{n}_{h}$ and the corresponding term in eq. (50). Furthermore, adding two volume terms (one for $c^{n}_{h}$ and another for $\mu^{n}_{h}$) makes the problem algebraically stable, i.e. the condition numbers of the resulting matrices are independent on how the surface cuts through the background mesh. Actually, the algebraic stabilization was the original motivation of introducing such volumetric terms in [7, 24] for unfitted surface finite element methods.

Remark 3.1** (Implementation).**

For the realization of the method, one uses the standard nodal basis functions for the bulk volumetric mesh $\mathcal{T}_{h}$. At time step $n$, only the degrees of freedom of the tetrahedra in the narrow band $\mathcal{O}(\Gamma^{n}_{h})$ are active. Since $\Gamma_{h}^{n}$ is the zero level of the $P_{1}$ finite element function $\phi_{h}^{n}$, the tetrahedra intersected by $\Gamma_{h}^{n}$ are those $T\in\mathcal{O}(\Gamma^{n}_{h})$ for which $\phi_{h}^{n}$ has a change in sign at different vertices. Then, $\Gamma_{h}^{n}\cap T$ is either a triangle or a flat rectangle, which can be further divided in two triangles. The vertexes of these triangles are immediately available from the nodal values of $\phi_{h}^{n}$ in $T$. So, it is straightforward to apply standard quadrature rules to compute the surface integrals in (50)–(51). We note that on hyperrectangles the implicit representation of $\Gamma_{h}$ by the bilinear level-set function can be treated with a marching cube [39] approximation.

Remark 3.2** (Higher order FEM).**

Although this paper discusses only the finite element method based on polynomials of degree 1, higher order trace finite elements are possible. To gain full benefits of using higher order elements, one should ensure that the geometric and finite element interpolation errors are of the same order. One way to achieve this, is to use higher order polynomials to define the discrete level set function $\phi_{h}$ that implicitly defines $\Gamma_{h}$. However, it is a non-trivial task to obtain a parametrized representation of the zero level of $\phi_{h}$ (for polynomial degree $\geq 2$) which allows for a straightforward application of numerical quadrature rules. The higher order case requires special approaches for the construction of quadrature rules as discussed, for example, in [21, 45, 53, 55, 62]. Isoparametric unfitted finite element method is an another elegant and efficient way to simultaneously leverage the geometric and interpolation accuracy of TraceFEM; see [34, 24].

Remark 3.3** (Methods based on surface PDEs extension).**

The finite element method (50)–(51) is different from a surface unfitted FEM based on the extension of a PDE from the surface to an ambient bulk domain [6, 26, 11, 52]. Both methods avoid surface triangulation and remeshing (if the surface evolves). However, in the method originally suggested in [6] a PDE is first extended and then solved in a one-dimension higher domain (the challenges of such approach are discussed and partially addressed in [26, 52]), while the present method extends only the solution, i.e. a function rather than a PDE is extended to the narrow band. In the TraceFEM, the PDE is formulated on the original surface.

4 Numerical examples

We assess the model introduced in Sec. 2 and the method presented in Sec. 3 on a series of benchmark tests in Sec. 4.1. In Sec. 4.2, we apply our approach to study spinodal decomposition and pattern formation on colliding surfaces. Finally, in Sec. 4.3 we consider the phase separation on a sphere splitting into two droplets.

For the numerical tests, we will consider mobility $M$ as in (15) with $\sigma=1$ (the only exception is test 2c, where we set $\sigma=16$) and free energy (16). Furthermore, we always use constant density $\rho=1$. Although, there are no reasons to assume $\rho=\mbox{const}$ for motions that violate the inextensibility condition (cf. (5)), we expect that moderate smooth variations of $\rho$ would lead to minor changes in lateral phase dynamics. In addition, $\rho=1$ is a fair assumption for the assessment of the numerical method and the model. We set $c_{\delta}=1$ in (48) to define the extension strip. Implementation of the method was been done using the FE package DROPS [13].

4.1 Validation

To study the accuracy of the finite element method described in Sec. 3, we consider two sets of benchmark tests. The first set features a sphere undergoing rigid body motion and is presented in Sec. 4.1.1. The second set, illustrated in Sec. 4.1.2, involves an ellipsoid that dilates and shrinks.

4.1.1 Sphere undergoing rigid body motion

We consider a unit sphere, initially centered at the origin. During the time interval of interest $[0,1]$, the sphere is translated with translational velocity $\mathbf{v}$ and rotated with angular velocity $\boldsymbol{\omega}$. We will report the numerical results for two tests:

Test 1a: $\mathbf{v}=(1,0,0)$ and $\boldsymbol{\omega}=(1,0,0)$. Its aim is to illustrate the methodology.

• -

Test 1b: $\mathbf{v}=(10,0,0)$ and $\boldsymbol{\omega}=(10,0,0)$. Its aim is to check the spatial and temporal accuracy.

In both test we take $\epsilon=0.1$.

We embed the sphere undergoing rigid body motion in outer domain

[TABLE]

The initial triangulation $\mathcal{T}_{h_{0}}$ of $\Omega$ consists of 12 sub-cubes, where each of the sub-cubes is further subdivided into 6 tetrahedra. In addition, the mesh is refined towards the surface, and $\ell\in\mathbb{N}$ denotes the level of refinement, with the associated mesh size $h_{\ell}=2^{-\ell-2}10/3$.

Test 1a. To better illustrate the methodology, we first consider a simple test case. We take the homogeneous equation Cahn–Hilliard equation. We set $c_{x_{2}}$ as the initial solution, where

[TABLE]

Here and in the following, $\mathbf{x}=(x_{1},x_{2},x_{3})^{T}$ denotes a point in $\mathbb{R}^{3}$. For this test we set $\delta=0$ (we will use $c_{x_{i}}$ with $\delta\neq 0$ and $i\neq 2$ later).

Fig. 3 displays the solution computed with mesh $\ell=5$ at times $t=0,0.5,1$, together with a view of the bulk mesh. Notice the mesh refinement in the neighborhood of the surface. This ensures that there are always enough elements to resolve the interface thickness $\epsilon$.

Test 1b. We consider the following synthetic solution to the Cahn–Hilliard equation:

[TABLE]

The above $c^{*}$ is the exact solution to the non-homogeneous equation

[TABLE]

The exact chemical potential $\mu^{*}$ can be readily computed from eq. (14). The non-zero right-hand side $g$ is calculated so that the solution follows the motion of the sphere without any change. To compute the right-hand side $g$ it is helpful to notice that $Y$ is a real spherical harmonic function.

We report in Fig. 4 (left) the discrete $L_{2}(0,T;L^{2}(\Gamma_{h}))$ norm:

[TABLE]

and the discrete $L_{2}(0,T;H^{1}(\Gamma_{h}))$ norm:

[TABLE]

of the error for order parameter $c$ (blue lines) and chemical potential $\mu$ (red lines) plotted against the refinement level $\ell$. The time step was refined together with the mesh size according to $\Delta t=4^{1-\ell}/10$. All the norms reported in Fig. 4 (left) are computed on the approximate surface $\Gamma_{h}$, where $c^{*}$ and $\mu^{*}$ were defined through their normal extensions from $\Gamma$. The observed second order convergence in the $L_{2}(0,T;L^{2}(\Gamma_{h}))$ for both concentration and chemical potential as well as the first order convergence in the $L_{2}(0,T;H^{1}(\Gamma_{h}))$ norm for the concentration are consistent with the well known error analysis of finite element methods for the Cahn–Hilliard equation [18, 19]. For the error in chemical potential we also observe the almost second order of convergence in $L_{2}(0,T;H^{1}(\Gamma_{h}))$ norm. We do not have an explanation for this apparent super-convergence.

4.1.2 Oscillating ellipsoid

The second series of tests is inspired by a numerical experiment in [19]. We consider time-dependent surface $\Gamma(t)=\{\mathbf{x}\in\mathbb{R}^{3}:\phi(\mathbf{x},t)=0\}$ with

[TABLE]

which is an ellipsoid centered at origin with variable length for the principal axis aligned with the $x_{1}$-axis. We embed this evolving surface in domain $\Omega=[-5/3,5/3]^{3}$. As initial solution to the Cahn–Hilliard equation, we take a small perturbation about 0.5 given by:

[TABLE]

We study the evolution of the discrete Lyapunov energy:

[TABLE]

for three datasets:

Test 2a: $k=1$, $\sigma=1$, $\epsilon=0.1$, $\Delta{}t=0.01$.

• -

Test 2b: $k=1$, $\sigma=1$, $\epsilon=0.01$, $\Delta{}t=0.01$.

• -

Test 2c: $k=5$, $\sigma=16$, $\epsilon=0.1$, $\Delta{}t=0.001$.

For all the tests, we set $T=1$. To mesh the domain we follow the same procedure used for the first sets of tests, i.e. we divide $\Omega$ into sub-cubes that are further divided into tetrahedra. We consider different levels of refinement $\ell$ associated with mesh size $h=2^{-\ell-2}10/3$.

Test 2a. We present the evolution of the discrete Lyapunov energy in Fig. 4 (right). We see that mesh $\ell=1$ is too coarse to correctly capture the evolution of $E^{L}_{h}(c_{h})$. The results obtained with mesh $\ell=2$ are very close to results obtained with meshes $\ell=3,4$, which are almost superimposed for the entire time interval under consideration. Fig. 4 (right) seems to suggest that that the solution converges to a time periodic solution. In order to verify that, we would need to run the simulation for a longer period of time. Instead, we prefer to decrease the value of $\epsilon$ in the next test. We recall that $\epsilon$ is a crucial modeling parameter, since it defines the thickness of the layer where phase transition takes place and also defines the intrinsic time scale. Thus, the case of smaller $\epsilon$ is numerically challenging.

Test 2b. We keep all the model and discretization parameters as in test 2a, but we reduce the value of $\epsilon$ to $0.01$. Fig. 5 (left) shows the evolution of $E^{L}_{h}(c_{h})$ over time. Since this test is more challenging, to observe convergence for the discrete Lyapunov energy it takes a higher level of refinement than test 2a. In fact, from Fig. 5 (left) we see that the results computed with mesh $\ell=4$ are still quite far from the results obtained with mesh $\ell=5$, which are almost superimposed to the results for mesh $\ell=6$. This indicates that both meshes $\ell=5$ and $\ell=6$ are sufficiently refined for $\epsilon=0.01$. Notice that the discrete Lyapunov energy computed with meshes $\ell=5,6$ decreases rapidly till around $t=0.3$ and then it gently decreases, although not monotonically. Moreover, we note that the energy may have locally increase driven by the evolution of the domain. This phenomenon will be even more pronounced in test 2c, which is presented next .

We show the evolution of the solution computed with mesh $\ell=6$ in Fig. 6. Around $t=0.3$, we observe one big pink domain (i.e., $c=1$) and several small black domains (i.e., $c=0$). For $t>0.3$, the aspect of the solution (one big pink domain and several smaller black domains) remains almost unchanged as the surface gets dilated and shrunk.

Test 2c. We present the evolution of the discrete Lyapunov energy in Fig. 5 (right). We see that the results computed with meshes $\ell=4,5$ are in excellent agreement. In [19], good convergence of the solution is obtained for short times, i.e. till around $t=0.2$. For the choice of parameters associated to test 2c, which are close to the parameter values in [19], we observe good convergence of the solution for the whole interval of time under consideration. We notice that the discrete Lyapunov energy computed with meshes $\ell=4,5$ decreases (not monotonically) till around $t=0.35$. For $t>0.35$, Fig. 5 (right) suggests that the solution converges to a time periodic solution.

We show the evolution of the solution computed with mesh $\ell=6$ in Fig. 7. Note that we use the truncated scale here for $c$ to better illustrate the solution, which is otherwise rather diffusive (since $\epsilon=0.1$ is used). From Fig. 7, we see that the solution configuration remains the same for $t>0.35$ as the surface gets dilated and shrunk. The patterns and the evolution look very similar to those in Figure 4 of [19].

4.2 Colliding spheres

We consider an evolving surface $\Gamma(t)$ that undergoes a topological change and experiences a local singularity. The dynamics of $\Gamma(t)$ is inspired by tests presented in [25, 54, 35]. The computational domain is $\Omega=[-10/3,10/3]\times[-5/3,5/3]^{2}$ and the evolving surface is the zero level set of level set function

[TABLE]

with

[TABLE]

for $t\in[0,1.5/|w|]$. The initial configuration $\Gamma(0)$ is close to two balls of radius 1, centered at $\mathbf{x}_{c}^{\pm}(0)$. Parameter $w$ is the collision speed. For $t>0$, $w>0$ the two spheres approach each other until time $\tilde{t}\approx 0.235/w$, when they touch at the origin; then, for $t\in(\tilde{t},1.5/w)$, surface $\Gamma(t)$ is simply connected.

In the vicinity of $\Gamma(t)$, the gradient $\nabla\phi$ and the time derivative $\partial_{t}\phi$ are well defined and given by simple algebraic expressions. The normal velocity field of $\Gamma(t)$ can be computed to be

[TABLE]

We assume the following material velocity:

[TABLE]

which models the parallel advection away from the merging line and nearly normal motion near the $x_{1}=0$ plane.

We consider different meshes for $\Omega$ with mesh size $h=2^{-\ell-2}10/3$, where $\ell$ is the refinement level.

The setup is a simple model of two plasma membrane fusion. In all further experiments we take $\epsilon=0.01$. We remark that this is a realistic value of interface thickness for applications related to phase separation in lipid bilayers. In fact, if we consider a typical giant vesicle with an average diameter of 30 $\mu$m, on which phase separation can be visualized using fluorescence microscopy [66] with a resolution of about 300 nm, the thickness of transition region between the phases is approximately 1% of vesicle diameter.

4.2.1 Colliding spheres with pre-separated phases

For this test, we set $w=1$. The configuration at $t=3/2$ is the ball centered around 0 with radius $2^{1/3}$. We end the simulation at $t=1$, when the two balls have merged but they have not evolved into one sphere yet.

As initial solution, we take $c_{x_{2}}$ with $\delta=0$ for the ball initially centered at $\mathbf{x}_{c}^{-}(0)$ and $c_{x_{1}}$ with $\delta=0$ for the ball initially centered at $\mathbf{x}_{c}^{+}(0)$. See (53) for the definition of $c_{x_{i}}$. We set $\Delta t=0.001$.

Fig. 8 displays the evolution of the surface and the solution for $t\in[0.05,1]$. Before the two balls touch at $t=\tilde{t}$, we see that the ball centered at $\mathbf{x}_{c}^{-}(t)$ (resp., $\mathbf{x}_{c}^{+}(t)$) has a horizontal (resp., vertical) interface separating a pink and black domain of equal size. For $t>\tilde{t}$, the black and pink domain of the ball centered at $\mathbf{x}_{c}^{-}(t)$ get in contact with the black domain of the ball centered at $\mathbf{x}_{c}^{+}(t)$. Between the pink domain of the ball centered at $\mathbf{x}_{c}^{-}(t)$ and the black domain of the ball centered at $\mathbf{x}_{c}^{+}(t)$ an interface gets formed. Such interface is virtually located on the curve of minimal length on the simply connected surface $\Gamma(t)$, for $t>\tilde{t}$. This is consistent with the well-known limiting (as $\epsilon\to 0$) behavior of the stationary Cahn-Hilliard problem that corresponds to solving the isoperimetric problem [44, 61].

4.2.2 Pattern formation on colliding sphere

Let $t_{s}=0.23/w<\tilde{t}$ be a time shift in order to have phase separation occur close to the surface collision. All the simulations whose results are reported in this section start at time $t_{s}$. We consider two tests:

Slow collision: $w=1$, $T=1+t_{s}$, $\Delta t=10^{-4}$ for $t\in(t_{s},t_{s}+0.01]$ and $\Delta t=10^{-2}$ for $t\in(t_{s}+0.01,t_{s}+1]$.

• -

Fast collision: $w=10$, $T=0.1+t_{s}$, $\Delta t=10^{-4}$ for $t\in(t_{s},t_{s}+0.01]$ and $\Delta t=10^{-3}$ for $t\in(t_{s}+0.01,t_{s}+0.1]$.

The two initial spheres touch each other at the origin at time $\tilde{t}\approx{}0.235$ for the slow collision test and $\tilde{t}\approx{}0.0235$ for the fast collision test. Let $\text{rand}()$ be a uniformly distributed random number between 0 and 1. For both tests, we take $c_{0}(\mathbf{x})=\text{rand}()$ for the ball initially centered at $\mathbf{x}_{c}^{-}(0)$ and $c_{0}=0$ for the ball initially centered at $\mathbf{x}_{c}^{+}(0)$.

As we have observed in our previous work [69] and is well known (see, e.g.,[28, 57]), the evolution of the solution to the Cahn–Hilliard problem goes through an initial fast phase, during which a pattern gets formed, followed by a slowdown in the process of dissipation of the interfacial energy. In order to capture the different time scales, we prescribe different time steps for the different stages. It would be less intrusive to use some time-adaptivity strategy, which is not addressed in this paper.

Fig. 9 shows the evolution of the surface and the solution for the slow collision test for $t\in[t_{s},t_{s}+1]$. We see the quick formation of a pattern, that is visible already at $t=t_{s}+0.005$. Thus, when the two balls touch at $\tilde{t}\approx{}0.235$, the already phase separated solution of the ball centered at $\mathbf{x}_{c}^{-}(t)$ enters in contact with the black domain that covers the entire ball centered at $\mathbf{x}_{c}^{+}(t)$. Similarly to what observed in Sec. 4.2.1, for $t>\tilde{t}$ we see that the interfaces between the pink domains and the large black domain tend to align with the curve of minimal length on surface $\Gamma(t)$. Hence, the pink domains remain confined on the portion of the surface in the $x_{1}<0$ semi-space. Possible factors that come into play to determine the appearance of the pattern on $\Gamma(t)$ are:

The droplet shape of the balls before merging and the bottleneck shape, with a short perimeter cross-section, after merging. 2. 2.

The contact with the black domain covering the ball centered at $\mathbf{x}_{c}^{+}(t)$. 3. 3.

The colliding dynamics.

To understand the role of factors 1. and 2., we will study phase separation on two surfaces frozen in time (not evolving) relative to the slow collision test (i.e., $w=1$ and $t_{s}=0.23$): half of domain $\Gamma(t_{s})$, hereafter called single ball, and surface $\Gamma(t_{s}+0.02)$, hereafter called hourglass. The fast collision test, i.e. $w=10$, is meant to understand the role of factor 3.

We start by comparing phase separation on the single ball and on the hourglass. The initial condition for the single ball is simply $c_{0}(\mathbf{x})=\text{rand}()$ for the entire domain. The initial condition for simulation on the hourglass mimics the initial condition for the slow collision test, i.e. $c_{0}(\mathbf{x})=\text{rand}()$ for $\mathbf{x}$ with $x_{1}\leq 0$ and $c_{0}=0$ for $\mathbf{x}$ with $x_{1}>0$. We consider the same meshes described in Sec. 4.2.

In Fig. 10, we see the results for both simulations for $t\in[0,1]$. We observe elongated pink domains aligned with the $x_{1}$-axis, i.e. the horizontal axis in the figure, for both cases. Since both surfaces are not evolving, we believe this patter is solely due to the shapes (droplet and bottleneck). On the other hand, the interaction with the large black domain on the hourglass does not seem to significantly influence the pattern for $t<1$.

Next, we present the results for the fast collision test. We report in Fig. 11 the evolution of the surface and the solution for $t\in[t_{s},t_{s}+0.1]$. We see that, due to the fast collision dynamics, the two balls touch before phase separation occurs on the balls centered at $\mathbf{x}_{c}^{-}(t)$. Upon phase separation, we see the same elongated pink domains aligned with the $x_{1}$-axis already observed in Fig. 9 and 10. The notable difference between the corresponding solutions in Fig. 9 and 11 (taking into account the 10 times faster dynamics in the latter) are the thicker interfaces along the bottleneck in Fig. 11, where we see larger blue regions in the $t=t_{s}+0.02,t_{s}+0.04,t_{s}+0.08$ panels, and the longer portions of interface on the bottleneck in Fig. 11, $t=t_{s}+0.1$ panel, which shows a ‘mushroom’ shape of the pink fingers comparing to a more stretched, oval shape for slow collision motion.

For further comparison, we report in Fig. 12 a side view of the solution at two different times for the single ball (top), slow collision (center), and fast collision (bottom) tests. We selected times $t=0.4,1$ for the single ball and slow collision tests (same stage of phase separation) and times $t=0.04,0.1$ for the fast collision test, which feature the same surface as the slow collision test for $t=0.4,1$. Notice that the solution for the fast collision test is at an earlier stage of phase separation. Given this and the fact that the simulations are started from a random initial condition, we conclude that the collision speed does not have a substantial influence on the pattern appearance at the back side of colliding drops.

4.3 Splitting spheres

We consider an evolving surface undergoing a reverse dynamics with respect to the one considered in Sec. 4.2: the initial configuration is a sphere centered at the origin, which splits into two droplets, in turn evolving towards spheres centered at $\mathbf{x}_{c}^{\pm}(T)$. The computational domain is again $\Omega=[-10/3,10/3]\times[-5/3,5/3]^{2}$. The evolving surface $\Gamma(t)$ is the zero level set of $\phi(\mathbf{x},1.5/|w|+t)$ with $w=-1$, where $\phi$ is given by (57). The velocity vector field is defined by (59) with $w=-1$. Surface $\Gamma(t)$ is simply connected for $t<1.5-\tilde{t}$, where $\tilde{t}=0.235$. Obviously, this evolving surface experiences a local singularity too.

We consider the same meshes used for the results reported in Sec. 4.2. Below we present the results for 2 numerical experiments:

Rotated interface: with initial solution is given by $c_{0}^{x_{3}}$ from (53), with $\delta=0$ and rotated by $\pi/4$ in $x_{1}x_{2}$ plane.

• -

Horizontal interface: with initial solution $c_{0}^{x_{3}}$ from (53), with $\delta=20$.

Fig. 13 and 14 show the evolution of the interface and the solution for $t\in[0,1.5]$ for the rotated and horizontal interface tests, respectively. We see that as the two hemispheres are being pulled in opposite directions along the $x_{1}$-axis, a horizontal interface forms between the pink domain and the black domain in both figures. In the case of the rotated interface test, the horizontal portion of the interface is located symmetrically with respect to the $x_{1}$-axis. Thus, it does not get affected by the bottleneck formation. On the other hand, the horizontal interface in Fig. 14 is located non-symmetrically with respect to the $x_{1}$-axis. Hence, when the bottleneck forms the interface curves. As a result, the two balls into which the surface has evolved at $t=1.5$ present a curved interface. See Fig. 14, bottom right panel.

5 Conclusions

We presented a formulation of the Cahn–Hilliard equation on a time-dependent surface $\Gamma(t)$ that uses tangential calculus induced by the embedding of $\Gamma(t)$ in $\mathbb{R}^{3}$. This makes the formulation particularly suitable for the development of fully Eulerian numerical techniques for the solution of the problem. The fully Eulerian technique proposed in this paper is based on the unfitted trace finite element method for spatial discretization and the finite difference method combined with an implicit extension procedure for time discretization. The degrees of freedom are tailored to a background time-independent mesh and we used the usual nodal basis to build the systems of algebraic equations at each time step; space-time integration is avoided. All these ingredients lead to a rather straightforward implementation of the method in a standard finite element software. Numerical experiments demonstrated the following main properties of the proposed approach: (i) optimal second order accuracy for $P_{1}$ approximation of the concentration, chemical potential, and the level set function that defines the surface; (ii) numerical stability even in the case of fast and large deformations; (iii) the ease of handling phase transition on a surface undergoing topological changes. The third feature was previously available only for methods based on phase-field representation of the surface itself.

There are many directions to further develop the proposed method. Some of them are: employing $P^{k}$ ($k>1$) trace finite elements together with higher order isoparametric surface recovery; the coupling of lateral phase separation with the surface evolution through the line tension forces; accounting for the surface fluidity and bulk phenomena. We believe that all these extensions fit well within the proposed framework and we plan to address some of them in the future.