{'fraction_non_alphanumeric': 0.0615807058642805, 'fraction_numerical': 0.021218172740888667, 'mean_word_length': 4.023051697530864, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Error estimates for general non-linear Cahn–Hilliard equations on evolving surfaces

In this paper, we consider the Cahn–Hilliard equation on evolving surfaces with prescribed velocity and a general non-linear potential. High-order evolving surface finite elements are used to discretise the weak equation system in space, and a modified matrix–vector formulation for the semi-discrete problem is derived. The anti-symmetric structure of the weak equation system is preserved by the matrix–vector formulation and it is utilised to prove optimal-order and uniform-in-time error estimates. An extension of the convergence results is given for general non-linear Cahn–Hilliard equations on evolving surfaces. The paper is concluded by a variety of numerical experiments.


Introduction
This paper studies the Cahn-Hilliard equation on evolving surfaces with prescribed surface velocity. The non-linear potential is only assumed to satisfy locally Lipschitz-type assumptions. The Cahn-Hilliard equation is formulated as a system of second order equations, exhibiting an anti-symmetric structure. The semi-discretisation of the system by higher order evolving surface finite elements, cf. [10,18], preserves this anti-symmetric structure, which is utilised to prove a new convergence result. Optimal-order uniform-in-time error estimates in the L 2 and H 1 norms for both solution variables are proven. Cahn and Hilliard first described an equation modelling phase separation processes in [5]. Which since found many applications, e.g. in foam modelling, dendritic flow, phase separation models, and image processing.
The Cahn-Hilliard equation on a stationary surface with boundary was first investigated by Du, Ju and Tian in [9]. They study a full discretisation of the Cahn-Hilliard equation with homogeneous Dirichlet boundary conditions, and prove optimal-order error estimates in the L 2 norm for u, using linear finite elements. Elliott and Ranner were the first to consider the Cahn-Hilliard equation on a closed evolving surface with a prescribed velocity in [15]. They proved optimal-order uniform-in-time error estimates in the L 2 and H 1 norms in the first variable and optimal-order L 2 -in-time error estimates in the L 2 and H 1 norms in the second variable, using a discretisation by linear evolving surface finite elements. Using a new stability proof, the results of this paper improve the error estimates for the second variable from optimal-order L 2in-time to optimal-order uniform-in-time estimates. Furthermore, the proof is built, such that it is extendable to general non-linear Cahn-Hilliard equations in a straightforward way.
The main result of this paper, optimal-order uniform-in-time error estimates, is proven by studying the stability and consistency of the semidiscretisation.
In the stability analysis, the difference between the Ritz map of the exact solution and the numerical solution is estimated in terms of defects and their time derivatives. To account for initial errors in the second variable, a modification of the second-order system is required. The stability proof uses energy estimates, performed in the matrix-vector formulation, and utilises the anti-symmetric structure of the error equations, testing the error equations with the errors and also with their time derivatives. The stability analysis was first developed for Willmore flow in [20]. A uniform-in-time L ∞ bound for the numerical solution is key to estimate the non-linear term. It is obtained from the time-uniform H 1 norm error bounds using an inverse estimate and exists for a small time due to a continuous initial function. The stability proof is independent of geometric errors.
We also give details to an extension of the stability proof to rather general non-linear Cahn-Hilliard equations, which include proliferation terms [27, equation (3.1)], where an additional non-linearity appears in both variables in both equations of the system. The generality of the stability proof can also be seen through the related results in [20] and [17]. A further advantage of this stability proof, is that it is expected to be generalisable to a full discrete method based on linearly implicit backward difference methods, on which we intend to report in a subsequent work.
In the consistency analysis the L 2 norms of the defects and their time derivatives are estimated. The bounds use geometric error estimates, including interpolation and Ritz map error estimates, bounds on the discrete surface velocity, and geometric approximation errors for high-order evolving surface finite elements, see [18].
The paper is structured as follows. In Section 2, based on the works [10] and [15], the weak formulation for the Cahn-Hilliard equation on evolving surfaces is derived as a system of equations. In Section 3 the evolving surface finite element method is used to discretise this system of equations in space. The obtained semi-discrete problem is written as a matrix-vector formulation. In Section 4 the novel error estimates proved in this work are stated and discussed in comparison to the existing results by Elliott and Ranner [15]. Section 5 contains the stability part of the proof. Section 6 treats the consistency part of the proof. In Section 7 the two parts are combined to prove the main result. In Section 8 the result is extended to general non-linear Cahn-Hilliard equations. In Section 9 a full discretisation to the problem is given, cf. [2,1]. In Section 10 the theoretical results are complemented by numerical experiments.

Cahn-Hilliard equation on evolving surfaces
In the following we consider a smoothly evolving closed hypersurface Γ (t) ⊂ R d+1 , with d ≤ 3, for 0 ≤ t ≤ T . The initial surface Γ (0) = Γ 0 is given (and at least C 2 ), and it evolves with the given and sufficiently smooth velocity v. The surface Γ (t) is given as the image of a smooth mapping X : Γ 0 ×[0, T ] → R d+1 , by Γ (t) = {X(p, t) | p ∈ Γ 0 }. The embedding X and the velocity v satisfy the ordinary differential equation (ODE): Let ν denote the unit outward normal vector to Γ (t). Then the surface (or tangential) gradient on Γ (t), of a function u : Γ (t) → R, is denoted by ∇ Γ (t) u, and is given by ∇ Γ (t) u = ∇ū − (∇ū · ν)ν (the surface gradient is independent of the extensionū), while the Laplace-Beltrami operator on Γ (t) is given by ∆ Γ (t) u = ∇ Γ (t) · ∇ Γ (t) u. Moreover, ∂ • u denotes the material derivative of u, i.e. ∂ • u(·, t) = d/dt(u(X(·, t), t)) = ∂ t u(·, t) + v · ∇u(·, t). The space-time manifold will be denoted by G T = ∪ t∈[0,T ] Γ (t) × {t}. For more details on these notions we refer to [10,11]. In this paper we consider the Cahn-Hilliard equation on evolving surfaces. It is the fourth-order equation for a scalar function u : G T → R given by on Γ (t), (2.2) with continuous (and sufficiently regular) initial conditions u(·, 0) = u 0 on the initial surface Γ 0 . The parameter ε > 0 describes the length of the transition regions, we are however not interested in taking the limit ε → 0. The non-linear scalar function W : R → R is a chemical potential and here it is only assumed to have locally Lipschitz continuous derivatives W and W . A typical example is a double well potential, i.e. W (u) = (u 2 − 1) 2 . The solution u ∈ [−1, 1] models the concentration of surfactant fluids, with u = ±1 indicating the pure occurrences of each, cf. [5]. The Cahn-Hilliard equation on a stationary surface Γ can be derived as the H −1 (Γ )-gradient flow of the Ginzburg-Landau energy cf. [15,Remark 2.1]. In [15] it is stated, that to obtain a gradient flow on an evolving surface, a model for the surface velocity v is needed, leading to a coupled system for u and v. In the evolving surface case, w = −ε∆ Γ (t) u + ε −1 W (u) is the variation of the evolving surface Ginzburg-Landau energy, see [29].

Second order system and weak formulation
By introducing the auxiliary function w : G T → R we rewrite the Cahn-Hilliard equation (2.2) into a system of second order partial differential equations: For u, w : with initial data u(0, ·) = u 0 on Γ 0 . On the evolving surface Γ (t) we recall the definition of standard Sobolev spaces L 2 (Γ (t)), and H 1 (Γ (t)) and its higher order variants, endowed with their usual norms, see [10,12].
Assuming that the initial condition u 0 ∈ H 2 (Γ 0 ) has bounded energy (2.3) then, there exists a unique weak solution pair (u, w) to the Cahn-Hilliard equation (2.5), whose energy stays bounded by the energy of the initial data. Furthermore, u and w have H 2 (Γ (t)) regularity, and satisfy the estimate, for

Abstract formulation
We will use the time-dependent bilinear forms, cf. [12,11], for any u, ϕ ∈ H 1 (Γ (t)): We further define a * (t; ·, ·) = a(t; ·, ·)+m(t; ·, ·). All bilinear forms are symmetric in u and ϕ, m and a * are positive definite, while a is positive semi-definite. Whenever it is possible, without confusion, we will omit the omnipresent timedependence of the bilinear forms and write m(·, ·) instead of m(t; ·, ·). We note here that the bilinear forms directly generate the (semi-)norms: The weak formulation (2.5) can then be rewritten, using the bilinear forms from above, as

and (2.6) is rewritten as
The transport formula for the above bilinear forms, [12,Remark 3.3], will be used later on, and reads (2.10)

Semi-discretisation of Cahn-Hilliard equation on evolving surfaces
For the numerical solution of the above examples we consider a high-order evolving surface finite element method. In the following, from [10,11,7,18], we will briefly recall the construction of the discrete evolving surface, the highorder evolving surface finite element space, the lift operation, and the discrete bilinear forms, etc., which will be used to discretise the Cahn-Hilliard equation of Section 2.

Evolving surface finite elements
The smooth initial surface Γ (0) is approximated by a k-order interpolating discrete surface denoted by Γ h (0) := Γ k h (0), with vertices p j , j = 1, . . . , N , and is given by the (high-order) triangulation, with maximal mesh width h. More details and the properties of such a discrete initial surface can be found in Section 2 of [7].
The triangulation of the surface Γ (t), denoted by Γ h (t) := Γ k h (t), is obtained by integrating the ODE (2.1) (with the known velocity v) from time 0 to t for all the nodes p j of the initial triangulation. The nodes x j (t) are on the exact surface Γ (t) for all times. The discrete surface Γ h (t) remains to be an interpolation of Γ (t) for all times. We always assume that the evolving (high-order) triangles are forming an admissible triangulation of the surface Γ (t), which includes quasi-uniformity, and that the discrete surface is not a global double covering, cf. Section 5.1 of [10]. The discrete tangential gradient on the discrete surface Γ h (t) is given by: understood in an elementwise sense, with ν h denoting the normal to Γ h (t).
The high-order evolving surface finite element space S h (t) H 1 (Γ (t)) on Γ h (t) is spanned by continuous, piecewise linear nodal basis functions on Γ h (t) satisfying for each node (x j ) N j=1 φ k (x j , t) = δ kj , for j, k = 1, . . . , N.
Then the finite element space is given as The discrete velocity V h of the surface Γ h (t) is the evolving surface finite element interpolation of the surface velocity v of Γ (t), i.e. (3.1) The discrete material derivative is then given by The key transport property of basis functions derived in Proposition 5.4 in [10], is

Lift
Following [10,7], we define the lift operator · to compare functions on Γ h (t) with functions on Γ (t). For functions ϕ h : Γ h (t) → R, we define the lift as Finally, the lifted finite element space is denoted by S h (t), and is given as

Discrete bilinear forms
The time-dependent discrete bilinear forms on S h (t), i.e. the discrete counterparts of m, a and g, are given, for u h , ϕ h ∈ S h (t), by (3.5) As in the continuous case we let a * h (t; ·, ·) = a h (t; ·, ·) + m h (t; ·, ·). The discrete bilinear forms, clearly inherit the properties of their continuous counterparts, such as the transport formula (2.10), see, e.g. [12,18].
As in the continuous case, the discrete bilinear forms directly generate the discrete (semi-)norms, for u h ∈ S h (t), According to [10,7], the discrete norms and their continuous counterparts are h-uniformly equivalent, for ϕ h ∈ S h (t), (3.6)

Semi-discrete problem
The semi-discrete problem corresponding to the Cahn-Hilliard equation (2.5) reads: Find a solution u h (·, t) ∈ S h (t) with time-continuous discrete material with given initial data u h (·, 0) = u 0 h on Γ 0 h , while w h (·, 0) = w 0 h is obtained by solving the elliptic problem (3.7b) (with a known u 0 h ) once. Equivalently, the semi-discrete problem corresponding to the weak form (2.6), using the discrete version of the transport formula (2.10) for (3.7a), reads: Find a solution u h (·, t) ∈ S h (t) with time-continuous discrete material derivative ∂ • h u h (·, t) ∈ S h (t) and w h (·, t) ∈ S h (t) such that for all test functions again, with given initial data u h (·, 0) = u 0 h on Γ 0 h .

Matrix-vector formulation
We collect the nodal values of , the solution pair of the semi-discrete problem (3.7), into the vectors u(t) = (u 1 (t), . . . , u N (t)) ∈ R N and w(t) = (w 1 (t), . . . , w N (t)) ∈ R N . We define the time-dependent matrices, the mass and stiffness matrix, corresponding to the bilinear forms m h and a h , respectively, and the semi-linear term involving W : (3.9) We further define the matrix corresponding to the bilinear form a * h : We also note that, via the transport property (3.3), the time derivative of the mass matrix is given bẏ The discrete material derivative of a surface finite element function u h (·, t) ∈ S h (t), with nodal values u(t), again by using the transport property (3.3) of the basis functions and the product rule, is given by Thus, the nodal values of ∂ • h u h are given by the vectoru(t). The finite element semi-discretisation of the Cahn-Hilliard equation (3.7) then reads: The anti-symmetric structure of (3.10), which is shared with (2.4) and (3.7), can be best seen in the rewritten form: In order to exploit this favourable structure, the stability analysis will use the matrix-vector system (3.10).
For computations, it is however more advantageous to use the equivalent matrix-vector formulation where the surface velocity V h does not appear directly, as compared to the term withṀ(t) in (3.10).

A modified problem
The initial value u(0) can be chosen suitably, on the other hand the initial value w(0) is obtained from the second equation of the system (3.10) (or equivalently (3.11)). Our error analysis will require that the errors in both initial values to be O(h k+1 ) in the H 1 (Γ h ) norm. For u this can be achieved using the Ritz map (in which case the initial error in u vanishes), however, such an error estimate is still not feasible for w.
To obtain optimal-order error estimates we modify the equation (3.10b) (and equivalently (3.11b) as well) using a time-independent correction term. Letw(0) ∈ R N denote the solution obtained from (3.10b) at time t = 0, and let w * (0) ∈ R N contain the nodal value of the Ritz map of w(0), and set The second equation is then modified, such that the system (3.10) reads: Similarly, the equivalent system (3.11) is modified to: The semi-discrete problems (3.7) and (3.8) are modified accordingly: The initial value w(0) is obtained by solving the elliptic problem (3.13b) at t = 0, which, via (3.12) and (3.10b), yields  [15] In [15] Elliott and Ranner proved the first error estimates for the evolving surface finite element discretisation of the Cahn-Hilliard equation on evolving surfaces. In particular, they proved optimal-order error estimates in the L ∞ (L 2 ) and L ∞ (H 1 ) norms for u, and in the L 2 (L 2 ) and L 2 (H 1 ) norms for w, using linear surface finite elements. , and assume that they are sufficiently regular, e.g. satisfy (4.1).

Error estimates by Elliott and Ranner
Then there exists an h 0 > 0 such that for all h ≤ h 0 the errors between the solutions u and w and the evolving surface finite element solutions u h and w h of degree 1 of (3.7), satisfy the optimal-order error estimates in both variables: where the constant C > 0 is independent of h, but depends on ε, on the bounds of the Sobolev norms of the solution u, w, and ∂ • u, on the surface evolution, and on the length T of the time interval.
The constant C is explicitly given in the cited theorem.
The analysis presented in [15] is not restricted to linear evolving surface finite elements. By a straightforward combination of the proofs used for [15, Theorem 5.1] with the relevant high-order approximation results shown in [18], the above theorem stays valid for degree k evolving surface finite elements with h k+1 and h k in place of h 2 and h. For a detailed proof we refer to [3, Section 4.1].

Uniform-in-time error estimates
We next state a new convergence result for the evolving surface finite element semi-discretisation of order k ≥ 1. Then there exists an h 0 > 0 such that for all h ≤ h 0 the errors between the solutions u and w and the evolving surface finite element solutions u h and w h of degree k, with nodal vectors solving (3.14), satisfy the optimal-order uniform-in-time error estimates in both variables, for 0 ≤ t ≤ T , whereas the time derivatives of the error in u satisfies, for where the constant C > 0 is independent of h and t, but depends on ε, on the bounds of the Sobolev norms of the solution u and w, on the surface evolution, and on the length of the time interval T .
Sufficient regularity conditions on u and w required by Theorem 4.2 are: and for the surface velocity: 2) The main differences between Theorem 4.1 by Elliott and Ranner [15] and Theorem 4.2 from above are the following: The former result provides a better scaling in ε between the L 2 norm and H 1 semi-norm than the latter, note, however, that in both theorems the constant C (unfavourably) depends on ε. Our result proves uniform-in-time error estimates in the H 1 norm not only for the error in u but also for the error in w, which is not available through the analysis of Elliott and Ranner [15]. Theorem 4.2 also provides error estimates in the errors of the material derivatives for u. Naturally, the additional error estimates require more regularity assumptions on u, w and v.
Theorem 4.2 will be proven by studying the questions of stability and consistency. The consistency of the algorithm is shown by proving high-order estimates for the defects (the error obtained by inserting the Ritz map of the exact solutions into the method), which are obtained by using geometric and approximation error estimates for high-order evolving surface finite elements from [18], which combines techniques of [10,12] and [7].
The main issue in the proof is stability, i.e. a mesh independent, uniformin-time bound of the errors in terms of the defects. The main idea of the stability proof was originally developed for Willmore flow [20], and it relies on energy estimates that exploit the anti-symmetric structure of the Cahn-Hilliard equation, see (2.4), (3.7), and (3.13). The basic idea of the stability proof is concisely sketched in Figure 5.1. In order to estimate the non-linear term, a key issue in the stability proof is to ensure that the L ∞ norm of the error in u remains bounded. The uniform-in-time H 1 norm error bounds together with an inverse estimate provide a bound in the L ∞ norm. Similarly, it is also possible to show such an L ∞ norm bound for the error in w, provided by our uniform-in-time H 1 norm bounds in both u and w. This enables to generalise our stability proof, and hence also the convergence proof, to quite general non-linear problems, which contain a non-linear expression in both equations in (2.4), possibly involving the gradient of u as well. The generalisations only require very minor changes in the stability proof, as elaborated in Section 8. Such examples include Cahn-Hilliard equation with a proliferation term [27]. More details to these problems and to their stability analysis are given in Section 8, along with Theorem 8.1 which generalises Theorem 4.2 to these more general non-linear Cahn-Hilliard equations.

