Numerical modelling of phase separation on dynamic surfaces
Vladimir Yushutin Annalisa Quaini Maxim Olshanskii
math.NA cs.NA physics.bio-ph physics.comp-ph 65M60, 35Q79, 35Q92



Numerical modelling of phase separation on dynamic surfaces
    
Vladimir YushutinDepartment of Mathematics, University of Houston, Houston, Texas 77204 (yushutin@math.uh.edu).
Annalisa QuainiDepartment of Mathematics, University of Houston, Houston, Texas 77204 (quaini@math.uh.edu); Partially supported by NSF through grant  DMS-1620384 and by the  Bridge Funding Program of University of Houston.
Maxim OlshanskiiDepartment of Mathematics, University of Houston, Houston, Texas 77204 (molshan@math.uh.edu); Partially supported by NSF through grant  DMS-1717516.

    
=========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================



 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 ℝ^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.


Surface Cahn–Hilliard equation;  Evolving interfaces; TraceFEM; Membrane fusion.
 






§ 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 <cit.>.
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 <cit.>,  dynamics of  plasma membranes and  multicomponent vesicles
<cit.>, and
lipid rafts formation <cit.>. 
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
<cit.> 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 <cit.> and on more general
stationary domains in <cit.>. A finite difference method for a diffuse volumetric representation of the surface
CH equation was introduced in <cit.>.
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 <cit.>, where numerical examples
for a sphere and saddle surface are provided; results obtained by FEM on more general surfaces
are presented in <cit.>.
In <cit.>, 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 <cit.>, 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 <cit.>, 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 <cit.>.
The FEM used for the space discretization in <cit.>, 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 <cit.>.
The authors of <cit.> 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 <cit.>, 
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
<cit.>.
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 <cit.>, 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 <cit.>:
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 <cit.>, which is explained in Remark <ref>.
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 <cit.>
and evolving surfaces in <cit.>.
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 <cit.>.
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 <cit.>.
The numerical experiments presented in Sec. <ref> mimic phase separation occurring on a two component
lipid vesicle, leading to fusion with another vesicle. These and other numerical results in
Sec. <ref> 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 <ref>,
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 <ref>.
After the assessment of the numerical method for a set of model problems, in Section <ref>
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 <ref> closes the paper with a few concluding remarks.



§ MATHEMATICAL MODEL




 §.§ Preliminaries


Let Γ(t) be a closed, smooth, evolving surface in ℝ^3
for t ∈ [0, T], where T > 0 is a final time.
Let the 3D bulk domain Ω be such that Γ(t)⊂Ω, for all t ∈ [0, T].


We assume surface Γ(t) evolves according to a given smooth velocity field ,  i.e., (,t) is the velocity of point ∈Γ(t).  Consider the decomposition of  into normal (_N=u_N) and tangential (_T) components:

    = _T + u_N,  _T·=0,

where  is the outward normal vector on Γ.
While the normal velocity u_N completely defines the geometric evolution of the surface, we are interested in the motion
of Γ(t) as a material surface, which also includes tangential deformations. The material deformation  can be defined through the Lagrangian mapping Ψ(t, ·) from Γ(0) to Γ(t), i.e. for ∈Γ(0), Ψ(t,) solves the ODE system

    Ψ(0,)=,  ∂Ψ(t,)/∂ t=(t,Ψ(t,)),   t∈[0,T].



For any sufficiently smooth function f in a neighborhood of Γ(t) its tangential gradient is
defined as f=∇ f-(·∇ f). The tangential (surface) gradient f
then depends only on values of f restricted to Γ(t) and · f=0 holds.
For a vector field  on Γ(t) we define  componentwise.
The surface divergence operator for  and the Laplace–Beltrami operator for f
are given by:


    :=  (),     f :=  ( f),

where (·) is the trace of a matrix.
For a fixed t∈[0,T], L^2(Γ(t)) is the Lebesgue space of square-integrable functions on Γ(t) and
H^1(Γ(t)) is the Sobolev space of all functions f∈ L^2(Γ(t)) such that ∇_Γ f∈ L^2(Γ(t))^3.
For a subdomain S(t)⊂Γ(t), we recall the integration by parts identity:
 
    ∫_S(t) f÷_Γ  ds=∫_∂ S(t) f·  dγ -∫_S(t)·∇_Γ f  ds+∫_S(t)κ f·  ds,

where κ is the sum of principle curvatures and  is the normal vector on ∂ S(t) that is tangential to Γ(t) and outward for S(t);  vector field  and scalar function f are such that all
the quantities in (<ref>) exist.

Let  S(t)⊂Γ(t) be a material area evolving according to (<ref>). Then, the following surface analogue of the Reynolds transport theorem holds (see, e.g., Lemma 2.1 in
<cit.>):

    d/ dt∫_S(t) f    ds =  ∫_S(t)(ḟ + f )   ds,

where ḟ is the material derivative of  f.
One can write the material derivative in terms of partial derivatives

    ḟ = ∂ f/∂ t + ·∇ f.

The terms on the right-hand side of (<ref>) are well defined if one identifies f with its arbitrary smooth extension
from the surface to  a neighborhood of the space-time manifold ⊂ℝ^4,

    =⋃_t∈[0,T]Γ(t)×{t}.

Of course, ḟ is an intrinsic surface quantity and the quantity on right-hand side of (<ref>) 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 Γ.





 §.§ Model setup