Preliminaries
This section is dedicated to the definition of a few concepts, such as the comparison of various quantities on different discrete surfaces and a generalised Ritz map, which will be all used throughout the stability analysis.
The finite element matrices M(t), A(t), and K(t) induce (semi-)norms which correspond to discrete Sobolev (semi-)norms: for a vector w ∈ R N corresponding to the finite element function w h ∈ S h (t). From [21, Lemma 4.6] we recall the following estimates for the time derivatives of the mass and stiffness matrix, and, additionally, we prove that they also hold for the second order time derivatives.
where the constant c > 0 is independent of h, but depends on the surface velocity v.
Proof The first two estimates were shown in Lemma 4.6 of [21]. We prove the estimate (5.2c) for the second derivative of the mass matrix.
To estimate the first integral we recall how to interchange surface differential operators with the material derivative [13, Lemma 2.6] , for u : and analogously for discrete differential operators. Then the second formula from above is used to estimate the integral, together with the bounds on the discrete velocity V h . The latter is implied by the sufficient regularity of the velocity v, and recalling that V h is the interpolation of v, cf. (3.1), see Lemma 6.2 or [3, Lemma 3.1.6]. We altogether obtain The second integral is directly bounded by The estimate for the stiffness matrix is shown by the exact same arguments, cf. [3].

Error equations and defects
Before turning to the stability analysis, let us define a Ritz map of the exact solution onto the evolving surface finite element space, from [25,18] we recall the definition of a time-dependent Ritz map on evolving surfaces: R h : The Ritz map is then defined as the lift of We will often suppress the omnipresent time-dependency of the Ritz map. In [25] it was shown that the above Ritz map is well-defined, optimalorder error estimates for high-order evolving surface FEM were shown in [18], and are recalled in Lemma 6.4. Let us consider now the (unlifted) Ritz map of the exact solutions u and w of (2.4), which are denoted by whose nodal values are collected into the vectors The nodal vectors of the Ritz maps of the exact solutions satisfy the system (3.10) only up to some defects d u (t) and : The errors between the nodal values of the semi-discrete solutions and of the Ritz maps of the exact solutions are denoted by e u (t) = u(t) − u * (t) and e w (t) = w(t) − w * (t) in R N . By subtracting (5.5) from (3.13) we obtain that the errors e u and e w (corresponding to the functions e u h and e w h ∈ S h (t)) satisfy the following error equations: Then there exists h 0 > 0 such that the following error bound holds for h ≤ h 0 and 0 ≤ t ≤ T : The constant C > 0 is independent of t and h, but depends on 0 < < 1, and exponentially on the final time T .
Proof The proof is based on energy estimates, and its basic idea is very similar to that of [20]. Proving uniform-in-time H 1 norm error estimates is essential for handling the non-linear term, which is done by deriving an L ∞ norm bound for the errors using an inverse estimate.
In order to achieve a uniform-in-time stability bound, two sets of energy estimates are needed. These energy estimates strongly exploit the anti-symmetric structure of (2.4). (i) In the first, an energy estimate is proved for e u , but comes with a critical term involvingė u . (ii) The second estimate uses the time derivative of (5.6b), leads to a bound of this critical term and also to a uniform-in-time bound for e w . The combination of these two energy estimates will give the above stability bound. The structure and basic idea of the proof is sketched in Figure 5.1. In order to handle the semi-linear term we first Stability estimate Proposition 5.1 combining (i) and (ii) Sketch of the structure of the energy estimates for the stability proof. In the diagram r 1 and r 2 denote the right-hand sides of (5.6a) and (5.6b). (Note that, after time differentiation, the term R 2 not only contains the time derivative of r 2 , but other terms involving derivatives of matrices as well.) prove the stability bound on a time interval where the L ∞ norm of e u h is small enough, and then show that this time interval can be enlarged up to T .
In the following c and C are generic constants that take different values on different occurrences. Whenever it is possible, without confusion, we omit the argument t of time-dependent vectors but not of time-dependent matrices. By > 0 we will denote a small number, used in Young's inequality, and hence we will often incorporate h independent multiplicative constants into those, yet unchosen, factors.
We start by stating that there exists a maximal time 0 < t * ≤ T such that, for all t ≤ t * , Since e u h (·, 0) = 0 and since u h and u * h are continuous in time, we directly infer that t * > 0.
Thus, by the assumption that the Ritz map of the exact solution satisfies u * h (t) L ∞ (Γ h (t)) ≤ M , we obtain the following bound for the numerical solution: for all 0 ≤ t ≤ t * and for h ≤ h 0 sufficiently small. Thus, for f ∈ C(R) for all 0 ≤ t ≤ t * and h ≤ h 0 sufficiently small. We first prove the stated stability bound for 0 ≤ t ≤ t * , and then show that indeed t * coincides with T .
Energy estimate (i): We take the first error equation (5.6a) and test it with e u , while the second one (5.6b) is tested by ε −1 e w , to obtain By adding the two equations, and by the symmetry of A, we eliminate the mixed term e T u A(t)e w , and obtain Using the product rule and symmetry of M we rewrite the first term as which altogether yields Similarly, we test (5.6a) by e w and (5.6b) byė u , now a subtraction leads to cancelling the mixed term e T w M(t)ė u , and again by the product rule and the symmetry of A, we obtain Dividing the second equality by ε and taking their linear combination yields 13) The terms on the right-hand side are now estimated separately. The terms involving time derivatives of matrices are estimated using Lemma 5.1, by For the non-linear terms, using (5.11) and the local-Lipschitz property of W , we obtain 15) and similarly The defect terms are estimated by Cauchy-Schwarz inequality, as The terms involving the correction term are bounded similarly as the defect terms. Using equality (5.7) and the norm equivalence in time [14,Lemma 4.1] (to change the time from 0 to t), we obtain Altogether, by the combination of the estimates (5.14)-(5.17) with (5.13), by multiple Young's inequalities (with > 0 chosen later on) and by absorptions to the left-hand side, we obtain Integrating from 0 to t ∈ (0, t * ], and using that e u (0) = 0, we obtain the first energy estimate: (5.20) Note the critical term, with ė u (s) K(s) , on the right-hand side, which cannot be bounded or absorbed in any direct way.
Energy estimates (ii) To control the critical term on the right-hand side of (5.20) we will now derive an energy estimate, which includes this term on the left-hand side. To this end, we first differentiate the second equation of (5.6) with respect to time (note that the time-independent ϑ vanishes), and, after rearranging the terms, we obtain the following system: Testing the error equation system (5.21) twice, similarly as before in Part (i), would not lead to a feasible energy estimate, but to a bound which include a new critical termė u . The issue is avoided by separating the two estimates for the error equations, (ii.a) and (ii.b), and then taking their weighted combination in (ii.c).
(ii.a) We test (5.21a) byė u and (5.21b) by ε −1 e w , adding the two equations together to cancel the mixed termė T u A(t)e w , and using the product rule as before, we obtain (5.22) The right-hand side terms are again estimated separately. The ones in the first line are bounded, using Lemma 5.1, by The non-linear term occurs differentiated with respect to time. Therefore, with the help of the transport formula (2.10) and inserting ∓ Γ h (t) W (u * h (·, t)) ∂ • h u * h (·, t)φ j after the second equality, we obtain 24) where we used the local-Lipschitz continuity of W and W , together with (5.12), and the bounds on V h obtained by interpolation error estimates (for details, see [3, Lemma 3.1.6]).
The defect terms are bounded, similarly as before, by (5.25) Altogether, by plugging in (5.23)-(5.25) into (5.22), then using Young's inequalities (possibly with a small number > 0) and absorptions to the lefthand side, we obtain the first energy estimate of this part: (5.26) (ii.b) We now test (5.21a) byė T w and(5.21b) byė T u , then subtracting the second from the first equation to cancel the mixed termė T w M(t)ė u , then using the product rule again, after a division by ε we obtain 27) The terms are again estimated separately. Let us first highlight that it is not possible to directly estimate the terms containingė w (t) in their current form, because there is no term on the left-hand side to absorb them. Therefore, we first rewrite them using the product rule, and estimate them using Lemma 5.1, to obtaiṅ (5.28) The terms with matrix derivatives are bounded, using Lemma 5.1, by