On Γ(t) we consider a heterogeneous mixture of two species with densities ρ_i, i = 1, 2.
Let ρ=ρ_1+ρ_2 be the total density. The conservation of mass  written for an arbitrary material area  S(t)⊂Γ(t) and (<ref>) yield

    0= d/ dt∫_S(t)ρ  ds =∫_S(t)(ρ̇ + ρ)   ds.

Since the above identity holds for  arbitrary   S(t)⊂Γ(t), it implies

    ρ̇ + ρ=0  on Γ(t).

Following <cit.>,
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∈ [0,1]. 
Mass conservation  for one component on S(t)⊂Γ(t) takes the form

    d/ dt∫_S(t)ρ c    ds =    -∫_∂ S(t)·  dγ.

where  is a mass flux.

By applying the transport formula  (<ref>) to the left-hand side of (<ref>)  and the integration by
parts formula (<ref>) to the right-hand side of (<ref>), we obtain


    ∫_S(t) (ρ̇ c+ ċρ  +cρ)   ds = - ∫_S(t)  ds

Thanks to (<ref>), the above equality simplifies to

    ∫_S(t)ċρ  ds = - ∫_S(t)  ds

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

    ċ   +  1/ρ = 0   on Γ(t),



The classical assumption for the flux  is  Fick's law:

    = - M μ  on Γ,   μ = δ f/δ c,

where M is the so-called mobility coefficient (see <cit.>) and
μ 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 <cit.>, we set

    f(c) = 1/ϵf_0(c) + ϵ/2  |  c |^2.

In eq. (<ref>), 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 (<ref>)–(<ref>),
the interface between the two components is a layer of size ϵ.