(5.29)
The non-linear term is bounded, similarly to (5.24), by (5.30) The defect terms are bounded, similarly as before, by (5.31) Altogether, by plugging in (5.28)-(5.31) into (5.27), then using Young's inequalities (possibly with a small number > 0) and absorptions to the lefthand side, we obtain the second energy estimate of this part: 32) with a particular constant c 0 > 0 (independent of h, but depending on W , viz. on the constant in (5.24)).
(ii.c) We now take the weighted combination of the energy estimates from (ii.a) and (ii.b): multiplying the estimate (5.26) by 4c 0 and adding it to the estimate (5.32). Collecting the terms, (and since > 0 is yet to be chosen, further constants are incorporated into it), we obtain (5.33) The first term on the right-hand side of (5.33) can directly be absorbed to the left-hand side, which yields Integrating the above inequality from 0 to t ≤ t * , and then dividing by min{c 0 , 1/2}, yields We estimate the newly obtained non-integrated terms on the right-hand side using Lemma 5.1, Cauchy-Schwarz and Young's inequalities, a further absorption, and using that e u (0) and e w (0) are zero, we then obtain This energy estimate now contains the (previously) critical termė u on the left-hand side.
Combining the energy estimates: We now take again a c 1 -weighted linear combination (in order to absorb the term c 1 e u 2 M(t) ) of the two energy estimates (5.20) and (5.34), to obtain By choosing small enough, the first term (previously the critical term) on the left-hand side can now be absorbed. This enables us to use Gronwall's inequality, which then yields the stated stability estimate on [0, t * ]. Now, it only remains to show that, in fact, t * = T , for h sufficiently small. The proven stability bound (for 0 ≤ t ≤ t * ) together with the assumed defect bounds (5.8) imply By an inverse estimate, see, e.g. [4,Theorem 4.5.11], we have, for 0 ≤ t ≤ t * , for sufficiently small h. Therefore the bound (5.10) can be extended to t ≥ t * , which contradicts the maximality of t * unless we already have t * = T . We hence proved the stability bound (5.9) over [0, T ], and completed the proof.

Consistency
Before we turn to proving consistency of the spatial semi-discretisation, and to the proof of Theorem 4.2, we collect some preparatory results: error estimates of the nodal interpolations on the surface, for the Ritz map, and some results which estimate various geometric errors. Most of these results were shown in [12,7,18]. Let us briefly recall our assumptions on the evolving surface and on its discrete counterpart, from Section 2 and 3.1: Γ (t) is a closed smooth (at least C 2 ) surface in R d+1 with d ≤ 3, evolving with the surface velocity v with v, ∂ • v ∈ L ∞ (0, T ; W k+1,∞ (Γ (t))). The discrete surface Γ h (t) is a k-order interpolation of Γ (t) at each time, and therefore its velocity V h is the nodal interpolation of v.

Interpolation error estimates
The following result gives estimates for the error in the interpolation. Our setting follows that of Section 2.5 of [7].
Let us assume that the surface Γ (t) is approximated by the interpolation surface Γ h (t) of order k. Then for any u ∈ H k+1 (Γ (t)), there is a unique k-order surface finite element interpolation I h u ∈ S h (t), furthermore we set ( I h u) = I h u. Lemma 6.1 Let u(·, t) ∈ H k+1 (Γ (t)) for all 0 ≤ t ≤ T . The surface interpolation operator of I h of order k satisfies the following error estimates, for u = u(·, t) and for 0 ≤ t ≤ T , with a constant c > 0 independent of h and t, but depending on G T .