By combining eqs. (<ref>), (<ref>)–(<ref>), we obtain
the Cahn–Hilliard equation posed on evolving material surface Γ(t):









    ρ̇ + ρ   =0   on Γ(t), 
    ċ  -  ρ^-1(M (1/ϵf_0' - ϵΔ_Γ c))    = 0  on Γ(t),

for all t ∈ (0, T).
We close system (<ref>)–(<ref>) with initial conditions

    ρ(·, t) = ρ_0,  and   c(·, t) = c_0    on Γ(0).

There are no boundary conditions since the boundary of Γ(t) is empty.

Note that eq. (<ref>) is decoupled from eq. (<ref>),
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. _T=0, eq. (<ref>)  simplifies to
 
    ρ̇ = -κ u_Nρ.

For inextensible membranes, it holds =0 and eq. (<ref>) further simplifies  to ρ̇=0.


The system (<ref>)–(<ref>) is similar to the Cahn-Hilliard equation on evolving surface found in <cit.>,
but different. In <cit.>, the equation is written for the conserved variable u=ρ  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  μ = δ f/δ u.
Formally, system (<ref>)–(<ref>) and the surface Cahn–Hilliard equation in <cit.> are
the same problem for inextensible membranes and if one assumes ρ_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., <cit.> for further discussion and references.





 §.§ Weak form of the equations

Eq. (<ref>) is a transport problem for the total density.
It can be integrated independently of eq. (<ref>) as long as the evolution of the surface
is given. Since ρ is completely determined by the initial data and , but not by the phase composition of Γ(t),
for the purposes of this paper we shall assume that ρ is given.
Eq. (<ref>) 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. (<ref>) in mixed form, i.e. as two coupled second-order equations:

    ċ -  ρ^-1(M μ)  = 0   on Γ(t) , 
       μ = 1/ϵf_0' - ϵΔ_Γ c   on Γ(t).


System (<ref>)–(<ref>) 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

    M = M(c) = σ c(1 - c),

with σ>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

    f_0(c) = 1/4 c^2(1 - c)^2.



We are interested in a finite element numerical method for the evolving surface Cahn–Hilliard problem
(<ref>)–(<ref>). We first  need a weak  formulation of system
(<ref>)–(<ref>).
To devise it, we start with multiplying (<ref>) by ρ
and a smooth v:→, while we multiply (<ref>)
by a smooth q:→. Next, we integrate  over Γ(t) and employ
the integration by parts identity (<ref>) with S(t)=Γ(t) (this implies ∂ S(t)=∅). The integration by parts is applied to the diffusion terms in (<ref>) and (<ref>).
The curvature terms vanish  since the fluxes M ∇_Γμ and ∇_Γ c are tangential. We obtain the following equalities:

    ∫_Γ(t)ρċ v   ds +  ∫_Γ(t) M μ  v   ds = 0, 
       ∫_Γ(t)μ q   ds - 1/ϵ∫_Γ(t) f_0'(c)  q   ds -  ϵ∫_Γ(t) c   q   ds = 0.


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 : L^q(), 1≤ q≤∞ and H^1() are standard Lebesgue and Sobolev spaces on . We also need

    H^1,Γ   ={v∈ L^2() :  v_L^2()≤∞},  
    
      L^∞_1,Γ   ={v∈ L^∞() : sup_t∈[0,T] v_L^2()≤∞}.

Note that  ∇_  and ∫_ can be written in terms of tangential calculus on Γ and the geometric motion of Γ, i.e. u_N and . See  <cit.> for details and properties of the above spaces.

We introduce the week formulation of (<ref>)–(<ref>):


 Find c∈ L^∞_1,Γ∩ H^1() and μ∈ H^1,Γ satisfying (<ref>)–(<ref>) for almost every t∈(0,T) and for all (v,q) ∈ H^1,Γ× H^1,Γ.

For a closely related problem (which coincides with (<ref>)–(<ref>) if ÷_Γ=0, ρ=1, M=), the above weak formulation was shown in <cit.> to have the unique solution.  Extending the well-posedness analysis for (<ref>)–(<ref>) is out of the scope for this paper.




§ A HYBRID FINITE DIFFERENCE / FINITE ELEMENT NUMERICAL METHOD


In this section, we present a completely Eulerian numerical method for solving problem (<ref>)–(<ref>).
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 <cit.>.
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
(<ref>)–(<ref>), 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 <cit.>
and suggest a  hybrid,  finite-difference in time and finite element in space discretization of
problem (<ref>)–(<ref>).

We start with several important observations about the differential problem  (<ref>)–(<ref>) and its alternative integral formulation.



 §.§ Extending c and μ off the surface
 We are  able to decompose the material derivative into the sum of Eulerian terms,

    ċ= ∂ c/∂ t + ·∇ c,

if we assume  an arbitrary smooth extension of c to  a neighborhood of , denoted further by (). Let us define ()=⋃_t∈(0,T)_δ(Γ(t))×{t}, where

    _δ(Γ(t))={∈^3 : (,Γ(t))<δ}.

We let δ be sufficiently small such that _δ(Γ(t))⊂Ω for all t∈(0,T).

Among infinitely many  possible smooth extensions of c to (),
it will be convenient to assume the one given by the close point projection on Γ(t):
Fix t∈(0,T), then for  ∈_δ(Γ(t)) we denote its close point projection on Γ(t) by (). For smoothly evolving Γ(t),
the projection () is well defined in (), and we consider the extensions given by
 
 c^e(,t)=c((),t) and μ^e(,t)=μ((),t) for (,t)∈().
 
We also introduce the extension of the spatial normal vector field, ^e(,t)=((),t).
Recall that (·,t) is the outward normal for Γ(t).
Now c^e and μ^e can be characterized as extensions along the directions given by the  normal field,
i.e. the so called normal extensions:
  
    ∂ c^e/∂^e= ∂μ^e/∂^e=0  in ().

If Γ(t) is a C^2 surface, then its normal field is C^1-smooth and c∈ H^1(Γ(t))
implies c^e∈ H^1(_δ(Γ(t))). In turn, the regularity of Γ(t) together
with assumptions on its evolution (e.g., in terms of smoothness of the mapping Ψ)
yield the space-time regularity of c^e (if c is smooth); see more details in, e.g.,
<cit.>.

We shall identify functions defined on  with their normal extensions and skip the upper index e. Once we do this,  the weak formulation  (<ref>)–(<ref>) yields the following identities for c and μ:

    ∫_Γ(t)ρ (∂ c/∂ t + ·∇ c)  v   ds +  ∫_Γ(t) M μ  v   ds
    
       = 0, 
    ∫_Γ(t)μ q   ds - 1/ϵ∫_Γ(t) f_0'(c)  q   ds -  ϵ∫_Γ(t) c   q   ds
    
       = 0,
    ∂μ/∂ = ∂ c/∂   =0  in _δ(Γ(t)),

for all smooth v and q defined in Γ(t). 
We shall base our discretization method on equalities  (<ref>)–(<ref>). The crucial idea is that our numerical solution should approximate the concentration and chemical potential on Γ(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.



 §.§ 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 <cit.>  and often leads at the end to the same fully discrete scheme.
To this end, consider the uniform time step Δ t=T/N, and let  t_n=nΔ t and I_n=[t_n-1,t_n).
Denote by c^n, μ^n an approximation of c(t_n) and μ(t_n), define Γ_n := Γ(t_n), n=0,…,N.
We assume that () is a sufficiently large neighborhood of  such that

    Γ_n⊂_δ(Γ_n-1)  for n=1,…,N;

see Fig. <ref>.
In this case, c^n-1 is well-defined on Γ_n and we can approximate the time derivative by simple finite difference

    ∂ c/∂ t(t_n)≃[c]_t^n:=c^n-c^n-1/Δ t.




This brings us to the implicit Euler method for problem (<ref>)–(<ref>):


    ∫_Γ_nρ^n ([c]_t^n + ^n ·∇ c^n)  v   ds +  ∫_Γ^n M^nμ^n   v   ds
       = 0, 
    ∫_Γ_nμ^n    q   ds  -  ϵ∫_Γ_n c^n   q   ds-
    1/ϵ∫_Γ_n f_0'(c^n)  q   ds    = 0,
    ∂μ^n/∂ = ∂ c^n/∂ =0  in    _δ(Γ_n),

It is important that all terms in (<ref>)–(<ref>) are defined on the
`current' surface Γ_n and its neighborhood.
In particular, c^n-1 is well-defined on Γ_n thanks to condition (<ref>),
which gives the meaning to [c]_t^n in (<ref>)–(<ref>),



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

    ∫_Γ_nρ^n ([c]_t^n + ^n ·∇ c^n)  v   ds +  ∫_Γ^n M^n-1μ^n   v   ds 
       = 0, 
    ∫_Γ_n  (μ^n - β_s Δ t[c]_t^n)   q   ds  -  ϵ∫_Γ_n c^n   q   ds
    
       = 1/ϵ∫_Γ_n f_0'(c^n-1)  q   ds,
    ∂μ^n/∂ = ∂ c^n/∂   =0  in _δ(Γ_n).

Following the idea of first order stabilization in <cit.>, we included
term -β_s Δ t[c]_t^n in eq. (<ref>).
 On a stationary domain, for constant mobility and a slightly modified double-well potential f_0
 than (<ref>), the scheme above was shown in <cit.> to be
 stable for large enough β_s=O(ϵ^-1).
In the semi-implicit Euler method both nonlinear terms,
i.e. the mobility coefficient in eq. (<ref>) and the potential
in eq. (<ref>), are extrapolated from the previous time step.
So, technically at each time step the method requires solving a linear system
(<ref>)–(<ref>) and further extending the solution along normals,
i.e. solving (<ref>) with given data on Γ_n.
We will show in Sec. <ref>
that both steps (solving for μ^n and c^n on Γ_n and extension)
can be naturally combined in one linear solve on the finite element level.




 §.§ 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 <ref>,
we shall illustrate by numerical examples that the evolution
of Γ(t) can actually produce a (local) increase of the system energy.

It is natural to say that (<ref>)–(<ref>) or
(<ref>)–(<ref>) are numerically stable
if the solution satisfy an energy bound  uniform in the discretization parameter Δ t.
In this Section, we demonstrate such a bound for (<ref>)–(<ref>)
subject to the following simplifications. First, we assume constant density and mobility: ρ^n=1, M^n=1.
We also assume a tangential incompressibility condition:
 
    ÷_Γ(,t)=0,  ∀ ∈Γ(t),

 which corresponds to the motion of an inextensible membrane.


We start by taking v=Δ t μ^n in (<ref>) and q=-Δ t([c]_t^n+ ^n ·∇ c^n) in (<ref>) and sum the resulting equalities.
After few cancellations, we obtain:

    Δ tμ^n^2_Γ_n   +β_s c^n-c^n-1^2_Γ_n
     + β_sΔ t ∫_Γ_n (c^n-c^n-1)(^n ·∇ c^n)  ds    + ϵ/2( c^n^2_Γ_n - c^n-1^2_Γ_n+ (c^n-c^n-1)^2_Γ_n)    +ϵΔ t∫_Γ_n c^n   (^n ·∇ c^n)   ds   = -1/ϵ∫_Γ_nf_0'(c^n)  (c^n-c^n-1)   ds -Δ t/ϵ∫_Γ_nf_0'(c^n)  (^n ·∇ c^n)   ds,

where ·_Γ_n stands for the L^2(Γ_n) norm.
First, we treat the terms with f_0' on the right-hand side.
Here and throughout the rest of this section, we use the fact that (<ref>)
yields ∇ c^n=∇_Γ c^n, ∇f_0(c^n)=f_0(c^n), etc.
Using integration by parts, formula (<ref>), and condition (<ref>), we have:

    ∫_Γ_nf_0'(c^n)  (^n ·∇ c^n)   ds=∫_Γ_n^n ·f_0(c^n)    ds = -∫_Γ_n (÷_Γ^n)f_0(c^n)  ds=0.

With the help of trunceted Taylor expansion and f_0”≥-1/4, we get

    -∫_Γ_nf_0'(c^n)  (c^n-c^n-1)   ds   = -∫_Γ_n (f_0(c^n) - f_0(c^n-1))   ds       -1/2∫_Γ_nf_0”(ξ^n)  |c^n-c^n-1|^2   ds   ≤  -∫_Γ_n (f_0(c^n) - f_0(c^n-1))   ds + 1/4c^n-c^n-1^2_Γ_n.





Next, we consider the terms on the left-hand side of (<ref>).
Using Cauchy–Schwarz inequality and estimate
ab≤1/2δa^2 + δ/2b^2, ∀ a,b∈ℝ with δ∈(0,+∞),
we obtain:

    β_sΔ t ∫_Γ_n (c^n-c^n-1)(^n ·∇ c^n)  ds ≤β_s/2c^n-c^n-1^2_Γ_n + C |Δ t|^2/ϵ c^n^2_Γ_n,

 where  constant C depends on 𝐮(,t).
The last term on the left-hand side of  (<ref>) is handled by the integration by parts

    ϵΔ t∫_Γ_n c^n  (^n ·∇ c^n)  ds
    =ϵΔ t∫_Γ_n c^n· [^n] c^n  ds       + ϵΔ t∫_Γ_n c^n (^n · )  c^n  ds
    = ϵΔ t∫_Γ_n c^n· [^n] c^n  ds      ≤  CϵΔ t c^n_Γ^n^2.

Combining (<ref>)–(<ref>), for β_s>1/2ϵ we get


    ϵ/2 c^n^2_Γ_n + 1/ϵ∫_Γ_nf_0(c^n)   ds +Δ tμ^n^2_Γ_n≤ϵ/2 c^n-1^2_Γ_n+1/ϵ∫_Γ_nf_0(c^n-1)  ds
    
         +CΔ t(ϵ+|Δ t|/ϵ) c^n_Γ^n^2,

where the constant C depends on 𝐮(,t).

Before applying the discrete Gronwall inequality to pass from (<ref>) to a priori stability estimate, we need to relate
the quantities c^n-1^2_Γ_n and ∫_Γ_nf_0(c^n-1)  ds to c^n-1^2_Γ_n-1 and ∫_Γ_n-1f_0(c^n-1)  ds, respectively.
To this end, consider the closest point projection ^n-1:
_δ(Γ_n-1)→Γ_n-1. Since c^n-1 satisfies (<ref>), we can write c^n-1()=c^n-1(^n-1()) in _δ(Γ_n-1). In particular, due to (<ref>)
this representation holds  for ∈Γ_n.
For the surface measures on Γ_n-1 and Γ_n, we have <cit.>:

    μ^n()ds^n()= ds^n-1(^n-1()),   |1-μ^n()|≤ C Δ t,   ∈Γ_n,

with some C>0 depending only on surface velocity . Therefore, for Δ t small enough
we obtain

    ∫_Γ_nf_0(c^n-1)  ds^n = ∫_Γ_n-1f_0(c^n-1)  [μ^n]^-1 ds^n-1≤ (1+C Δ t)∫_Γ_n-1f_0(c^n-1)  ds^n-1.

In turn, the surface gradients on Γ_n and Γ_n-1 are related through (cf., e.g., <cit.>)

    ∇_Γ_n c^n-1()= _Γ_n()(1-d()())∇_Γ_n-1 c^n-1(^n-1()),  ∈Γ_n,

where _Γ_n is the orthogonal projector on the tangential space on Γ_n,
d() a signed distance function for Γ_n-1, and ()=∇^2 d() is the shape operator.  Our assumptions on Γ(t) and its smooth evolution imply the uniform in n bounds _L^∞(_δ(Γ_n-1))≤ C and
d_L^∞(Γ_n)≤ CΔ t. Together with (<ref>), this leads to estimate

    c^n-1^2_Γ_n≤ (1+C Δ t)  c^n-1^2_Γ_n-1.


By employing (<ref>) and (<ref>) in (<ref>), we are led to the following bound

    ϵ/2 c^n^2_Γ_n + 1/ϵ∫_Γ_nf_0(c^n)   ds + Δ tμ^n^2_Γ_n     ≤(1+C Δ t) (ϵ/2 c^n-1^2_Γ_n-1+1/ϵ∫_Γ_n-1f_0(c^n-1)  ds) 
             +CΔ t(ϵ+Δ t/ϵ) c^n_Γ^n^2.

We sum the inequalities (<ref>) over n=1,…,N  to get

    ϵ/2 c^N^2_Γ_N + 1/ϵ∫_Γ_Nf_0(c^N)   ds +Δ t∑_n=1^Nμ^n^2_Γ_n     ≤C(1+|Δ t|/ϵ^2) ∑_n=0^NΔ t (  ϵ/2 c^n^2_Γ_n+1/ϵ∫_Γ_nf_0(c^n)  ds)         +ϵ/2 c^0^2_Γ_0+1/ϵ∫_Γ_0f_0(c^0)  ds.

Now assume that Δ t is small enough such that 1-Δ tC(1+|Δ t|/ϵ^2)=α>0
and apply the discrete Gronwall inequality  to obtain

    ϵ/2 c^n^2_Γ_n + 1/ϵ∫_Γ_nf_0(c^n)   ds +Δ t∑_k=1^nμ^k^2_Γ_k≤ C( data),  for all n=1,…,N,

where the constant on the right-hand side depends on given problem data, i.e., on , Γ(0), c^0, but does not depend on n and Δ t. The constant C( data) is also uniformly bounded with respect to ϵ if |Δ t|≲ϵ.




 §.§ Fully discrete method


The (quasi)-stationary system (<ref>)–(<ref>) involves only integrals
over Γ^n and _δ(Γ^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 <cit.>, to discretize problem
(<ref>)–(<ref>) in space.

We start by assuming a shape regular triangulation  _h
of the bulk computational domain Ω⊂^3.
At this point, we consider the same surface-independent mesh for all t∈(0,T).
However, in practice dynamic refinement is possible and the
numerical examples in Section <ref> demonstrate its
feasibility and use. We will address this option later.

The domain formed by all the elements in _h^Γ is denoted by .
On _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^ bulk:

    V_h^ bulk={v∈ C(Ω) :  v∈ P_1(T) for any T∈_h}.



Let the surface be given implicitly as the zero level of the level set function ϕ: Ω×[0,T]→:

    Γ(t)={∈^3 : ϕ(,t)=0},

with |∇ϕ|>0 in (). Note that ϕ can be defined as the solution to the level-set equation <cit.>:

    ϕ_t+^ bulk·∇ϕ=0,

where ^ bulk is a given velocity field in Ω such that ^ bulk·=· on Γ(t).
For the purpose of numerical integration, we approximate Γ with a “discrete” surface Γ_h, which is defined as the zero level set of a P_1 Lagrangian interpolant for ϕ on the one time refined mesh:

    Γ_h={∈Ω : ϕ_h()=0,  ϕ_h:=I_h/2(ϕ())∈ V^ bulk_h/2},   with _h=
    ∇ I^2_h/2(ϕ)/|∇ I^2_h/2(ϕ)|,

where I^2_h/2(ϕ) is a P_2 nodal interpolant of the level set function. See Fig. <ref>.
Note that _h is defined in a neighborhood of Γ_h.
The discrete surface Γ_h satisfies the following geometric approximation properties:

    (Γ,Γ_h)=O(h^2)  and  ^e-_h_L^∞(Γ_h)=O(h).

The analysis of the TraceFEM for scalar elliptic problem on a steady surface <cit.>
and parabolic problems on evolving surfaces <cit.> shows that the geometric consistency error (<ref>) 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 (<ref>) is an open problem that we shall address elsewhere.
When only the initial position of the surface and the bulk advection field ^ bulk are given,
the evolution of Γ can be recovered solving eq. (<ref>) numerically.
Then, the calculated solution ϕ_h defines Γ_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
Γ^n_h such that Γ^n+1_h stays inside this narrow band
and so all discrete quantities at t^n+1 are computable.
This narrow band  (Γ^n_h) is defined as the union  of tetrahedra within distance δ_n from
Γ^n_h:

    (Γ^n_h)=
    ⋃{T : T∈𝒯_h : (T,Γ_h)≤δ_n}.

To ensure

    Γ^n+1_h⊂(Γ^n_h),

 we set the minimum thickness of the extension layer to be

    δ_n := c_δΔ t sup_t∈(t_n,t_n+1)·_h_L^∞(Γ_h)

where c_δ≥ 1 is an O(1) mesh-independent constant.
If ϕ_h is an approximate distance function, then the condition (S,Γ_h)≤δ_n in (<ref>) can be replaced by | ϕ_h^n() |≤δ_n  for any vertex  of T.

We define finite element spaces

    V_h^n={v ∈ C((Γ^n_h)) :  v∈ P_1(T), ∀ T⊂(Γ^n_h)}.

These spaces are the restrictions of the time-independent bulk space V_h^ bulk on all tetrahedra from (Γ^n_h).




Our hybrid finite difference in time / finite element in space method is based on the semi-discrete formulation (<ref>)–(<ref>). It reads:
Given c_h^0 ∈ V_h^0, for n=1,…,N find  c_h^n∈ V_h^n and μ_h^n∈ V_h^n satisfying

    ∫_Γ^n_hρ^n (c_h^n-c_h^n-1/Δ t + ^e ·∇ c^n_h )  v_h   ds +  ∫_Γ^n_h M(c^n-1_h) ∇_Γ_h^nμ^n_h  ∇_Γ_h^n v_h   ds       +  ρ_μ∫_(Γ^n_h)∂μ^n_h/∂_h∂ v_h/∂_h  dx= 0, 
       ∫_Γ^n_h  (μ^n_h - β_s [c_h^n-c_h^n-1]   q_h   ds  -  ϵ∫_Γ^n_h∇_Γ_h^n c^n_h  ∇_Γ_h^n q_h   ds       -ρ_c ∫_(Γ^n_h)∂ c^n_h/∂_h∂ q_h/∂_h  dx = 1/ϵ∫_Γ^n_h f_0'(c^n-1_h)  q_h   ds.

for all v_h, q_h∈ V_h^n. Here  ^e()=(^n()), i.e. the lifted data on Γ^n_h from Γ^n. The first term in (<ref>) is well-defined thanks to condition (<ref>) with the index shifted n→ n-1.
In accordance to the analysis of the scalar advection-diffusion problems <cit.>,
parameters ρ_μ and ρ_c are set to be

    ρ_μ=ρ_c=h^-1.




The role of the  term ρ_c ∫_(Γ^n_h)∂ c^n_h/∂_h∂ q_h/∂_h  dx
in eq. (<ref>) is twofold. First, form

    a( c^n_h,  q_h)=ϵ∫_Γ^n_h∇_Γ_h^n c^n_h  ∇_Γ_h^n q_h   ds+ρ_c ∫_(Γ^n_h)∂ c^n_h/∂_h∂ q_h/∂_h  dx

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  (Γ^n_h).
This can be seen as an implicit extension procedure of finite element solution from Γ^n_h to the narrow band  (Γ^n_h). The same observation holds for μ^n_h
and the corresponding term in eq. (<ref>).
Furthermore, adding two volume terms (one for c^n_h and another for μ^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 <cit.> for unfitted surface finite element methods.

For the realization of the method,
one uses the standard nodal basis functions for the bulk volumetric mesh _h.
At time step n, only the degrees of freedom of the tetrahedra in the narrow band (Γ^n_h)
are active. Since Γ_h^n is the zero level of the P_1 finite element
function ϕ_h^n, the tetrahedra intersected by Γ_h^n are those T∈(Γ^n_h)
for which ϕ_h^n has a change in sign at different vertices.
Then, Γ_h^n∩ 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 ϕ_h^n in  T.
So, it is straightforward to apply standard quadrature rules to compute the surface
integrals in (<ref>)–(<ref>).
We note that on hyperrectangles the implicit representation of Γ_h by the bilinear level-set function
can be treated with a marching cube <cit.> approximation.



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 ϕ_h that implicitly defines Γ_h.
However, it is a non-trivial task to obtain a parametrized representation of the zero level
of ϕ_h (for polynomial degree ≥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  <cit.>. Isoparametric unfitted finite element method is an another elegant and efficient way to simultaneously leverage the geometric and interpolation accuracy of TraceFEM; see  <cit.>.



The finite element method (<ref>)–(<ref>) is different from a
surface unfitted FEM based on the extension of  a PDE from the surface to an ambient bulk
domain <cit.>.
Both methods avoid surface triangulation and remeshing (if the surface evolves).
However, in the  method  originally suggested in <cit.> 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  <cit.>), 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.




§ NUMERICAL EXAMPLES


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

For  the numerical tests, we will consider mobility M as in (<ref>) with σ=1
(the only exception is test 2c, where we set σ=16) and free energy (<ref>).
Furthermore, we always use constant density ρ=1.
Although, there are no reasons to assume ρ= for
motions that violate the inextensibility condition (cf. (<ref>)),
we  expect that moderate smooth variations of ρ would lead to minor changes
in lateral phase dynamics. In addition, ρ=1 is a fair assumption for the
assessment of the numerical method and the model.
We set c_δ=1 in (<ref>) to define the extension strip.
Implementation of the method was been done using the FE package DROPS <cit.>.



 §.§ Validation


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



  §.§.§ 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  and rotated with angular
velocity ω.  We will report the numerical results for two tests:


  - Test 1a: =(1,0,0) and ω=(1,0,0). Its aim is to illustrate
the methodology.

  - Test 1b: =(10,0,0) and ω=(10,0,0). Its aim is to
check the spatial and temporal accuracy.

In both test we take ϵ = 0.1.

We  embed the sphere undergoing rigid body motion in outer domain

    Ω= [-5/3,10/3]×[-5/3,5/3]×[-5/3,5/3].

The initial triangulation _h_0 of Ω 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 ℓ∈ℕ denotes the level of refinement, with the associated mesh size h_ℓ= 2^-ℓ-210/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

    c_x_i= 1/2(1+tanh(δ+x_i/2√(2)ϵ)) .

Here and in the following, = (x_1, x_2, x_3)^T denotes a point in ℝ^3.
For this test we set δ=0 (we will use c_x_i with δ≠0 and i≠2 later).




Fig. <ref> displays the solution computed with mesh ℓ = 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 ϵ.




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

    c^*(x_1,x_2)=1/2(Y(x_1,x_2)+ 1),    Y(x_1,x_2)=x_1 x_2,    t∈[0,0.1].

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

    ċ  -  ρ^-1(M (1/ϵf_0' - ϵΔ_Γ c)) = g.

The exact chemical potential μ^* can be readily computed from eq. (<ref>).
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. <ref> (left) the discrete  L_2(0,T; L^2(Γ_h)) norm:

    c^* - c_h _2,2= (1/T∑_k Δtc^*(t_k) -c_h(t_k)^2_L^2(Γ_h))^1/2

and the discrete  L_2(0,T; H^1(Γ_h)) norm:

    c^* - c_h _2,1= (1/T∑_k Δt(c^*(t_k) -c_h(t_k))^2_L^2(Γ_h))^1/2

of the error for order parameter c (blue lines) and chemical potential μ (red lines)
plotted against the refinement level ℓ. The time step was refined together with the mesh size according to Δ t=4^1-ℓ/10.
All the norms reported in Fig. <ref> (left) are
computed on the approximate surface Γ_h, where c^* and μ^*
were defined through their normal extensions from Γ.
The observed second order convergence in the L_2(0,T; L^2(Γ_h)) for both concentration and chemical potential as well as the first order convergence in the L_2(0,T; H^1(Γ_h)) norm for the concentration are consistent with the well known error analysis of finite element methods for the Cahn–Hilliard equation <cit.>.  For the error in chemical potential we also observe the almost second order of convergence in L_2(0,T; H^1(Γ_h)) norm. We do not have an explanation for this apparent super-convergence.




  §.§.§ Oscillating ellipsoid


The second series of tests is inspired by a numerical experiment in <cit.>.
We consider time-dependent surface Γ(t) = {∈ℝ^3: ϕ(, t) = 0 } with

    ϕ(, t) = (x_1/a(t))^2+x_2^2+x_3^2 - 1 ,   a(t)=1+0.2sin(2πkt),

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 Ω = [-5/3,5/3]^3.
As initial solution to the Cahn–Hilliard equation, we take a small perturbation
about 0.5 given by:

    c_0=0.5+0.05cos(2πx)cos(2πy)cos(2πz).

We study the evolution of the discrete Lyapunov energy:

    E^L_h(c_h) = ∫_Γ_h f(c_h) ds = ∫_Γ_h( 1/ϵf_0(c_h) + 1/2ϵ |  c_h |^2 )ds

for three datasets:


  - Test 2a: k=1, σ= 1, ϵ=0.1, Δt=0.01.

  - Test 2b: k=1, σ=1, ϵ=0.01, Δt=0.01.

  -  Test 2c: k=5, σ=16, ϵ=0.1, Δ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
Ω into sub-cubes that are further divided into tetrahedra.
We consider different levels of refinement ℓ associated with
mesh size h=2^-ℓ-210/3.



Test 2a. We present the evolution of the discrete Lyapunov energy in Fig. <ref> (right).
We see that mesh ℓ = 1 is too coarse to correctly capture the evolution of E^L_h(c_h).
The results obtained with mesh ℓ = 2 are very close to results obtained with meshes ℓ = 3,4,
which are almost superimposed for the entire time interval under consideration.
Fig. <ref> (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 ϵ
in the next test. We recall that ϵ
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 ϵ is numerically challenging.






Test 2b.
We keep all the model and discretization parameters as in test 2a, but we reduce the value
of ϵ to 0.01.
Fig. <ref> (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. <ref> (left) we see that the results computed with mesh ℓ = 4
are still quite far from the results obtained with mesh ℓ = 5, which are almost
superimposed to the results for mesh ℓ = 6. This indicates that both meshes
ℓ = 5 and ℓ = 6 are sufficiently refined for  ϵ = 0.01.
Notice that the discrete Lyapunov energy computed with meshes ℓ = 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 ℓ = 6
in Fig. <ref>. 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. <ref> (right).
We see that the results computed with meshes ℓ = 4, 5 are in excellent
agreement. 
In <cit.>, 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 <cit.>, 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 ℓ = 4, 5 decreases (not monotonically)
till around t = 0.35. For t > 0.35, Fig. <ref> (right) suggests that the solution converges to a time
periodic solution.


































We show the evolution of the solution computed with mesh ℓ = 6
in Fig. <ref>. Note that we use the truncated scale here for c to better illustrate
the solution, which is otherwise rather diffusive (since ϵ=0.1 is used).
From Fig. <ref>, 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 <cit.>.










 §.§ Colliding spheres

We consider an evolving surface Γ(t) that undergoes a topological change and experiences a
local singularity. The dynamics of Γ(t) is inspired by tests presented in <cit.>.
The computational domain is Ω = [-10/3,10/3]×[-5/3,5/3]^2 and
the evolving surface is the zero level set of level set function

    ϕ(, t) = 1 - 1/||  - _c^+(t)||^3 - 1/||  - _c^-(t)||^3,

with

    _c^±(t) = ± (3/2 - w  t, 0, 0 ),

for t ∈ [0, 1.5/|w|].
The initial configuration Γ(0) is close to two balls of radius 1, centered at
_c^±(0). Parameter w is the collision speed.
For t > 0, w>0 the two spheres approach each other until
time t̃≈ 0.235/w, when they touch at the origin; then, for t ∈ (t̃,  1.5/w), surface Γ(t) is simply connected.

In the vicinity of Γ(t), the gradient ∇ϕ and the time derivative ∂_t ϕ
are well defined and given by simple algebraic expressions.
The normal velocity field of Γ(t) can be computed to be

    = -  ∂_t ϕ/| ∇ϕ |^2∇ϕ.

We assume the following material velocity:

    =(1-tanh( |x_1 |/0.1))+ tanh( |x_1 |/0.1)( -sgn(x_1) w, 0, 0),

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 Ω with mesh size h=2^-ℓ-210/3, where ℓ is the refinement level.

The setup is a simple model of two plasma membrane fusion. In all further experiments we take ϵ = 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 μm, on which phase separation can be visualized using fluorescence microscopy
<cit.> with a resolution of about 300 nm, the thickness of transition region between the
phases is approximately 1% of vesicle diameter.






  §.§.§ 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 δ = 0 for the ball initially centered at _c^-(0)
and c_x_1  with δ = 0



for the ball initially centered at _c^+(0). See (<ref>)
for the definition of c_x_i.
We set Δ t = 0.001.

Fig. <ref> displays the evolution of the surface and the solution
for t ∈ [0.05,1]. Before the
two balls touch at t = t̃, we see that the ball centered at _c^-(t) (resp., _c^+(t)) has a horizontal
(resp., vertical) interface separating a pink and black domain of equal size. For t >  t̃,
the black and pink domain of the ball centered at _c^-(t) get in contact with the black domain
of the ball centered at _c^+(t). Between the pink domain of the ball centered at _c^-(t) and the black domain
of the ball centered at _c^+(t) an interface gets formed. Such interface is virtually located on the curve of minimal length
on the simply connected surface Γ(t), for t >  t̃.
This is consistent with the well-known limiting (as ϵ→0) behavior of
the stationary  Cahn-Hilliard problem that corresponds to solving the isoperimetric problem <cit.>.





  §.§.§ Pattern formation on colliding sphere














Let t_s=0.23/w < 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, Δ t = 10^-4 for t ∈ (t_s, t_s + 0.01] and Δ t = 10^-2 for t ∈ (t_s+0.01, t_s + 1].

  - Fast collision: w = 10, T = 0.1+t_s, Δ t = 10^-4 for t ∈ (t_s, t_s+ 0.01] and Δ t = 10^-3 for t ∈ (t_s+0.01, t_s+0.1].

The two initial spheres touch each other at the origin at time t̃≈0.235 for the slow collision test and t̃≈0.0235 for the fast collision test. 
Let rand() be a uniformly distributed random number between 0 and 1.
For both tests, we take c_0() =rand() for the ball initially centered at _c^-(0)
and c_0 = 0 for the ball initially centered at _c^+(0).

As we have observed in our previous work <cit.>
and is well known (see, e.g.,<cit.>),
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. <ref> shows the  evolution of the surface and the solution
for the slow collision test for t ∈ [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 t̃≈0.235, the already phase separated solution of the ball centered at _c^-(t)
enters in contact with the black domain that covers the entire ball centered at _c^+(t).
Similarly to what observed in Sec. <ref>, for t >  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 Γ(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 Γ(t) are:


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

  * The contact with the black domain covering the ball centered at _c^+(t).

  * 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 Γ(t_s), hereafter called single ball, and surface Γ(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() =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() =rand() for 
with x_1 ≤ 0 and c_0 = 0 for  with x_1>0.
We consider the same meshes described in Sec. <ref>.

In Fig. <ref>, we see the results for both simulations for t ∈ [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. <ref> the evolution of the surface and the solution for t ∈ [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 _c^-(t).
Upon phase separation, we see the same elongated pink domains
aligned with the x_1-axis already observed in Fig. <ref>
and <ref>.
The notable difference between the corresponding solutions
in Fig. <ref> and <ref>
(taking into account the 10 times faster dynamics in the latter)
are the thicker interfaces along the bottleneck in Fig. <ref>,
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. <ref>,
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. <ref> 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.






 §.§ Splitting spheres







We consider an evolving surface undergoing a reverse dynamics with respect to the
one considered in Sec. <ref>: the initial configuration is a sphere centered at the origin,
which splits into two droplets, in turn evolving towards spheres centered at
_c^±(T). The computational domain is again Ω = [-10/3,10/3]×[-5/3,5/3]^2.
The evolving surface Γ(t) is the zero level set of ϕ(, 1.5/|w|+t) with w=-1, where ϕ is given by (<ref>). The velocity vector field is defined by (<ref>) with w=-1.
Surface Γ(t) is simply connected for t < 1.5-t̃, where
t̃ = 0.235.
Obviously, this evolving surface experiences a local singularity too.


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


  - Rotated interface: with initial solution is given by c_0^x_3 from (<ref>), with δ=0 and rotated by π/4 in x_1x_2 plane.


  - Horizontal interface: with initial solution c_0^x_3 from (<ref>), with δ=20.



Fig. <ref> and <ref> show the evolution
of the interface and the solution for t ∈ [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. <ref> 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. <ref>, bottom right panel.




§ CONCLUSIONS


We presented a formulation of the Cahn–Hilliard equation
on a time-dependent surface Γ(t) that uses tangential calculus
induced by the embedding of Γ(t) in ℝ^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.

siam