Discrete surface velocities
This section gives a definition of a discrete velocity on the exact surface Γ (t) associated to V h , and explores approximation results for the discrete velocities.
The following result, recalled from [3, Lemma 3.1.6], shows boundedness of the discrete velocity V h , using the fact that it is the interpolation of v. The proof is based on the interpolation error estimate Lemma 6.1 and the interchange formulas (5.3).
Lemma 6.2 Assume that v and ∂ • v are in W k+1,∞ (Γ (t)). Then, for h ≤ h 0 sufficiently small, the following bounds hold: where the constant c > 0 is independent of h and t, but depends on G T .
We give here an explicit formula for v h the discrete surface velocity of Γ (t) associated to V h : for with y(x(t), t) ∈ Γ (t) denoting the unique solution to the lift, i.e. to x(t) = y(x(t), t) + d(x(t), t)ν(y(x(t), t), t) for x(t) ∈ Γ h (t). For more details see [12,Definition 4.3]. Apart from the original material derivative ∂ • on Γ (t), a discrete material derivative associated to the velocity v h is also defined, see [12, equation (4.9)], (element-wise) by That is we have the following three different material derivatives: We note here that it will be always clear from the context whether the discrete From [18,Lemma 5.4] we recall high-order error bounds between the velocity v h of the lifted material points and the surface velocity v (for the case k = 1, l = 0 we refer to [12]). Lemma 6.3 The difference between the continuous velocity v and the discrete velocity v h on Γ (t) is estimated by for l ≥ 0, with a constant c l > 0 independent of h and t, but depending on the surface velocity v.
Since we need to establish a bound for the discrete material derivatives of both defects d u and d w , we recall some transport formulas from [12,Lemma 4.2]: These formulas will help us to derive equations for ∂ • h d u and ∂ • h d w , which can be readily estimated, and will also be often used in the proofs in Section 6.1.4.

Error estimates for the generalised Ritz map
From [18, Theorem 6.3 and 6.4] we recall that the generalised Ritz map (5.4) satisfies the following optimal high-order error estimates. Lemma 6.4 Let u : G T → R such that u(·, t) and (∂ • ) (j) u(·, t) ∈ H k+1 (Γ (t)) for all 0 ≤ t ≤ T and j = 1, . . . , l. Then the error in the generalised Ritz map (5.4) satisfies the bounds, for 0 ≤ t ≤ T and for h ≤ h 0 with sufficiently small h 0 , where the constant c > 0 is independent of h and t, but depends on G T .

Geometric approximation errors
The time dependent bilinear forms m, g and their discrete counterparts m h , g h , from (2.7) and (3.5), respectively, satisfy the following high-order geometric approximation estimates, see [18,Lemma 5.6].
Then, for all h ≤ h 0 with h 0 sufficiently small, the following estimates hold where the constant c > 0 is independent of h and t, but depends on G T .
Similar results hold for the errors in the bilinear form a, cf. [18,Lemma 5.6], but these are not used herein. The above results furthermore imply . (6.4) Below we present and prove a new geometric approximation estimate which relates time derivatives of g and g h .
, with their corresponding lifts in S h (t). Then, for all h ≤ h 0 with h 0 sufficiently small, the following estimate holds for all h ≤ h 0 with h 0 sufficiently small, where the constant c > 0 is independent of h and t, but depends on the surface velocity v.
Proof Although, this lemma was first proved in [3, Lemma 3.1.8], due to its importance we present it here in full detail. Let µ h denote the quotient of the measures on Γ (t) and Γ h (t), we then start by differentiating the integral transformation with respect to time using the transport formulae (6.3), to obtain Using ∂ • using the geometric error estimate (∂ • h ) (2) µ h L ∞ ≤ ch k+1 , cf. [18,Lemma 5.2]. Altogether by triangle inequalities and by combining the above estimates, we obtain where we have used the bounds on the discrete velocity from Lemma 6.2, and the geometric estimate µ h − 1 L ∞ ≤ ch k+1 , from [18, Lemma 5.2].

Defect bounds
In this section we prove bounds for the defects and for their time derivatives, i.e. we prove that condition (5.8) of Proposition 5.1 is indeed satisfied.
where the constant c > 0 is independent of h and t, but depends on the bounds on Sobolev norms of u, w and the surface velocity v.
Proof The Ritz map (5.4) of the exact solutions u and w satisfies the discrete problem only up to some defects, d u (·, t) ∈ S h (t) and d w (·, t) ∈ S h (t), defined in (5.5). Rewriting these equations using the bilinear form notation from (2.7), we thus have, for an arbitrary ϕ h ∈ S h (t), Upon subtracting the corresponding equations for the exact solution (2.6) with ϕ = ϕ h and applying the transport formula (6.3a) (with ∂ • h ϕ h = 0), from the equations in (6.7), and then adding and subtracting some terms in order to apply the definition of the Ritz map R h (5.4), we obtain the following two equations satisfied by the defects d u and d w : We now estimate the defects and their material derivatives in the L 2 (Γ (t)) norm by bounding each pair on the right-hand sides of the above equations separately, using the geometric estimates from the previous subsection, and using similar techniques as in [12,18]. Since throughout the proofs most norms are on Γ (t), we will omit these below, and write L 2 , H k+1 instead of L 2 (Γ (t)), H k+1 (Γ (t)) etc. Bound for d u : For the pair in the first line, we add and subtract terms to obtain where we have used Lemma 6.5 together with the fact that ∂ • h (η ) = (∂ • h η) , the Ritz map error bound Lemma 6.4. The Ritz map error estimate is again used to show the bound ∂ By the same techniques, we prove the following bound for II u : The third term III u is estimated using similar arguments as before, by Lemma 6.5, Lemma 6.4, and Lemma 6.3, as follows The estimates (6.9)-(6.11) together, using the norm equivalence (3.6), and the definition of the L 2 (Γ h (t)) norm yields We start by differentiating the defect equation for d u (6.8a) with respect to time.
The first term is immediately bounded, using Lemma 6.2, the Cauchy-Schwarz inequality and (6.12), by The terms differentiated in time are estimated separately, using analogous techniques as before.
For the first term, by the transport formulas (6.3a) and (6.3b), we obtain 14) where for the inequality we used the arguments used to show (6.9) and (6.11).
By the same arguments, for the second term we obtain the bound By the time differentiation of the third term, using the transport formulas (6.3a) and (6.3b), we obtain The pair in the third line is estimated by previous arguments just as before, byİ The remaining pair in the rectangular brackets is estimated by similar ideas as above, adding and subtracting intermediate terms, using the geometric approximation estimate from Lemma 6.6, Ritz map error estimates Lemma 6.4, and Lemma 6.3 (for l = 0 and 1) and Lemma 6.2: (6.17) The combination of the estimates (6.13)-(6.17), using the norm equivalence (3.6), yields The L 2 norm of the defect d w (6.8b) is estimated by the same techniques by which the bound (6.10) was shown. Due to the (locally Lipschitz continuous) non-linear term W we need an L ∞ bound on the Ritz map, which we obtain by 19) using an inverse estimate [4,Theorem 4.5.11], interpolation error bounds Lemma 6.1, and for the last term the (sub-optimal) interpolation error estimate of [7, Proposition 2.7] (with p = ∞).
By similar techniques as before, and using (6.4) together with (6.19) the pairs for d w are estimated analogously. The bounds for I w and III w are straightforward, while II w is bounded, using the local Lipschitz continuity of W and (6.19), by (6.20) We altogether obtain the estimate Bound for ∂ • h d w : Just as for ∂ • h d u , we differentiate the expression (6.8b) with respect to time. Using again ∂ The first term is estimated using (6.21), while the remaining terms are bounded similarly to (6.15) and (6.20) (using (6.19)).
Altogether, we obtain upon recalling that u * h = R h u and w * h = R h w. The first terms in each error are directly and similarly bounded by error estimates for the Ritz map Lemma 6.4 -uniformly in time -by The second terms are the errors e u h , e w h and ∂ • h e u h , therefore bounded by the combination of the stability estimate (5.9) and the consistency estimates (6.6). In Proposition 5.1 the L ∞ norm assumption on u * h = R h u was proved in (6.19). Altogether, we obtain By combining the above estimates we obtain the stated error estimates.

Extensions to non-linear problems
In this section we will expand on Cahn-Hilliard equations with a more general non-linearity, such as the ones briefly discussed after Theorem 4.2, which includes Cahn-Hilliard equations with a proliferation term [27, equation (3.1)], and which is just slightly less general than the fourth-order non-linear parabolic PDE in [8].
We will now allow the non-linearities to depend on the gradients of u and w, and consider the general non-linear Cahn-Hilliard equation: The corresponding modified evolving surface finite element semi-discretisation of the above problem can be written in the matrix-vector formulation: with the same correction term ϑ as before, see (3.12). The non-linear terms g(u(t)) and f (u(t)) are then defined by The stability analysis of this semi-linear problem only requires the following modifications in the proof of Proposition 5.1: -Since the non-linear terms now include the gradient of u as well, we need to establish a uniform-in-time bound of the W 1,∞ (Γ h ) norm of the errors in u. This is achieved analogously to (5.10), though the definition of t * should now include a W 1,∞ norm bound for the errors e u h , but otherwise the exact same as (5.10). It is key to note that the arguments of (5.36) remain valid, in particular the first inequality holds between · W 1,∞ (Γ h ) and · H 1 (Γ h ) . -This W 1,∞ (Γ h ) norm bound allows to estimate the non-linear terms f and g, similarly as it was done for W in (5.15)-(5.16). For instance, the estimates for the non-linear term f then take the form: and similarly for g, while in Part (ii) the time differentiated terms with f are bounded as For these estimates we refer to the general estimates in Part (v) of the proof of Proposition 7.1 in [19]. In Part (ii), terms tested withė w are first rewritten using the arguments of (5.28), and then estimated via the inequalities above, while the time derivatives of the errors e u need to be absorbed during the combination of (ii.a) and (ii.b). Note that the K norms of the errors already appear on the right-hand side of (5.35), and hence cause no further issues during the Gronwall argument. -Note, however, that it is not possible to include w in the non-linear term f . In Part (ii), such a term would lead to an estimate which cannot be bounded without a new critical term ė w K(t) , even with the argument of (5.28). On the other hand, with a few further tricks it would be possible to consider non-linearities g(u, w) which are linear in w, see, e.g. [20], but we do not elaborate on such situations here.
Altogether, the following general stability estimate holds.
The constant C > 0 is independent of t and h, but depends on 0 < < 1, and exponentially on the final time T .
Therefore, subject to a straightforward revision of the consistency analysis, the above general stability result Proposition 8.1 implies optimal-order error estimates for the general non-linear Cahn-Hilliard equation, analogous to Theorem 4.2. Then there exists an h 0 > 0 such that for all h ≤ h 0 the error between the solutions u and w and the degree k evolving surface finite element solutions u h and w h of the modified system satisfy the optimal-order uniform-in-time error estimates in both variables, for 0 ≤ t ≤ T , w h (·, t) − w(·, t) L 2 (Γ (t)) + h w h (·, t) − w(·, t) H 1 (Γ (t)) ≤ Ch k+1 .

Full discretisation via linearly implicit backward difference formulae
We recall the matrix-vector formulation from (3.14): As a time discretisation, we consider the linearly implicit s-step backward differentiation formulae (BDF). For a step size τ > 0, and with t n = nτ ≤ T , the discretised time derivative is determined bẏ while the non-linear term uses an extrapolated value, and reads as: We determine the approximations to the variables u n to u(t n ) and w n to w(t n ) by the fully discrete system of linear equations, for n ≥ s, which is used for the upcoming numerical experiments. The starting values u i and w i (i = 0, . . . , s − 1) are assumed to be given. They can be precomputed using either a lower order method with smaller step sizes, or an implicit Runge-Kutta method. The method is determined by its coefficients, given by The classical BDF method is known to be zero-stable for s ≤ 6 and to have order s; see [16,Chapter V]. This order is retained by the linearly implicit variant using the above coefficients γ j ; cf. [2,1].
The anti-symmetric structure of the system is preserved, and can be observed in (9.3). Since the idea of energy estimates, using the G-stability theory of Dahlquist [6] and the multiplier technique of Nevanlinna & Odeh [28], can be transferred to linearly implicit BDF full discretisations (up to order 5), we strongly expect that Proposition 5.1 translates to the fully discrete case, and so does the convergence result Theorem 4.2; and also the generalisations in Section 8. This is strengthened by the successful application of these techniques to the analogous linearly implicit backward difference methods applied to evolving surface PDEs: [26,24,23] showing optimal-order error bounds for various problems on evolving surfaces. The method was also analysed for various geometric surface flows, for H 1 -regularised surface flows [22], and for mean curvature flow [19], both proving optimal-order error bounds for full discretisations.

Numerical experiments
In this section we report on numerical experiments, using (9.3), for the nonlinear Cahn-Hilliard equation on an evolving surface: -We perform a convergence test for the non-linear Cahn-Hilliard equation with the linear evolving surface FEM and BDF methods of various order. We would like to note here that [15] only presents errors and EOCs for a linear problem (using the linearly implicit Euler method). -We perform the same experiment as Elliott and Ranner in [15, Section 6.2], i.e. we report on the evolution of the Ginzburg-Landau energy along the surface evolution for the non-linear Cahn-Hilliard equation using the first and second order BDF methods.
In the numerical experiments we use the non-linear Cahn-Hilliard equation on an evolving surface (2.4) with the double-well potential, hence the nonlinear term is W (u) = u 3 − u. Formulated as a system the problem reads: with an extra inhomogeneity f (·, t) : Γ (t) → R, chosen such that the exact solution is known to be u(x, t) = e −6t x 1 x 2 , while w is also explicitly known through the second equation of (10.1). The surface Γ (t) evolves timeperiodically from a sphere into an ellipsoid and back. In particular the surface is given as the zero level set of a distance function: with a(t) = 1 + 0.25 sin(2πt). The initial surface Γ (0) = Γ 0 is the unit sphere. The surface evolution is computed using the ODE for the positions (2.1), with v = −V ν, with V = ∂ t d |∇d| and ν = ∇d |∇d| .
For the numerical experiments the ODE was solved numerically by the classical 4th order Runge-Kutta method with the smallest time step size present in the experiment.
Various numerical experiments have been carried out using the same evolving surface, in particular also for the Cahn-Hilliard equation by Elliott and Ranner [15], and for other problems as well, see, for instance [10,26].
The initial value u 0 h is the interpolation of the exact initial value u 0 . For high-order BDF methods the required additional starting values u i h (for i = 1, . . . , q − 1) are taken as the interpolation of the exact values, if they exist, as well or are otherwise computed using a cascade of steps performed by the preceding lower order method.

Convergence experiments
In these experiments we have used the parameter ε = 0.5. The final time is T = 1, the time discretisations use a sequence of time step sizes τ = 0.2 · 2 −i for i = 1, . . . , 7, and a sequence of initial meshes with (roughly quadrupling) degrees of freedom as reported in the figures.
In Figures 10.1-10.4 we report on the L ∞ (L 2 ) norm errors (left) and L ∞ (H 1 ) norm errors (right) between the numerical and exact solution for both variables u and w, i.e. the plots show the errors where the norms are understood as u − u h L ∞ (L 2 ) = max 0≤nτ ≤T u(·, nτ ) − (u n h ) L 2 (Γ (nτ )) .
For the first order BDF method, Figure 10.1 shows logarithmic plots of the errors against the mesh width h, the lines marked with different symbols correspond to different time step sizes. We also report on temporal convergence in Figure 10.2, where the roles are reversed, the errors are plotted against the time step size τ , and the lines with different markers correspond to different mesh refinements.
In Figure 10.1 we can observe two regions: a region where the spatial discretisation error dominates, matching to the order of convergence of our theoretical results of Theorem 4.2 (note the reference lines), and a region, with small mesh widths, where the temporal discretisation error dominates (the error curves flatten out). For the H 1 norm we observe better convergence rates as predicted. For Figure 10.2, the same description applies, but with reversed roles. Although, we do not study convergence of full discretisations, the classical order of the BDF methods can be observed. We note here, that flat error curves, which were completely dominated by a discretisation error, were not plotted.  The numerical experiments in [15, Section 6.2] reporting on the Ginzburg-Landau energy were repeated here for high-order BDF methods. We again consider the non-linear Cahn-Hilliard equation (10.1), with ε = 0.1 and with f = 0 on the same evolving surface Γ (t) as before, but with a(t) = 1 + 0.25 sin(10πt), and with initial value u 0 (x) = 0.1 cos(2πx 1 ) cos(2πx 2 ) cos(2πx 3 ).
This setting is the same as in [15,Section 6.2].
In Figure 10.5 and 10.6 we report on the time evolution of the Ginzburg-Landau energy (until T = 0.2 and T = 1) of the BDF2 / linear ESFEM discretisation. In both plots we have used the time step size τ = 10 −4 (the same as [15, Section 6.2]), and eight different mesh refinement levels (higher numbering denotes finer meshes). The meshes are not nested refinements of a single coarse grid, the coarsest has 54 while the finest has 10146 nodes.
As it was pointed out by Elliott and Ranner [15] "the energy does not decrease monotonically along solutions", see Figure 10.5, and as they predicted the solutions converge to a time-periodic solution, the periodicity in their energies can be nicely observed in Figure 10