# High-order evolving surface finite element method for parabolic problems on evolving surfaces

High-order evolving surface finite element method for parabolic problems on evolving surfaces Abstract High-order spatial discretizations and full discretizations of parabolic partial differential equations on evolving surfaces are studied. We prove convergence of the high-order evolving surface finite element method by showing high-order versions of geometric approximation errors and perturbation error estimates and by the careful error analysis of a modified Ritz map. Furthermore, convergence of full discretizations using backward difference formulae and implicit Runge–Kutta methods are also shown. 1. Introduction Numerical methods for partial differential equations (PDEs) on stationary and evolving surfaces and for coupled bulk–surface PDEs have been under intensive research in recent years. Surface finite element methods are all based on the fundamental article of Dziuk (1988), further developed for evolving surface parabolic problems by Dziuk & Elliott (2007, 2013b). High-order versions of various finite element methods for problems on a stationary surface have received attention in a number of publications previously. We give a brief overview of this literature here: The surface finite element method of Dziuk (1988) was extended to higher-order finite elements on stationary surfaces by Demlow (2009) Some further important results for higher-order surface finite elements were shown by Elliott & Ranner (2013). Discontinuous Galerkin methods for elliptic surface problems were analysed by Dedner et al. (2013) and then extended to high-order discontinuous Galerkin methods in Antonietti et al. (2015). Recently, unfitted (also called trace or cut) finite element methods have been investigated intensively (see, e.g., Olshanskii et al., 2009; Burman et al., 2015; Reusken, 2015). A higher-order version of the trace finite element method was analysed by Grande & Reusken (2014). However, to our knowledge, there are no articles showing convergence of the high-order evolving surface finite element method for parabolic partial differential equations on evolving surfaces. In this article, we extend the $$H^1$$- and $$L^2$$-norm convergence results of Dziuk & Elliott (2007, 2013b) to high-order evolving surface finite elements applied to parabolic problems on evolving surfaces with prescribed velocity. Furthermore, convergence results for full discretizations using Runge–Kutta methods, based on Dziuk et al. (2012), and using backward difference formulae (BDF), based on Lubich et al. (2013), are also shown. To prove high-order convergence of the spatial discretization, three main groups of errors have to be analysed: Geometric errors, resulting from the appropriate approximation of the smooth surface. Many of these results carry over from Demlow (2009) by careful investigation of time dependencies, while others are extended from Dziuk et al. (2012), Mansour (2013) and Lubich et al. (2013). Perturbation errors of the bilinear forms, whose higher-order version can be shown by carefully using the core ideas of Dziuk & Elliott (2013b). High-order estimates for the errors of a modified Ritz map, which was defined in Lubich & Mansour (2015). These projection error bounds rely on the nontrivial combination of the mentioned geometric error bounds and on the well-known Aubin–Nitsche trick. We further show convergence results for full discretizations using Runge–Kutta and BDF methods. The error estimates, based on energy estimates, for Runge–Kutta methods shown in Dziuk et al. (2012) and for BDF methods in Lubich et al. (2013) are applicable without any modifications, because the semidiscrete problem can be written in a matrix–vector formulation (cf. Dziuk et al., 2012), where the matrices have exactly the same properties as in the linear finite element case. Therefore, the fully discrete convergence results transfer to high-order evolving surface finite elements using the mentioned error estimates of the Ritz map. The implementation of the high-order method is also a nontrivial task. The matrix assembly of the time-dependent mass and stiffness matrices is based on the usual reference element technique. Similar to isoparametric finite elements, the approximating surface is parameterized over the reference element, and therefore all the computations are done there. It was pointed out by Grande & Reusken (2014) that the approach of Demlow (2009) requires explicit knowledge of the exact signed distance function to the surface$${\it{\Gamma}}$$. However, the signed distance function is used only in the analysis, but it is not required for the numerical computations away from the initial time level. This is possible since the high-order element is uniquely determined by its elements, and by using the usual reference element technique. It is used only for generating the initial surface approximation. Here, we consider only linear parabolic PDEs on evolving surfaces; however, we believe our techniques and results carry over to other cases, such as to the Cahn–Hilliard equation (Elliott & Ranner, 2015), to wave equations (Lubich & Mansour, 2015), to ALE methods (Elliott & Venkataraman, 2014; Kovács & Power Guerra, 2014), nonlinear problems (Kovács & Power Guerra, 2016) and to evolving versions of bulk–surface problems (Elliott & Ranner, 2013). For more details, see a remark later on. Furthermore, while, in this article, we consider only evolving surfaces with prescribed velocity, many of the high-order geometric estimates of this article are essential for the numerical analysis of parabolic problems where the surface velocity depends on the solution (cf. Kovács et al., 2016). The article is organized in the following way. In Section 2, we recall the basics of linear parabolic problems on evolving surfaces together with some notation based on Dziuk & Elliott (2007). Section 3 deals with the description of higher-order evolving surface finite elements based on Demlow (2009)Section 4 contains the time discretizations and states the main results of this article: semidiscrete and fully discrete convergence estimates. In Section 5, we turn to the estimates of geometric errors and geometric perturbation estimates. In Section 6, the errors in the generalized Ritz map and its material derivatives are estimated. Section 7 contains the proof of the main results. Section 8 briefly describes the implementation of the high-order evolving surface finite elements. In Section 9, we present some numerical experiments illustrating our theoretical results. 2. The problem Let us consider a smooth-evolving compact surface, given by a smooth signed distance function $$d$$, $${\it{\Gamma}}(t) = \{ x\in \mathbb{R}^{m+1} \mid d(x,t) = 0 \} \subset \mathbb{R}^{m+1}$$ ($$m\leq 3$$), $$0 \leq t \leq T$$, which moves with a given smooth velocity $$v$$. Let $$\partial^{\bullet} u = \partial_{t} u + v \cdot \nabla u$$ denote the material derivative of the function $$u$$, where $$\nabla_{{\it{\Gamma}}}$$ is the tangential gradient given by $$\nabla_{{\it{\Gamma}}} u = \nabla u -\nabla u \cdot \mathrm{n} \mathrm{n}$$, with unit normal $$\mathrm{n}$$. We are sharing the setting of Dziuk & Elliott (2007) and Dziuk & Elliott (2013b). We consider the following linear problem on the above surface, for $$u=u(x,t)$$:   ∂∙u+u∇Γ(t)⋅v−ΔΓ(t)u =f on Γ(t),u(⋅,0) =u0 on Γ(0), (2.1) where for simplicity we set $$f=0$$, but all our results hold with a nonvanishing inhomogeneity as well. Let us briefly recall some important concepts used later on, whereas, in general, for basic formulae we refer to Dziuk & Elliott (2013a). An important tool is the Green’s formula on closed surfaces,   ∫Γ(t)∇Γ(t)z⋅∇Γ(t)ϕ=−∫Γ(t)(ΔΓ(t)z)ϕ. We use Sobolev spaces on surfaces: for a smooth surface $${\it{\Gamma}}(t)$$, for fixed $$t\in[0,T]$$ or for the space–time manifold given by $$\mathscr{G}_T= \cup_{t\in[0,T]} {\it{\Gamma}}(t)\times\{t\}$$, we define   H1(Γ(t))= {η∈L2(Γ(t))∣∇Γ(t)η∈L2(Γ(t))m+1},H1(GT)= {η∈L2(GT)∣∇Γ(t)η∈L2(Γ(t))m+1,∂∙η∈L2(Γ(t))} and analogously for higher-order versions $$H^{k}({\it{\Gamma}}(t))$$ and $$H^{k}(\mathscr{G}_T)$$ for $$k\in \mathbb{N}$$ (cf. Dziuk & Elliott, 2007, Section 2.1). The variational formulation of this problem reads as follows: find $$u\in H^1(\mathscr{G}_T)$$ such that   ddt∫Γ(t)uφ+∫Γ(t)∇Γ(t)u⋅∇Γ(t)φ=∫Γ(t)u∂∙φ (2.2) holds for almost every $$t\in(0,T)$$ for every $$\varphi(\cdot,t) \in H^1({\it{\Gamma}}(t))$$ with $$\partial^{\bullet}\varphi(\cdot,t) \in L^2({\it{\Gamma}}(t))$$ and $$u(\cdot,0)=u_0$$ holds. For suitable $$u_{0}$$, existence and uniqueness results for (2.2) were obtained in Dziuk & Elliott (2007, Theorem 4.4). 3. High-order evolving surface finite elements We define the high-order evolving surface finite element method (ESFEM) applied to our problem following Dziuk (1988), Dziuk & Elliott (2007) and Demlow (2009) We use simplicial elements and continuous piecewise polynomial basis functions of degree $$k$$. 3.1 Basic notions The smooth initial surface $${\it{\Gamma}}(0)$$ is approximated by a $$k$$-order interpolating discrete surface constructed in (a) below. This high-order approximation surface is then evolved in time by the a priori known surface velocity $$v$$, detailed in (b). The construction presented here is from Demlow (2009, Section 2). (a) Let $${\it{\Gamma}}_h^1(0)$$ be a polyhedron having triangular elements (denoted by $$E$$) with vertices lying on the initial surface $${\it{\Gamma}}(0)$$, forming a quasiuniform triangulation $$\mathscr{T}_h^1(0)$$, with mesh parameter $$h$$, and unit outward normal $$\mathrm{n}_h^1$$. Next, for $$k \geq 2$$, we define $${\it{\Gamma}}_h^k(0)$$, the $$k$$-order polynomial approximations to $${\it{\Gamma}}(0)$$. For a given fixed element $$E \in \mathscr{T}_h^1(0)$$, with nodes (numbered locally) $$x_1, x_2, \dotsc, x_{n_k}$$, the corresponding Lagrange basis functions of degree $$k$$ on $$E$$ are denoted by $$\tilde\chi^k_1, \tilde\chi^k_2, \dotsc, \tilde\chi^k_{n_k}$$. For arbitrary $$x \in E$$, we define the discrete projection   pk(x,0)=∑j=1nkp(xj,0)χ~jk(x,0), with Γ(0)∋p(xj,0)=xj−d(xj,0)n(p(xj,0)). This definition yields a continuous piecewise polynomial map on $${\it{\Gamma}}_h^1(0)$$. Then the $$k$$-order approximation surface is defined by   Γhk(0)={pk(x,0)∣x∈Γh1(0)}. The vertices of $${\it{\Gamma}}_h^k(0)$$ are denoted by $$a_i(0)$$, $$i=1,2,\dotsc,N$$ (here numbered globally). Note that these vertices are sitting on the exact surface. The high-order triangulation is denoted by $$\mathscr{T}_h^k(0)$$. Further small details can be found in Demlow (2009). (b) The surface approximation $${\it{\Gamma}}_h^k(t)$$ at time $$t\in[0,T]$$ is given by evolving the nodes of the initial triangulation by the velocity $$v$$ along the space-time manifold. For instance, the nodes $$a_i(t)$$ are determined by the ordinary differential equation (ODE) system   ddtai(t)=v(ai(t),t). Then the base triangulation and the basis functions move along as well. Hence, for $$0 \leq t \leq T$$, the nodes $$(a_i(t))_{i=1}^N$$ define an approximation of $${\it{\Gamma}}(t)$$ of degree $$k$$,   Γhk(t)={pk(x,t)∣x∈Γh1(t)}, defined analogously as for $$t=0$$. Therefore, the discrete surface $${\it{\Gamma}}_h^k(t)$$ remains an interpolation of $${\it{\Gamma}}(t)$$ for all times. The high-order triangulation at $$t$$ is denoted by $$\mathscr{T}_h^k(t)$$. Remark 3.1 It is highly important to note here that computing $$p(\cdot,t)$$ or $$p^k(\cdot,t)$$, for $$t>0$$, is never used during the numerical computations, hence explicit knowledge of the distance function is avoided (except for $$t=0$$). Instead, the positions of the nodes appearing in an element of $${\it{\Gamma}}_h^k(t)$$ are used to construct a reference mapping over the reference triangle $$E_0$$. In fact, the nodes (locally numbered here) $$a_1(t), a_2(t), \dotsc, a_{n_k}(t)$$ of the high-order element uniquely determine an approximating surface of degree $$k$$ (being an image of a polynomial of $$m$$ variables over $$E_0$$). The ESFEM matrices are then assembled using the reference element and the reference mapping. The discrete tangential gradient on the discrete surface $${\it{\Gamma}}_h^k(t)$$ is given by   ∇Γhk(t)ϕ:=∇ϕ−∇ϕ⋅nhknhk, understood in an elementwise sense, with $$\mathrm{n}_h^k$$ denoting the normal to $${\it{\Gamma}}_h^k(t)$$. For every $$t\in[0,T]$$ and for the $$k$$-order approximation surface $${\it{\Gamma}}_h^k(t)$$, we define the finite element subspace of order $$k$$, denoted by $$S_h^k(t)$$, spanned by the continuous evolving basis functions $$\chi_j$$ of piecewise degree $$k$$, satisfying $$\chi_j(a_i(t),t) = \delta_{ij}$$ for all $$i,j = 1, 2, \dotsc, N$$:   Shk(t)= {ϕ∈C(Γhk(t))∣ϕ=ϕ~∘pk(⋅,t)−1 where ϕ~∈S~hk(t)}= span{χ1(⋅,t),χ2(⋅,t),…,χN(⋅,t)}, where   χj(⋅,t)=χ~j(pk(⋅,t)−1,t)(j=1,2,…,N), with $$\tilde\chi_j( \cdot ,t)$$ being the degree $$k$$ polynomial basis function over the base triangulation $${\it{\Gamma}}_h^1$$, spanning the space $$\tilde S_h^k(t) = \{ \tilde\phi \in C({\it{\Gamma}}_h^1(t)) \mid \tilde\phi \text{ is a piecewise polynomial of degree } k \}$$. By construction $$S_h^k(t)$$ is an isoparametric finite element space. We interpolate the surface velocity on the discrete surface using the basis functions and denote it by $$V_h$$. Then the discrete material derivative is given by   ∂h∙ϕh=∂tϕh+Vh⋅∇ϕh(ϕh∈Shk(t)). The key transport property derived by Dziuk & Elliott (2007, Proposition 5.4) is   ∂h∙χj=0forj=1,2,…,N. (3.1) This result translates from the linear ESFEM case by the original proof of Dziuk & Elliott (2007, Section 5.2) using barycentric coordinates and the chain rule. The spatially discrete problem (for a fixed degree $$k$$) then reads as follows: find $$U_h\in S_h^k(t)$$, with $$\partial^{\bullet}_h U_h\in S_h^k(t)$$ and temporally continuous, such that   ddt∫Γhk(t)Uhϕh+∫Γhk(t)∇ΓhUh⋅∇Γhϕh=∫Γhk(t)Uh∂h∙ϕh(∀ϕh∈Shk(t) with ∂h∙ϕh∈Shk(t)), (3.2) with the initial condition $$U_h^0\in S_h^k(0)$$ being a suitable approximation to $$u_0$$. Later on we will always work with a high-order approximation surface $${\it{\Gamma}}_h^k(t)$$ and with the corresponding high-order evolving surface finite element space $$S_h^k(t)$$ (with $$2\leq k\in\mathbb{N}$$); therefore from now on, we drop the upper index $$k$$, unless we would like to emphasize it or it is not clear from the context. 3.2 The ODE system Similarly to Dziuk et al. (2012), the ODE form of the above problem (3.2) can be derived by setting   Uh(⋅,t)=∑j=1Nαj(t)χj(⋅,t) in the semidiscrete problem, and by testing with $$\phi_h=\chi_j$$ ($$j=1,2,\dotsc,N$$), and using the transport property (3.1). The spatially semidiscrete problem (3.2) is equivalent to the following ODE system for the vector $$\alpha(t)=(\alpha_j(t))_{j=1}^N \in\mathbb{R}^N$$, collecting the nodal values of $$U_h(\cdot,t)$$:   ddt(M(t)α(t))+A(t)α(t) =0,α(0) =α0, (3.3) where the evolving mass matrix $$M(t)$$ and the stiffness matrix $$A(t)$$ are defined as   M(t)|kj=∫Γh(t)χjχkandA(t)|kj=∫Γh(t)∇Γhχj⋅∇Γhχk, for $$j,k = 1,2, \dotsc, N$$. 3.3 Lifts For the error analysis we need to compare functions on different surfaces. This is conveniently done by the lift operator, which was introduced in Dziuk (1988) and further investigated in Dziuk & Elliott (2007, 2013b). The lift operator maps a function on the discrete surface onto a function on the exact surface. Let $${\it{\Gamma}}_h(t)$$ be a $$k$$-order approximation to the exact surface $${\it{\Gamma}}(t)$$. Using the oriented distance function $$d$$ (cf. Dziuk & Elliott, 2007, Section 2.1), the lift of a continuous function $$\eta_h \colon {\it{\Gamma}}_h(t) \rightarrow \mathbb{R}$$ is defined as   ηhl(p,t):=ηh(x,t),x∈Γ(t), where for every $$x\in {\it{\Gamma}}_h(t)$$ the point $$p=p(x,t)\in{\it{\Gamma}}(t)$$ is uniquely defined via   p=x−n(p,t)d(x,t). (3.4) By $$\eta^{-l}$$ we denote the function on $${\it{\Gamma}}_h(t)$$ whose lift is $$\eta$$. In particular, we will often use the space of lifted basis functions   (Sh(t))l=(Shk(t))l={ϕhl | ϕh∈Shk(t)}. 4. Convergence estimates 4.1 Convergence of the semidiscretization We now formulate the convergence theorem for the semidiscretization using high-order evolving surface finite elements. This result is the higher-order extension of Dziuk & Elliott (2013b, Theorem 4.4). Theorem 4.1 Consider the ESFEM of order $$k$$ as a space discretization of the parabolic problem (2.1). Let $$u$$ be a sufficiently smooth solution of the problem, and assume that the initial value satisfies   ‖uh0−u(⋅,0)‖L2(Γ(0))≤C0hk+1. Then there exists $$h_0>0$$, such that for mesh size $$h \leq h_0$$, the following error estimate holds, for $$t \leq T$$:   ‖uh(⋅,t)−u(⋅,t)‖L2(Γ(t))+h(∫0t‖∇Γ(s)(uh(⋅,s)−u(⋅,s))‖L2(Γ(s))2ds)1/2≤Chk+1. The constant $$C$$ is independent of $$h$$ and $$t$$ but depends on $$T$$. The proof of this result is postponed to a later section, after we have shown some preparatory results. 4.2 Time discretization: BDF We apply a $$p$$-step BDF for $$p\leq5$$ as a discretization of the ODE system (3.3), coming from the ESFEM space discretization of the parabolic evolving surface PDE. We briefly recall the $$p$$-step BDF method applied to system (3.3) with step size $$\tau>0$$:   1τ∑j=0pδjM(tn−j)αn−j+A(tn)αn=0(n≥p), (4.1) where the coefficients of the method are given by $$\delta(\zeta)=\sum_{j=0}^p \delta_j \zeta^j=\sum_{\ell=1}^p \frac{1}{\ell}(1-\zeta)^\ell$$, while the starting values are $$\alpha_0, \alpha_1, \dotsc, \alpha_{p-1}$$. The method is known to be $$0$$-stable for $$p\leq6$$ and have order $$p$$ (for more details, see Hairer & Wanner, 1996, Chapter V). In the following result, we compare the fully discrete solution   Uhn=∑j=1Nαjnχj(⋅,tn), obtained by solving (4.1) and the Ritz map $$\widetilde{\mathscr{P}}_h : H^{1}({\it{\Gamma}}(t)) \rightarrow S_h^k(t)$$ ($$t\in[0,T]$$) of the sufficiently smooth solution $$u$$. The precise definition of the Ritz map is given later. The $$H_h^{-1}({\it{\Gamma}}_h(t))$$-norm of the ESFEM function $$R_h$$ is defined as   ‖Rh(⋅,t)‖Hh−1(Γh(t))=sup0≠ϕh∈Shk(t)mh(Rh(⋅,t),ϕh)‖ϕh‖H1(Γh(t)) . The following error bound was shown in Lubich et al. (2013) for BDF methods up to order 5 (see also Mansour, 2013). Theorem 4.2 (Lubich et al., 2013, Theorem 5.1) Consider the parabolic problem (2.1), having a sufficiently smooth solution for $$0\leq t \leq T$$. Couple the $$k$$-order ESFEM as space discretization with time discretization by a $$p$$-step backward difference formula with $$p\leq5$$. Assume that the Ritz map of the solution has continuous discrete material derivatives up to order $$p+1$$. Then there exists $$\tau_0>0$$, independent of $$h$$, such that for $$\tau \leq \tau_0$$, for the error $$E_h^n=U_h^n-\widetilde{\mathscr{P}}_h u(\cdot,t_n)$$ the following estimate holds for $$t_n=n\tau \leq T$$:   ‖Ehn‖L2(Γh(tn))+ (τ∑j=1n‖∇Γh(tj)Ehj‖L2(Γh(tj))2)1/2≤Cβ~h,pτp+ (τ∑j=1n‖Rh(⋅,tj)‖Hh−1(Γh(tj))2)1/2+Cmax0≤i≤p−1‖Ehi‖L2(Γh(ti)), where $$R_h$$ is the high-order ESFEM residual. The constant $$C$$ is independent of $$h, n$$ and $$\tau$$ but depends on $$T$$. Furthermore,   β~h,p2=∫0T∑ℓ=1p+1‖(∂h∙)(ℓ)(P~hu)(⋅,t)‖L2(Γh(t))dt. Remark 4.3 The most important technical tools in the proof of Theorem 4.2, and also in the proof of the BDF stability result in Lubich et al. (2013, Lemma 4.1), are the ODE formulation (3.3) and the key estimates first shown in Dziuk et al. (2012, Lemma 4.1), where the following estimates are shown: there exist $$\mu,\kappa>0$$ such that, for $$w,z \in \mathbb{R}^N$$,   wT(M(s)−M(t))z≤(eμ(s−t)−1)‖w‖M(t)‖z‖M(t),wT(A(s)−A(t))z≤(eκ(s−t)−1)‖w‖A(t)‖z‖A(t). The proof of these bounds involves only basic properties of the mass and stiffness matrices $$M(t)$$ and $$A(t)$$; hence it is independent of the order of the basis functions. Hence, the high-order ESFEM versions of these two inequalities also hold. Therefore, the original results of Lubich et al. (2013) hold here as well. Later on, when suitable results are at hand, we give some remarks on the smoothness assumption of the Ritz map. 4.3 Convergence of the full discretization We are now in a position to formulate one of the main results of this article, which yields optimal-order error bounds for high-order finite element semidiscretization coupled to BDF methods up to order 5 applied to an evolving surface PDE. Theorem 4.4 ($$k$$-order ESFEM and BDF-$$p$$) Consider the ESFEM of order $$k$$ as space discretization of the parabolic problem (2.1), coupled to the time discretization by a $$p$$-step backward difference formula with $$p\leq5$$. Let $$u$$ be a sufficiently smooth solution of the problem and assume that the starting values satisfy (with $$\mathscr{P}_h u=(\widetilde{\mathscr{P}}_h u)^l)$$)   max0≤i≤p−1‖uhi−(Phu)(⋅,ti)‖L2(Γ(ti))≤C0hk+1. Then there exist $$h_0>0$$ and $$\tau_0>0$$, such that for $$h\leq h_0$$ and $$\tau \leq \tau_0$$, the following error estimate holds for $$t_n=n\tau \leq T$$:   ‖uhn−u(⋅,tn)‖L2(Γ(tn))+h(τ∑j=pn‖∇Γ(tj)(uhj−u(⋅,tj))‖L2(Γ(tj))2)1/2≤C(τp+hk+1). The constant $$C$$ is independent of $$h, \ \tau$$ and $$n$$ but depends on $$T$$. The proof of this result is also postponed to a later section, after we have shown some preparatory results. Remark 4.5 We remark here that an analogous fully discrete convergence result is readily available for algebraically stable implicit Runge–Kutta methods (such as the Radau IIA methods), since the Runge–Kutta analogue of Theorem 4.2 has been proved in Dziuk et al. (2012). The combination of this result with our high-order semidiscrete error bounds (Theorem 7.1) proves the Runge–Kutta analogue of Theorem 4.4. 5. Geometric estimates In this section, we present further notation and some technical lemmas that will be used later on in the proofs leading to the convergence result. These estimates are high-order analogues of some previous results proved in Dziuk (1988), Dziuk & Elliott (2007), Demlow (2009), Mansour (2013) and Dziuk & Elliott (2013b). 5.1 Geometric approximation results In the following, we state and prove estimates for the errors resulting from the geometric surface approximation. Most of these estimates hold for a sufficiently small $$h$$, which here means for $$h\leq h_0$$ with a sufficiently small $$h_0>0$$. Lemma 5.1 (Equivalence of norms, Dziuk, 1988; Demlow, 2009) Let $$\eta_h : {\it{\Gamma}}_h(t) \rightarrow \mathbb{R}$$ with lift $$\eta_h^l : {\it{\Gamma}}(t) \rightarrow \mathbb{R}$$. Then, for a sufficiently small $$h$$, the discrete and continuous $$L^p$$ and Sobolev norms are equivalent, independently of the mesh size $$h$$. For instance, there is a constant $$c>0$$ such that for all $$h$$ sufficiently small,   c−1‖ηh‖L2(Γh(t))≤ ‖ηhl‖L2(Γ(t))≤c‖ηh‖L2(Γh(t)),c−1‖ηh‖H1(Γh(t))≤ ‖ηhl‖H1(Γ(t))≤c‖ηh‖H1(Γh(t)). We now turn to the study of some geometric concepts and their errors. By $$\delta_h$$ we denote the quotient between the continuous and discrete surface measures, $$\textrm{d} A$$ and $$\textrm{d} A_h$$, defined as $$\delta_h \textrm{d} A_h = \textrm{d} A$$. Further, we recall that $${\text{Pr}}$$ and $${\text{Pr}}_h$$ are the projections onto the tangent spaces of $${\it{\Gamma}}(t)$$ and $${\it{\Gamma}}_h(t)$$, respectively. We further set, from Dziuk & Elliott (2013b),   Qh=1δh(Id−dH)PrPrhPr(Id−dH), (5.1) where $$\mathscr{H}$$ ($$\mathscr{H}_{ij} = \partial_{x_j}\mathrm{n}_i$$) is the (extended) Weingarten map. Using this notation and (3.4), in the proof of Dziuk & Elliott (2013b, Lemma 5.5), it is shown that   ∇Γhϕh(x)⋅∇Γhϕh(x)=δhQh∇Γϕhl(p)⋅∇Γϕhl(p). (5.2) For these quantities we show some results analogous to their linear ESFEM version showed in Dziuk & Elliott (2013b, Lemma 5.4), Demlow (2009, Proposition 2.3) or Mansour (2013, Lemma 6.1). Later on the following estimates will play a key technical role. Lemma 5.2 For $${\it{\Gamma}}_h^k(t)$$ and $${\it{\Gamma}}(t)$$ as above, for a sufficiently small $$h$$, we have the following geometric approximation estimates:    ‖d‖L∞(Γh(t))≤chk+1, ‖1−δh‖L∞(Γh(t))≤chk+1, ‖n−nh‖L∞(Γh(t))≤chk, and ‖Id−δhQh‖L∞(Γh(t))≤chk+1, ‖(∂h∙)(ℓ)d‖L∞(Γh(t))≤chk+1, ‖(∂h∙)(ℓ)δh‖L∞(Γh(t))≤chk+1, ‖Pr((∂h∙)(ℓ)Qh)Pr‖L∞(Γh)≤chk+1, with constants depending only on $$\mathscr{G}_T$$ but not on $$h$$ or $$t$$. The first three bounds were shown in Demlow (2009, Proposition 2.3) for the stationary case. Noting that the constants depend only on $$\mathscr{G}_T$$ these inequalities are shown. The last four bounds are simply the higher-order extensions of the corresponding estimates of Mansour (2013, Lemma 6.1). They can be proved in the exact same way using the bounds of the first three estimates. These proofs are included in the Appendix. 5.2 Interpolation estimates for evolving surface finite elements The following result gives estimates for the error in the interpolation. Our setting follows that of Demlow (2009, Section 2.5). Let us assume that the surface $${\it{\Gamma}}(t)$$ is approximated by the interpolation surface $${\it{\Gamma}}_h^k(t)$$. Then for any $$w\in H^{k+1}({\it{\Gamma}}(t))$$, there is a unique $$k$$-order surface finite element interpolation $$\widetilde I_h^k w \in S_h^k(t)$$; furthermore we set $$(\widetilde I_h^k w)^l = I_h^k w$$. Lemma 5.3 (Demlow, 2009, Proposition 2.7) Let $$w:\mathscr{G}_T \rightarrow \mathbb{R}$$ such that $$w\in H^{k+1}({\it{\Gamma}}(t))$$ for $$0 \leq t \leq T$$. There exists a constant $$c>0$$ depending on $$\mathscr{G}_T$$, but independent of $$h$$ and $$t$$, such that for $$0 \leq t \leq T$$ and for a sufficiently small $$h$$,   ‖w−Ihkw‖L2(Γ(t))+h‖∇Γ(w−Ihkw)‖L2(Γ(t))≤ chk+1‖w‖Hk+1(Γ(t)). We distinguish the special case of a linear surface finite element interpolation on $${\it{\Gamma}}_h^k(t)$$. For $$w \in H^2({\it{\Gamma}}(t))$$, the linear surface finite element interpolant is denoted by $$I_h^{(1)} w$$, and it satisfies, with $$c>0$$,   ‖w−Ih(1)w‖L2(Γ(t))+h‖∇Γ(w−Ih(1)w)‖L2(Γ(t))≤ ch2‖w‖H2(Γ(t)). Note that for $$I_h^{(1)}$$ the underlying approximating surface $${\it{\Gamma}}_h^k(t)$$ is still of high order. For $$k=1$$, $$I_h^1$$ and $$I_h^{(1)}$$ simply coincide. The upper $$k$$ indices are again dropped later on. 5.3 Velocity of lifted material points and material derivatives Following Dziuk & Elliott (2013b) and Lubich & Mansour (2015) we define the velocity of the lifted material points, denoted by $$v_h$$, and the corresponding discrete material derivative $$\partial^{\bullet}_h$$. For arbitrary $$y(t) = p(x(t),t) \in {\it{\Gamma}}(t)$$, with $$x(t) \in {\it{\Gamma}}_h(t)$$, cf. (3.4), we have   ddty(t)=vh(y(t),t)=∂tp(x(t),t)+Vh(x(t),t)⋅∇p(x(t),t); (5.3) hence for $$y=p(x,t)$$ (see Dziuk & Elliott, 2013b),   vh(y,t)=(Pr−dH)(x,t)Vh(x,t)−∂td(x,t)n(x,t)−d(x,t)∂tn(x,t). Following Lubich & Mansour (2015, Section 7.3), we note that $$- \partial_t d(x,t)\mathrm{n}(x,t)$$ is the normal component of $$v(p,t)$$ and that the other terms are tangent to $${\it{\Gamma}}(t)$$; hence   v−vhis a tangent vector. (5.4) It is also important to note that $$v_h \neq V_h^l$$ (cf. Dziuk & Elliott, 2013b). The discrete material derivative of the lifted points on $${\it{\Gamma}}(t)$$ reads   ∂h∙φh=∂tφh+vh⋅φh(φh∈(Sh(t))l). The lifted basis functions also satisfy the transport property   ∂h∙χjl=0(j=1,2,…,N). To prove error estimates for higher-order material derivatives of the Ritz map, we need high-order bounds for the error between the continuous velocity $$v$$ and the discrete velocity $$v_h$$. We generalize here Dziuk & Elliott (2013b, Lemma 5.6) and Lubich & Mansour (2015, Lemma 7.3) to the high-order case. Lemma 5.4 For $$\ell \geq 0$$, there exists a constant $$c_\ell>0$$ depending on $$\mathscr{G}_T$$, but independent of $$t$$ and $$h$$, such that, for a sufficiently small $$h$$,   ‖(∂h∙)(ℓ)(v−vh)‖L∞(Γ(t))+h‖∇Γ(∂h∙)(ℓ)(v−vh)‖L∞(Γ(t))≤cℓhk+1. Proof. We follow the steps of the original proofs from Dziuk & Elliott (2013b) and Lubich & Mansour (2015). (a) For $$\ell=0$$. Using the definition (5.3) and the fact $$V_h = \widetilde I_h v$$, we have   |v(p,t)−vh(p,t)|=|Pr(v−Ihv)(p,t)+d(HIhv(p,t)+∂tn)|≤chk+1, where we have used the interpolation estimates, Lemma 5.2 and the boundedness of the other terms. For the gradient estimate we use the fact that $$\nabla_{{\it{\Gamma}}} d = \nabla_{{\it{\Gamma}}} \partial^{\bullet}_h d = 0$$ and the geometric bounds of Lemma 5.2:   |∇Γ(v−vh)| ≤c|v−Ihv|+c|∇Γ(v−Ihv)| +|(∇Γd)(HIhv+∂tn)|+|d(∇Γ(HIhv+∂tn))| ≤chk. (b) For $$\ell=1$$, we use the transport property and again Lemma 5.2:   |∂h∙(v−vh)| ≤|(∂h∙Pr)(v−Ihv)|+|Pr(∂h∙v−Ih∂h∙v)| +|(∂h∙d)(HIhv+∂tn)|+|d(∂h∙(HIhv+∂tn))| ≤chk+1. Again using $$\nabla_{{\it{\Gamma}}} d = \nabla_{{\it{\Gamma}}} \partial^{\bullet}_h d = 0$$ and the geometric bounds of Lemma 5.2,   |∇Γ∂h∙(v−vh)| ≤c|v−Ihv|+c|∇Γ(v−Ihv)| +c|∂h∙(v−Ihv)|+c|∇Γ(∂h∙v−Ih∂h∙v)| +|(∇Γ∂h∙d)(HIhv+∂tn)|+|d(∇Γ∂h∙(HIhv+∂tn))| ≤chk. (c) For $$\ell >1$$, the proof uses similar arguments. □ 5.4 Bilinear forms and their estimates We use the time-dependent bilinear forms defined in Dziuk & Elliott (2013b): for arbitrary $$z,\varphi \in H^1({\it{\Gamma}}(t))$$ and for their discrete analogues for $$Z_h, \phi_h \in S_h(t)$$,   m(t;z,φ) =∫Γ(t)zφ,a(t;z,φ) =∫Γ(t)∇Γz⋅∇Γφ,g(t;v;z,φ) =∫Γ(t)(∇Γ⋅v)zφ,b(t;v;z,φ) =∫Γ(t)B(v)∇Γz⋅∇Γφ,mh(t;Zh,ϕh) =∫Γh(t)Zhϕh,ah(t;Zh,ϕh) =∫Γh(t)∇ΓhZh⋅∇Γhϕh,gh(t;Vh;Zh,ϕh) =∫Γh(t)(∇Γh⋅Vh)Zhϕh,bh(t;Vh;Zh,ϕh) =∫Γh(t)Bh(Vh)∇ΓhZh⋅∇Γhϕh, where the discrete tangential gradients are understood in a piecewise sense, and with the tensors given by   B(v)|ij =δij(∇Γ⋅v)−((∇Γ)ivj+(∇Γ)jvi),Bh(Vh)|ij =δij(∇Γh⋅Vh)−((∇Γh)i(Vh)j+(∇Γh)j(Vh)i), for $$i,j=1,2,\dotsc,m+1$$. For more details, see Dziuk & Elliott (2013b, Lemma 2.1) (and the references in the proof) or Dziuk & Elliott (2013a, Lemma 5.2). We will omit the time dependency of the bilinear forms if it is clear from the context. 5.4.1 Discrete material derivative and transport properties The following transport equations, especially the one with discrete material derivatives on the continuous surface, are of great importance during the defect estimates later on. Lemma 5.5 (Dziuk & Elliott, 2013b, Lemma 4.2) Consider $${\it{\Gamma}}(t)$$ as the lift of the discrete surface $${\it{\Gamma}}_h^k(t)$$ (i.e., $${\it{\Gamma}}(t)$$ can be decomposed into curved elements that are lifts of the elements of $${\it{\Gamma}}_h^k(t)$$), moving with the velocity $$v_h$$ from (5.3). Then for any $$z,\varphi, \partial^{\bullet}_h z, \partial^{\bullet}_h \varphi \in H^1({\it{\Gamma}}(t))$$ we have   ddtm(z,φ) =m(∂h∙z,φ)+m(z,∂h∙φ)+g(vh;z,φ),ddta(z,φ) =a(∂h∙z,φ)+a(z,∂h∙φ)+b(vh;z,φ). The same formulae hold for $${\it{\Gamma}}(t)$$ when $$\partial^{\bullet}_h$$ and $$v_h$$ are replaced with $$\partial^{\bullet}$$ and $$v$$, respectively. Similarly, in the discrete case, for arbitrary $$z_h, \phi_h, \partial^{\bullet}_hz_h, \partial^{\bullet}_h\phi_h \in S_h(t)$$, we have   ddtmh(zh,ϕh) =mh(∂h∙zh,ϕh)+mh(zh,∂h∙ϕh)+gh(Vh;zh,ϕh),ddtah(zh,ϕh) =ah(∂h∙zh,ϕh)+ah(zh,∂h∙ϕh)+bh(Vh;zh,ϕh). 5.4.2 Geometric perturbation errors The following estimates are the most important technical results of this paper. Later on, they will play a crucial role in the defect estimates. We note here that these results extend the first-order ESFEM theory of Dziuk & Elliott (2013b) (for the first three inequalities) and Lubich & Mansour (2015) (for the last inequality) to the higher-order ESFEM case. The high-order version of the inequalities for time-independent bilinear forms $$a(\cdot,\cdot)$$ and $$m(\cdot,\cdot)$$ were shown in Elliott & Ranner (2013, Lemma 6.2). The following results generalizes these for the time-dependent case. Lemma 5.6 For any $$Z_h,\phi_h \in S_h^k(t)$$, and for their lifts $$Z_h^l,\phi_h^l \in H^1({\it{\Gamma}}(t))$$, we have the following bounds, with a sufficiently small $$h$$:   |m(Zhl,ϕhl)−mh(Zh,ϕh)| ≤chk+1‖Zhl‖L2(Γ(t))‖ϕhl‖L2(Γ(t)),|a(Zhl,ϕhl)−ah(Zh,ϕh)| ≤chk+1‖∇ΓZhl‖L2(Γ(t))‖∇Γϕhl‖L2(Γ(t)),|g(vh;Zhl,ϕhl)−gh(Vh;Zh,ϕh)| ≤chk+1‖Zhl‖L2(Γ(t))‖ϕhl‖L2(Γ(t)),|b(vh;Zhl,ϕhl)−bh(Vh;Zh,ϕh)| ≤chk+1‖∇ΓZhl‖L2(Γ(t))‖∇Γϕhl‖L2(Γ(t)), where the constants $$c>0$$ are independent of $$h$$ and $$t$$, but depend on $$\mathscr{G}_T$$. Proof. The proof of the first two estimates is a high-order generalization of Dziuk & Elliott (2013b, Lemma 5.5), while the proof of the last two estimates is a high-order extension of Lubich & Mansour (2015, Lemma 7.5). Their proofs follow these references. Both parts use the geometric estimates shown in Lemma 5.2. To show the first inequality, we estimate   |m(Zhl,ϕhl)−mh(Zh,ϕh)| =|∫Γ(t)ZhlϕhldA−∫Γh(t)ZhϕhdAh| =|∫Γ(t)(1−δh−1)ZhlϕhldA| ≤chk+1‖Zhl‖L2(Γ(t))‖ϕhl‖L2(Γ(t)). For the second inequality, similarly we have   |a(Zhl,ϕhl)−ah(Zh,ϕh)| =|∫Γ(t)∇ΓZhl⋅∇ΓϕhldA−∫Γh(t)∇ΓhZh⋅∇ΓhϕhdAh| =|∫Γ(t)(Id−δhQh)∇ΓZhl⋅∇ΓϕhldA| ≤chk+1‖∇ΓZhl‖L2(Γ(t))‖∇Γϕhl‖L2(Γ(t)). For the third estimate we start by taking the time derivative of the equality $$m(Z_h^l,\phi_h^l) = m_h(Z_h,\phi_h \delta_h)$$, using the first transport property from Lemma 5.5 to obtain   ddtm(Zhl,ϕhl) =m(∂h∙Zhl,ϕhl)+m(Zhl,∂h∙ϕhl)+g(vh;Zhl,ϕhl) =ddtmh(Zh,ϕhδh) =mh(∂h∙Zh,ϕhδh)+mh(Zh,(∂h∙ϕh)δh) +gh(Vh;Zh,ϕhδh)+mh(Zh,(∂h∙δh)ϕh). Hence, using $$\partial^{\bullet}_h w_h^l = (\partial^{\bullet}_h w_h)^l$$,   g(vh;Zhl,ϕhl)−gh(Vh;Zh,ϕhδh) = m((∂h∙Zh)l,ϕhl)−mh(∂h∙Zh,ϕhδh) +m(Zhl,(∂h∙ϕh)l)−mh(Zh,(∂h∙ϕh)δh)+mh(Zh,(∂h∙δh)ϕh) =mh(Zh,(∂h∙δh)ϕh). Hence, together with Lemma 5.2, the bound   |g(vh;Zhl,ϕhl)−gh(Vh;Zh,ϕh)| ≤|gh(Vh;Zh,ϕh(δh−1))|+|mh(Zh,(∂h∙δh)ϕh)| ≤c(‖∂h∙δh‖L∞(Γh(t))+‖1−δh‖L∞(Γh(t)))‖Zhl‖L2(Γ(t))‖ϕhl‖L2(Γ(t)) finishes the proof of the third estimate. For the last inequality we take a similar approach, using the second transport property from Lemma 5.5 and the relation (5.2) to obtain   ddtah(Zh,ϕh) =ah(∂h∙Zh,ϕh)+ah(Zh,∂h∙ϕh)+bh(Vh;Zh,ϕh) =ddt∫Γ(t)Qhl∇ΓZhl⋅∇Γϕhl =∫Γ(t)Qhl∇Γ∂h∙Zhl⋅∇Γϕhl+∫Γ(t)Qhl∇ΓZhl⋅∇Γ∂h∙ϕhl +∫Γ(t)(∂h∙Qhl)∇ΓZhl⋅∇Γϕhl+∫Γ(t)B(vh)Qhl∇ΓZhl⋅∇Γϕhl. Again, using $$\partial^{\bullet}_h w_h^l = (\partial^{\bullet}_h w_h)^l$$, together with (5.2) and the geometric estimates of Lemma 5.2 we obtain   |bh(Vh;Zh,ϕh)−b(vh;Zhl,ϕhl)| =|∫Γ(t)(∂h∙Qhl)∇ΓZhl⋅∇Γϕhl+∫Γ(t)B(vh)(Qhl−Id)∇ΓZhl⋅∇Γϕhl| ≤chk+1‖∇ΓZhl‖L2(Γ(t))‖∇Γϕhl‖L2(Γ(t)), completing the proof. □ 6. Generalized Ritz map and higher-order error bounds We recall the generalized Ritz map for evolving surface PDEs from Lubich & Mansour (2015). Definition 6.1 (Ritz map) For any given $$z\in H^{1}({\it{\Gamma}}(t))$$, there is a unique $$\widetilde{\mathscr{P}}_h z\in S_h^k(t)$$ such that for all $$\phi_h\in S_h^k(t)$$, with the corresponding lift $$\varphi_h=\phi_h^l$$, we have   ah∗(P~hz,ϕh)=a∗(z,φh), (6.1) where we let $$a^{\ast}=a+m$$ and $$a_h^{\ast}=a_h+m_h$$, so that the forms $$a$$ and $$a_h$$ are positive definite. Then $$\mathscr{P}_h z \in (S_h^k(t))^l$$ is defined as the lift of $$\widetilde{\mathscr{P}}_h z$$, i.e., $$\mathscr{P}_h z = (\widetilde{\mathscr{P}}_h z)^l$$. We note here that originally in Lubich & Mansour (2015, Definition 8.1) an extra term appeared involving $$\partial^{\bullet} z$$ and the surface velocity, which is not needed for the parabolic case. The Ritz map above is still well defined. Galerkin orthogonality does not hold in this case, just up to a small defect. Lemma 6.2 (Galerkin orthogonality up to a small defect) For any $$z\in H^{1}({\it{\Gamma}}(t))$$ and for all $$\varphi_{h}\in (S_h^k(t))^l$$ and for a sufficiently small $$h$$,   |a∗(z−Phz,φh)|≤chk+1‖Phz‖H1(Γ(t))‖φh‖H1(Γ(t)), (6.2) where $$c$$ is independent of $$\xi$$, $$h$$ and $$t$$. Proof. Using the definition of the Ritz map and Lemma 5.6, we estimate   |a∗(z−Phz,φh)| =|ah∗(P~hz,ϕh)−a∗(Phz,φh)|≤chk+1‖Phz‖H1(Γ(t))‖φh‖H1(Γ(t)). □ 6.1 Errors in the Ritz map Now we prove higher-order error estimates for the Ritz map (6.1) and also for its material derivatives; the analogous results for the first-order ESFEM case can be found in Dziuk & Elliott (2013b, Section 6) or in Mansour (2013, Section 7). 6.1.1 Error bounds for the Ritz map Theorem 6.3 Let $$z:\mathscr{G}_T \rightarrow \mathbb{R}$$ with $$z\in H^{k+1}({\it{\Gamma}}(t))$$ for every $$0 \leq t \leq T$$. Then the error in the Ritz map satisfies the bound, for $$0 \leq t \leq T$$ and for $$h \leq h_0$$ with a sufficiently small $$h_0$$,   ‖z−Phz‖L2(Γ(t))+h‖z−Phz‖H1(Γ(t))≤chk+1‖z‖Hk+1(Γ(t)), where the constant $$c>0$$ is independent of $$h$$ and $$t$$. Proof. To ease the presentation, we suppress all time arguments $$t$$ appearing in the norms within the proof (except for some special occasions). (a) We first prove the gradient estimate. Starting with the definition of the $$H^1({\it{\Gamma}}(t))$$-norm, then using the estimate (6.2), we have   ‖z−Phz‖H1(Γ)2 =a∗(z−Phz,z−Phz) =a∗(z−Phz,z−Ihz)+a∗(z−Phz,Ihz−Phz) ≤‖z−Phz‖H1(Γ)‖z−Ihz‖H1(Γ)+chk+1‖Phz‖H1(Γ)‖Ihz−Phz‖H1(Γ) ≤chk‖z−Phz‖H1(Γ)‖z‖Hk+1(Γ) +chk+1(2‖z−Phz‖H1(Γ)2+‖z‖H1(Γ)2+ch2k‖z‖Hk+1(Γ)2), using interpolation error estimate, and for the second term, we used the estimate   ‖Phz‖H1(Γ)‖Ihz−Phz‖H1(Γ) ≤(‖Phz−z‖H1(Γ)+‖z‖H1(Γ))(‖Ihz−z‖H1(Γ)+‖z−Phz‖H1(Γ)) ≤2‖z−Phz‖H1(Γ)2+‖z‖Hk+1(Γ)2+ch2k‖z‖Hk+1(Γ)2. Now using Young’s inequality, and for a sufficiently small (but $$t$$ independent) $$h\leq h_0$$, we have the gradient estimate   ‖z−Phz‖H1(Γ(t))≤chk‖z‖Hk+1(Γ(t)). (b) The $$L^2$$-estimate follows from the Aubin–Nitsche trick. Let us consider the problem   −ΔΓ(t)w+w=z−PhzonΓ(t); then the usual elliptic theory (see, e.g., Dziuk & Elliott, 2013a, Section 3.1; Aubin, 1998) yields the following: the solution $$w\in H^2({\it{\Gamma}}(t))$$ satisfies the bound, with $$c>0$$ independent of $$t$$,   ‖w‖H2(Γ(t))≤c‖z−Phz‖L2(Γ(t)). By testing the elliptic weak problem with $$z-\mathscr{P}_h z$$, using (6.2) again, and using the linear finite element interpolation $$I_h^{(1)}$$ on $${\it{\Gamma}}_h^k(t)$$, we obtain   ‖z−Phz‖L2(Γ)2 =a∗(z−Phz,w) =a∗(z−Phz,w−Ih(1)w)+a∗(z−Phz,Ih(1)w) ≤‖z−Phz‖H1(Γ)‖w−Ih(1)w‖H1(Γ)+chk+1‖Phz‖H1(Γ)‖Ih(1)w‖H1(Γ) ≤chk‖z‖Hk+1(Γ)ch‖w‖H2(Γ)+chk+1‖Phz‖H1(Γ)‖Ih(1)w‖H1(Γ), where for the second term we have now used   ‖Phz‖H1(Γ)‖Ih(1)w‖H1(Γ) ≤(‖z−Phz‖H1(Γ)+‖z‖H1(Γ))(‖w−Ih(1)w‖H1(Γ)+‖w‖H1(Γ)) ≤(1+chk)‖z‖Hk+1(Γ)(1+ch)‖w‖H2(Γ). Then the combination of the gradient estimate for the Ritz map and the interpolation error yields   ‖z−Phz‖L2(Γ(t))1c‖w‖H2(Γ(t))≤‖z−Phz‖L2(Γ(t))2≤chk+1‖z‖Hk+1(Γ(t))‖w‖H2(Γ(t)), which completes the proof. □ 6.1.2 Error bounds for the material derivatives of the Ritz map Since, in general, $$\partial^{\bullet}_h \mathscr{P}_h z = \mathscr{P}_h \partial^{\bullet}_h z$$ does not hold, we need the following result. Theorem 6.4 The error in the material derivatives of the Ritz map, for any $$\ell\in\mathbb{N}$$, satisfies the following bounds, for $$0 \leq t \leq T$$ and for $$h \leq h_0$$ with a sufficiently small $$h_0$$:   ‖(∂h∙)(ℓ)(z−Phz)‖L2(Γ(t))+h ‖∇Γ(∂h∙)(ℓ)(z−Phz)‖L2(Γ(t))≤cℓhk+1∑j=0ℓ‖(∂∙)(j)z‖Hk+1(Γ(t)), where the constant $$c_\ell>0$$ is independent of $$h$$ and $$t$$. Proof. The proof is a modification of Mansour (2013, Theorem 7.3). Again, to ease the presentation, we suppress the $$t$$ argument of the surfaces norms (except for some special occasions). For $$\ell=1$$: (a) We start by taking the time derivative of the definition of the Ritz map (6.1), use the transport properties Lemma 5.5 and apply the definition of the Ritz map once more; we arrive at   a∗(∂h∙z,φh) =−b(vh;z,φh)−g(vh;z,φh) +ah∗(∂h∙P~hz,ϕh)+bh(Vh;P~hz,ϕh)+gh(Vh;P~hz,ϕh). Then we obtain   a∗(∂h∙z−∂h∙Phz,φh) =−b(vh;z−Phz,φh)−g(vh;z−Phz,φh) +F1(φh), (6.3) where   F1(φh) =(ah∗(∂h∙P~hz,ϕh)−a∗(∂h∙Phz,φh)) +(bh(Vh;P~hz,ϕh)−b(vh;Phz,φh)) +(gh(Vh;P~hz,ϕh)−g(vh;Phz,φh)). Using the geometric estimates of Lemma 5.6, $$F_1$$ can be estimated as   |F1(φh)|≤chk+1(‖∂h∙Phz‖H1(Γ(t))+‖Phz‖H1(Γ(t)))‖φh‖H1(Γ(t)). (6.4) The velocity error estimate Lemma 5.4 yields   ‖∂h∙z‖H1(Γ)≤‖∂∙z‖H1(Γ)+chk‖z‖H2(Γ). Then using $$\partial^{\bullet}_h \mathscr{P}_h z$$ as a test function in (6.3), and using the error estimates of the Ritz map, together with the estimates above, with $$h\leq h_0$$, we have   ‖∂h∙Phz‖H1(Γ)2 = a∗(∂h∙Phz,∂h∙Phz) =b(vh;z−Phz,∂h∙Phz)+g(vh;z−Phz,∂h∙Phz)+a∗(∂h∙z,∂h∙Phz)−F1(∂h∙Phz) ≤ chk‖z‖Hk+1(Γ)‖∂h∙Phz‖H1(Γ)+‖∂h∙z‖H1(Γ)‖∂h∙Phz‖H1(Γ) +chk+1(‖∂h∙Phz‖H1(Γ)+‖Phz‖H1(Γ))‖∂h∙Phz‖H1(Γ) ≤chk‖z‖Hk+1(Γ)‖∂h∙Phz‖H1(Γ)+(‖∂∙z‖H1(Γ)+chk‖z‖H2(Γ))‖∂h∙Phz‖H1(Γ) +chk+1(‖∂h∙Phz‖H1(Γ)+‖z−Phz‖H1(Γ)+‖z‖H1(Γ))‖∂h∙Phz‖H1(Γ) ≤(chk‖z‖Hk+1(Γ)+‖∂∙z‖H1(Γ))‖∂h∙Phz‖H1(Γ)+chk+1‖∂h∙Phz‖H1(Γ)2, absorption using an $$h\leq h_0$$, with a sufficiently small $$h_0>0$$, and dividing through yields   ‖∂h∙Phz‖H1(Γ)≤c‖∂∙z‖H1(Γ)+chk‖z‖Hk+1(Γ). Combining all the previous estimates and using Young’s inequality, the Cauchy–Schwarz inequality and Theorem 6.3, for a sufficiently small $$h\leq h_0$$, we obtain   a∗(∂h∙z−∂h∙Phz,φh) ≤ c‖z−Phz‖H1(Γ)‖φh‖H1(Γ) +chk+1(‖∂h∙Phz‖H1(Γ)+‖Phz‖H1(Γ))‖φh‖H1(Γ) ≤ chk‖z‖Hk+1(Γ)‖φh‖H1(Γ) +chk+1(‖∂∙z‖H1(Γ)+(1+chk)‖z‖Hk+1(Γ))‖φh‖H1(Γ) ≤ chk(‖z‖Hk+1(Γ)+h‖∂∙z‖H1(Γ))‖φh‖H1(Γ). Then, similarly to the previous proof, we have   ‖∂h∙z−∂h∙Phz‖H1(Γ)2 ≤a∗(∂h∙z−∂h∙Phz,∂h∙z−∂h∙Phz) =a∗(∂h∙z−∂h∙Phz,∂h∙z−Ih∂∙z)+a∗(∂h∙z−∂h∙Phz,Ih∂∙z−∂h∙Phz) ≤‖∂h∙z−∂h∙Phz‖H1(Γ)‖∂h∙z−Ih∂∙z‖H1(Γ) +chk(‖z‖Hk+1(Γ)+h‖∂∙z‖H1(Γ))‖Ih∂∙z−∂h∙Phz‖H1(Γ). Finally, the interpolation estimates, Young’s inequality and absorption using a sufficiently small $$h\leq h_0$$, yields the gradient estimate. (b) The $$L^2$$ estimate again follows from the Aubin–Nitsche trick. Let us now consider the problem   −ΔΓ(t)w+w=∂h∙z−∂h∙PhzonΓ(t), together with the usual elliptic estimate, for the solution $$w\in H^2({\it{\Gamma}}(t))$$,   ‖w‖H2(Γ(t))≤c‖∂h∙z−∂h∙Phz‖L2(Γ(t)); again, $$c$$ is independent of $$t$$ and $$h$$. Following the proof of Dziuk & Elliott (2013b, Theorem 6.2), let us first bound   −b(vh;z−Phz,Ih(1)w) =b(vh;z−Phz,w−Ih(1)w)−b(vh;z−Phz,w) ≤chk‖z‖Hk+1(Γ) ch‖w‖H2(Γ)−b(vh;z−Phz,w) =chk+1‖z‖Hk+1(Γ)‖w‖H2(Γ)+b(v;z−Phz,w) +b(v;z−Phz,w)−b(vh;z−Phz,w), where again $$I_h^{(1)}$$ denotes the linear finite element interpolation operator on $${\it{\Gamma}}_h^k(t)$$. The pair in the last line can be estimated, using Lemma 5.4, by   b(v;z−Phz,w)−b(vh;z−Phz,w)≤ ∫Γ(t)|B(v)−B(vh)||∇Γ(z−Phz)||∇Γw|≤ ch2k‖z‖Hk+1(Γ)‖w‖H1(Γ). Finally, for the remaining term, the proof of Dziuk & Elliott (2013b, Theorem 6.2) yields   b(v;z−Phz,w)≥−c‖z−Phz‖L2(Γ)‖w‖H2(Γ)≥−chk+1‖z‖Hk+1(Γ)‖w‖H2(Γ). For the other bilinear form in (6.3), we have   g(vh;z−Phz,φh)≤chk+1‖z‖Hk+1(Γ)‖w‖H1(Γ). The combination of all these estimates with (6.4) yields   a∗(∂h∙z−∂h∙Phz,Ihw)≤chk+1‖z‖Hk+1(Γ)‖w‖H2(Γ). By testing the above elliptic weak problem with $$z-\mathscr{P}_h z$$, and using the above bound and the gradient estimate from (a), we obtain   ‖∂h∙z−∂h∙Phz‖L2(Γ)2= a∗(∂h∙z−∂h∙Phz,w)= a∗(∂h∙z−∂h∙Phz,w−Ih(1)w)+a∗(∂h∙z−∂h∙Phz,Ih(1)w)≤ chk(‖z‖Hk+1(Γ)+h‖∂∙z‖H1(Γ)) ch‖w‖H2(Γ(t))+chk+1‖z‖Hk+1(Γ)‖w‖H2(Γ)≤ chk+1(‖z‖Hk+1(Γ)+h‖∂∙z‖H1(Γ))‖w‖H2(Γ). For $$\ell>1,$$ the proof is analogous. □ 7. Error bounds for the semidiscretization and full discretization 7.1 Convergence proof for the semidiscretization By combining the error estimates in the Ritz map and in its material derivatives and the geometric results of Section 5, we prove convergence of the high-order ESFEM semidiscretization. Proof of Theorem 4.1. The result is simply shown by repeating the arguments of Dziuk & Elliott (2013b, Section 7) for our setting but using the high-order versions for all results: geometric estimates Lemma 5.2, perturbation estimates of bilinear forms Lemma 5.6 and Ritz map error estimates Theorems 6.3 and 6.4. □ 7.2 Convergence proof for the full discretization 7.2.1 Bound of the semidiscrete residual We follow the approach of Mansour (2013, Section 8.1) and Lubich et al. (2013, Section 5) by defining the ESFEM residual $$R_h(\cdot,t) = \sum_{j=1}^N r_j(t) \chi_j(\cdot,t)\in S_h^k(t)$$ as   ∫Γhk(t)Rhϕh=ddt∫Γhk(t)P~huϕh+∫Γhk(t)∇Γh(P~hu)⋅∇Γhϕh−∫Γhk(t)(P~hu)∂h∙ϕh, (7.1) where $$\phi_h\in S_h^k(t)$$, and $$\widetilde{\mathscr{P}}_h u(\cdot,t)$$ is the Ritz map of the smooth solution $$u$$. We now show the optimal-order $$H_h^{-1}$$-norm estimate of the residual $$R_h$$. Theorem 7.1 Let the solution $$u$$ of the parabolic problem be sufficiently smooth. Then there exist $$C>0$$ and $$h_0>0$$, such that for all $$h\leq h_0$$ and $$t\in[0,T]$$, the finite element residual $$R_h$$ of the Ritz map is bounded as   ‖Rh‖Hh−1(Γh(t))≤Chk+1. Proof. (a) We start by applying the discrete transport property to the residual equation (7.1):   mh(Rh,ϕh)= ddtmh(P~hu,ϕh)+ah(P~hu,ϕh)−mh(P~hu,∂h∙ϕh)= mh(∂h∙P~hu,ϕh)+ah(P~hu,ϕh)+gh(Vh;P~hu,ϕh). (b) We continue by the transport property with discrete material derivatives from Lemma 5.5 but for the weak form, with $$\varphi:=\varphi_h=(\phi_h)^l$$:   0= ddtm(u,φh)+a(u,φh)−m(u,∂∙φh)= m(∂h∙u,φh)+a(u,φh)+g(vh;u,φh)+m(u,∂h∙φh−∂∙φh). (c) Subtraction of the two equations, using the definition of the Ritz map (6.1) and using that   ∂h∙φh−∂∙φh=(vh−v)⋅∇Γφh holds, we obtain   mh(Rh,ϕh) =mh(∂h∙P~hu,ϕh)−m(∂h∙u,φh) +gh(Vh;P~hu,ϕh)−g(vh;u,φh) +m(u,φh)−mh(P~hu,ϕh) +m(u,(vh−v)⋅∇Γφh). All the pairs can be easily estimated separately as $$c h^{k+1} \|\varphi_h\|_{H^1({\it{\Gamma}}(t))}$$ by combining the geometric perturbation estimates of Lemma 5.6, the velocity estimate of Lemma 5.4 and the error estimates of the Ritz map from Theorems 6.3 and 6.4. The proof is finished using the definition of the $$H_h^{-1}$$-norm and the equivalence of norms Lemma 5.1. □ 7.2.2 Proof of Theorem 4.4. Using the error estimate for the BDF methods Theorem 4.2 and using the bounds for the semidiscrete residual Theorem 7.1, we give here a proof for the fully discrete error estimates of Theorem 4.4. Proof of Theorem 4.4. The global error is decomposed into two parts,   uhn−u(⋅,tn)=(uhn−(Phu)(⋅,tn))+((Phu)(⋅,tn)−u(⋅,tn)), and then the terms are estimated separately by results from above. The first term is estimated, analogously to Thomée (2006) or exactly as in Lubich et al. (2013) and Mansour (2013) as follows. The vectors collecting the nodal values of the error $$u_h^n - (\mathscr{P}_h u)(\cdot,t_n)$$ satisfy the fully discrete problem (4.1) perturbed by the semidiscrete residual (7.1) (cf. Mansour, 2013). Then applying results for BDF methods Theorem 4.2, together with the residual bound Theorem 7.1 (and by the assumption on the initial value approximation), gives the desired bound $${\mathscr {O}}(h^{k+1} + \tau^p)$$. The second term is directly estimated by the error estimates for the Ritz map and for its material derivatives, Theorems 6.3 and 6.4. □ Remark 7.2 In Remark 4.5, we noted that the analogous fully discrete results can be shown for algebraically stable implicit Runge–Kutta methods. To be more precise, for the Runge–Kutta analogue, instead of Theorem 4.2 one has to use the error bounds of Dziuk et al. (2012, Theorem 5.1), but otherwise the proof above remains the same. Remark 7.3 The statement in the introduction on various extensions is supported by the following facts. The first two points in the introduction, geometric errors and perturbation errors of the bilinear forms (Sections 5.1 and 5.4), and the basic high-order evolving surface finite element setting (Section 3.1), are independent of the considered problem. The velocity estimates (Section 5.3) depend only on the surface evolution as well. Furthermore, these general high-order results could be used in analogous ways, as their linear counterparts in the cited articles; or missing ones can be easily shown based on the ideas presented here (e.g., involving the additional term for ALE maps (Elliott & Venkataraman, 2014; Kovács & Power Guerra, 2014). The proof of the error bounds of a modified Ritz map (although perhaps defined slightly differently; see, e.g., wave equations Lubich & Mansour, 2015, or quasilinear problems Kovács & Power Guerra, 2016), rely on these geometric approximation results and some general techniques, such as the almost Galerkin orthogonality and the Aubin–Nitsche trick. Naturally, the estimates for the semidiscrete residual greatly depend on the problem itself, but in the above-mentioned articles, similar (if not the same) ideas and techniques were used. Finally, convergence of time discretizations, proved by energy estimates, carry over to high-order discretizations of various problems, from Dziuk et al. (2012), Lubich et al. (2013), Kovács & Power Guerra (2014), Lubich & Mansour (2015) and Kovács & Power Guerra (2016). Remark 7.4 Theorem 4.2 requires sufficient temporal regularity of the Ritz map. By having sufficient regularity of the solution and the surface evolution on $$[0,T]$$, using Theorem 6.4 and equivalence of norms, we obtain   ‖(∂h∙)(ℓ)(P~hu)(⋅,t)‖L2(Γh(t))≤ c‖(∂∙)(ℓ)(u−Phu)(⋅,t)‖L2(Γ(t))+c‖(∂∙)(ℓ)u‖L2(Γ(t))≤ c(cℓhk+1+1)∑j=0ℓ‖(∂∙)(j)u‖Hk+1(Γ(t)). Here, the $$H^{k+1}({\it{\Gamma}}(t))$$-norm on the right-hand side could be replaced by $$H^{2}({\it{\Gamma}}(t))$$ (also the power in $$h^{k+1}$$ would be reduced to $$2$$) by modifiying the proofs of Theorems 6.3 and 6.4, using the linear interpolation $$I_h^{(1)}$$ on $${\it{\Gamma}}_h^k(t)$$ instead of $$I_h$$. Alternatively, a weaker condition can be obtained in the following way. By having sufficient regularity of the solution at $$t=0$$ and having smooth evolution of the surface, by repeating the proof of Dziuk et al. (2012, Theorem 9.1) for the Ritz map instead of the semidiscrete solution, the following estimate is obtained:   supt∈[0,T]‖(∂h∙)(ℓ)P~hu(⋅,t)‖L2(Γhk(t))+∫0T‖∇Γh((∂h∙)(ℓ)P~hu(⋅,s))‖L2(Γhk(s))ds≤c∑j=1ℓ‖(∂h∙)(j)u(⋅,0)‖L2(Γhk(0)). The result could even be obtained directly by using Dziuk et al. (2012, Theorem 9.1) and modifying our results using the semidiscrete solution instead of the Ritz map. 8. Implementation The implementation of the high-order ESFEM code follows the typical method of finite element mass matrix, stiffness matrix and load vector assembly, mixed with techniques from isoparametric FEM theory. This was also used for linear ESFEM (see Dziuk & Elliott, 2007, Section 7.2). Similarly to the linear case, a curved element $$E_h$$ of the $$k$$-order interpolation surface $${\it{\Gamma}}_h^k(t)$$ is parameterized over the reference triangle $$E_0$$, chosen to be the unit simplex. Then the polynomial map of degree $$k$$ between $$E_h$$ and $$E_0$$ is used to compute the local matrices. All the computations are done on the reference element, using the Dunavant quadrature rule (see Dunavant, 1985). Then the local values are summed to their correct places in the global matrices. In a typical case, the surface is evolved by solving a series of ODEs, hence only the initial mesh is created based on $${\it{\Gamma}}(0)$$. Naturally, the problem of velocity-based grid distortion is still present. Possible ways to overcome this are methods using the DeTurck trick (see Elliott & Fritz, 2016) or using ALE finite elements (see Elliott & Venkataraman, 2014; Kovács & Power Guerra, 2014). 9. Numerical experiments We performed various numerical experiments with quadratic approximation of the surface $${\it{\Gamma}}(t)$$ and using quadratic ESFEM to illustrate our theoretical results. 9.1 Example 1: Parabolic problem on a stationary surface Let us briefly report on numerical experiments for parabolic problems on a stationary surface, as a benchmark problem. Let $${\it{\Gamma}}\subset\mathbb{R}^3$$ be the unit sphere, and let us consider the parabolic surface PDE   ∂tu−ΔΓu=f, with given initial value and inhomogeneity $$f$$ chosen such that the solution is $$u(x,t)=e^{-6t}x_1x_2$$. Let $$(\mathscr{T}_k)_{k=1,2,\dotsc,n}$$ and $$(\tau_k)_{k=1,2,\dotsc,n}$$ be series of meshes and time steps, respectively, such that $$2 h_k = h_{k-1}$$ and $$2 \tau_k = \tau_{k-1}$$, with $$h_1=\sqrt2$$ and $$\tau_1=0.2$$. For each mesh $$\mathscr{T}_k$$ with corresponding step size $$\tau_k$$, we numerically solve the surface PDE using second-order ESFEM combined with third-order BDF methods. Then, by $$e_k$$ we denote the error corresponding to the mesh $$\mathscr{T}_k$$ and step size $$\tau_k$$. We then compute the errors between the lifted numerical solution and the exact solution using the following norm and seminorm:   L∞(L2): max1≤n≤N‖uhn−u(⋅,tn)‖L2(Γ(tn)),L2(H1): (τ∑n=1N‖∇Γ(tn)(uhn−u(⋅,tn))‖L2(Γ(tn))2)1/2. Using the above norms, the experimental order of convergence rates (EOCs) are computed by   EOCk=ln⁡(ek/ek−1)ln⁡(2)(k=2,3,…,n). In Table 1, we report on the EOCs for the second-order ESFEM coupled with the BDF3 method; theoretically, we expect EOC$$\approx3$$ in the $$L^\infty(L^2)$$-norm, and EOC$$\approx2$$ in the $$L^2(H^1)$$-seminorm. Table 1. Errors and EOCs in the$$L^\infty(L^2)$$- and$$L^2(H^1)$$-norms for the stationary problem Level  dof  $$L^\infty(L^2)$$  EOC  $$L^2(H^1)$$  EOC  $$1$$  $$6$$  $$5.3113\cdot 10^{-3}$$  —  $$8.0694\cdot 10^{-3}$$  —  $$2$$  $$18$$  $$2.9257\cdot 10^{-3}$$  $$0.9511$$  $$4.3162\cdot 10^{-3}$$  $$0.99805$$  $$3$$  $$66$$  $$9.2303\cdot 10^{-4}$$  $$1.7122$$  $$1.9050\cdot 10^{-3}$$  $$1.2139$$  $$4$$  $$258$$  $$1.7285\cdot 10^{-4}$$  $$2.4338$$  $$5.4403\cdot 10^{-4}$$  $$1.8207$$  $$5$$  $$1026$$  $$2.6463\cdot 10^{-5}$$  $$2.7124$$  $$1.2265\cdot 10^{-4}$$  $$2.1529$$  $$6$$  $$4098$$  $$3.6845\cdot 10^{-6}$$  $$2.8457$$  $$2.4772\cdot 10^{-5}$$  $$2.3088$$  Level  dof  $$L^\infty(L^2)$$  EOC  $$L^2(H^1)$$  EOC  $$1$$  $$6$$  $$5.3113\cdot 10^{-3}$$  —  $$8.0694\cdot 10^{-3}$$  —  $$2$$  $$18$$  $$2.9257\cdot 10^{-3}$$  $$0.9511$$  $$4.3162\cdot 10^{-3}$$  $$0.99805$$  $$3$$  $$66$$  $$9.2303\cdot 10^{-4}$$  $$1.7122$$  $$1.9050\cdot 10^{-3}$$  $$1.2139$$  $$4$$  $$258$$  $$1.7285\cdot 10^{-4}$$  $$2.4338$$  $$5.4403\cdot 10^{-4}$$  $$1.8207$$  $$5$$  $$1026$$  $$2.6463\cdot 10^{-5}$$  $$2.7124$$  $$1.2265\cdot 10^{-4}$$  $$2.1529$$  $$6$$  $$4098$$  $$3.6845\cdot 10^{-6}$$  $$2.8457$$  $$2.4772\cdot 10^{-5}$$  $$2.3088$$  In Figs 1 and 2, we report on the errors   ‖u(⋅,Nτ)−uhN‖L2(Γ)and‖∇Γ(u(⋅,Nτ)−uhN)‖L2(Γ), at time $$N\tau=1$$. The logarithmic plots show the errors against the mesh width $$h$$ (in Fig. 1) and against time-step size $$\tau$$ (in Fig. 2). Fig. 1. View largeDownload slide Spatial convergence of the BDF3 / quadratic SFEM discretization for the stationary surface PDE. Fig. 1. View largeDownload slide Spatial convergence of the BDF3 / quadratic SFEM discretization for the stationary surface PDE. Fig. 2. View largeDownload slide Temporal convergence of the BDF3/quadratic SFEM discretization for the stationary surface PDE. Fig. 2. View largeDownload slide Temporal convergence of the BDF3/quadratic SFEM discretization for the stationary surface PDE. The different lines correspond to different time-step sizes and to different mesh refinements, respectively in Figs 1 and 2. In both figures we can observe two regions. In Fig. 1, a region where the spatial discretization error dominates, matching the convergence rates of our theoretical results, and a region, with fine meshes, where the time discretization error dominates (the error curves flatten out). In Fig. 2, the same description applies but with reversed roles. First, the time discretization error dominates, while for smaller step sizes the spatial error dominates. The convergence in space (Fig. 1) and in time (Fig. 2) can both be observed to be nicely in agreement with the theoretical results (note the reference lines). 9.2 Example 2: Evolving surface parabolic problem In the following experiment, we consider the parabolic problem (2.1) on the evolving surface given by   Γ(t)={x∈R3 | a(t)−1x12+x22+x32−1=0}, where $$a(t)=1+\frac{1}{4}\sin(2\pi t)$$ (see, e.g., Dziuk & Elliott, 2007; Dziuk et al., 2012; Mansour, 2013), with given initial value and inhomogeneity $$f$$ chosen such that the solution is $$u(x,t)=e^{-6t}x_1x_2$$. Similarly to the stationary surface case, we again report on the experimental orders of convergence and similar spatial and temporal convergence plots. They are all produced exactly as described above. The EOCs for the evolving surface problem solved with BDF method of order $$3$$ and evolving surface finite elements of second order can be seen in Table 2. Table 2. Errors and EOCs in the$$L^\infty(L^2)$$- and$$L^2(H^1)$$-norms for the evolving surface problem Level  dof  $$L^\infty(L^2)$$  EOC  $$L^2(H^1)$$  EOC  1  6  $$5.7898\cdot 10^{-3}$$  —  $$8.4446\cdot 10^{-3}$$  —  2  18  $$9.5840\cdot 10^{-4}$$  2.8688  $$1.8173\cdot 10^{-3}$$  2.4504  3  66  $$4.0725\cdot 10^{-4}$$  1.2704  $$1.6549\cdot 10^{-3}$$  0.13895  4  258  $$9.1096\cdot 10^{-5}$$  2.1755  $$5.4513\cdot 10^{-4}$$  1.6133  5  1026  $$1.3847\cdot 10^{-5}$$  2.7226  $$1.2774\cdot 10^{-4}$$  2.0971  6  4098  $$1.9534\cdot 10^{-6}$$  2.8267  $$2.6135\cdot 10^{-5}$$  2.2901  Level  dof  $$L^\infty(L^2)$$  EOC  $$L^2(H^1)$$  EOC  1  6  $$5.7898\cdot 10^{-3}$$  —  $$8.4446\cdot 10^{-3}$$  —  2  18  $$9.5840\cdot 10^{-4}$$  2.8688  $$1.8173\cdot 10^{-3}$$  2.4504  3  66  $$4.0725\cdot 10^{-4}$$  1.2704  $$1.6549\cdot 10^{-3}$$  0.13895  4  258  $$9.1096\cdot 10^{-5}$$  2.1755  $$5.4513\cdot 10^{-4}$$  1.6133  5  1026  $$1.3847\cdot 10^{-5}$$  2.7226  $$1.2774\cdot 10^{-4}$$  2.0971  6  4098  $$1.9534\cdot 10^{-6}$$  2.8267  $$2.6135\cdot 10^{-5}$$  2.2901  The errors at time $$N\tau=1$$ in different norms can be seen in the following plots: the different lines again correspond to different time-step sizes and to different mesh refinements in Figs 3 and 4, respectively. Fig. 3. View largeDownload slide Spatial convergence of the BDF3/quadratic ESFEM discretization for the evolving surface PDE. Fig. 3. View largeDownload slide Spatial convergence of the BDF3/quadratic ESFEM discretization for the evolving surface PDE. Fig. 4. View largeDownload slide Temporal convergence of the BDF3/quadratic ESFEM discretization for the evolving surface PDE. Fig. 4. View largeDownload slide Temporal convergence of the BDF3/quadratic ESFEM discretization for the evolving surface PDE. Acknowledgements The author is grateful for the valuable discussions with Prof. Christian Lubich during the preparation of the manuscript. We would also like to thank the anonymous referees for their extremely constructive suggestions and comments. Funding Deutsche Forschungsgemeinschaft, SFB 1173 to B.K. Appendix. Proof of the geometric approximation results of Lemma 5.2 For clarity we recall Lemma 5.2. Lemma A.1 For $${\it{\Gamma}}_h^k(t)$$ and $${\it{\Gamma}}(t)$$ as above, for $$h\leq h_0$$ with a sufficiently small $$h_0>0$$, we have the geometric approximation estimates:    ‖d‖L∞(Γh(t))≤chk+1, ‖1−δh‖L∞(Γh(t))≤chk+1, ‖n−nh‖L∞(Γh(t))≤chk, and ‖Id−δhQh‖L∞(Γh(t))≤chk+1, ‖(∂h∙)(ℓ)d‖L∞(Γh(t))≤chk+1, ‖(∂h∙)(ℓ)δh‖L∞(Γh(t))≤chk+1, ‖Pr((∂h∙)(ℓ)Qh)Pr‖L∞(Γh)≤chk+1, with constants depending only on $$\mathscr{G}_T$$, but not on $$h$$ or $$t$$. Proof of Lemma 5.2. The proofs follow the ideas of Dziuk & Elliott (2007, Lemma 5.1) and Mansour (2013, Lemma 6.1), in combination with the ideas and techniques of the proof of Demlow (2009, Proposition 2.3). Let $$E=E(t) \subset {\it{\Gamma}}_h^k(t)$$ be an element of the discrete surface. By $$I_h$$ we denote the $$k$$-order interpolation operator of Section 5.2. (i) Since the nodes of $$E$$ lie on the exact surface $${\it{\Gamma}}(t)$$ we have that the interpolate $$\widetilde I_h d$$ vanishes on $$E$$. Then by using standard interpolation estimates (from Lemma 5.3 or from Brenner & Scott, 2007) we obtain   ‖d‖L∞(E)=‖d−I~hd‖L∞(E)≤chk+1‖d‖Wk+1,∞(E)≤chk+1‖d‖Ck,1(GT). Higher-order norm estimates are shown analogously:   ‖p−pk‖Wi,∞≤chk+1−i. (A.1) (iii) For the normal vector estimate,   |n(x^)−nhk(x^)|≤ |n(pk(x~))−n(p(x~))|+|n(p(x~))−nhk(pk(x~))|≤ c(GT)hk+1+c(GT)hk, (A.2) where for the last estimate we used the smoothness of $$\mathscr{G}_T$$ and the above bound on $$d$$, and the bounds (A.1) and the Gram–Schmidt orthonormalization algorithm (cf. Demlow, 2009). (ii) The second estimate is shown by recalling, from Demlow (2009, (2.10)), that (for a fixed $$t\in[0,T]$$)   δhk(x)=n(x)⋅nhk(x)∏j=1m(1−d(x,t)Kj(x))(x∈Γhk(t)), (A.3) where $$K_j(x) = {\kappa_j(p(x,t))}/{1+d(x,t)\kappa_j(p(x,t))}$$ with $$\kappa_j$$ being the principle curvatures; cf. Demlow (2009) Then, following the proof of Demlow (2009, Proposition 4.1), using $$\|d\|_{L^\infty} = {\mathscr {O}}(h^{k+1})$$ and (A.2), we obtain   |1−δhk|≤ c(GT)hk+1+c(GT)|1−n⋅nhk|≤ c(GT)hk+1+c(GT)|n−nhk|2≤c(GT)hk+1. (A.4) (iv) To show the fourth estimate, we use the idea of the linear ESFEM case. Using the previous estimates and (5.1) the definition of $$Q_h$$,   |Id−δhQh|≤|Pr−PrPrhPr|+chk+1. Then, for an arbitrary unit vector $$z$$,   |(Pr−PrPrhPr)z|=|z⋅(nh−(nh⋅n)n)(nh−(nh⋅n)n)|≤ch2k, where the estimate follows, using (A.2), from   |nh−(nh⋅n)n|≤ |(n⋅n)nh−(nh⋅n)nh|+|(nh⋅n)nh−(nh⋅n)n|≤ |((n−nh)⋅n)nh|+|(nh⋅n)(nh−n)|≤chk. See also the proof of Demlow (2009, Proposition 4.1). The proofs of the estimates with material derivatives are similar to their nondifferentiated versions; we follow the ideas of Mansour (2013). (v) Again, since $${\left( {\partial _h^ \bullet } \right)^{(l)}} d$$ vanishes at the nodes of $$E$$, hence the interpolant $$\widetilde I_h {\left( {\partial _h^ \bullet } \right)^{(l)}} d$$ vanishes on $$E$$ completely. Again by interpolation estimates we obtain   ‖(∂h∙)(l)d‖L∞(E)=‖(∂h∙)(l)d−I~h(∂h∙)(l)d‖L∞(E)≤chk+1‖(∂h∙)(l)d‖Wk+1,∞(E)≤chk+1‖(∂h∙)(l)d‖Ck,1(GT). (vi) The sixth estimate is shown by taking the material derivative of (A.3) and by a similar argument to (A.4),   ‖(∂h∙)(ℓ)δh‖L∞(Γh(t))≤ c(GT)hk+1+c(GT)|∂h∙(n−nhk)|2≤c(GT)hk+1, where the last inequality follows by the chain rule and since the combination of (A.1) and (A.2), together with $$\mathrm{n}=\nabla d$$, yields $$|\partial^{\bullet}_h(\mathrm{n}-\mathrm{n}_h^k)| = {\mathscr {O}}(h^k)$$. (vii) Let us show the last estimate for $$\ell=1$$. The higher-order version follows by a similar argument recursively. It is clear from the definition of $$Q_h$$ (see (5.1)) that with a smooth remainder function $$R$$ we can write   Qh=1δhPrPrhPr+dR(δh,Pr,Prh,H). Now using the facts $$\|d\| , \|\partial^{\bullet}_h d\| ={\mathscr {O}}(h^{k+1})$$, $$\delta_h=1+{\mathscr {O}}(h^{k+1})$$ and $$\|\partial^{\bullet}_h \delta_h\|={\mathscr {O}}(h^{k+1})$$ we bound   Pr(∂h∙Qh)Pr=Pr∂h∙(PrPrhPr)Pr+O(hk+1). The first term here is estimated separately. Using $$\partial^{\bullet}_h \mathrm{n} \cdot \mathrm{n} = 0$$, in Mansour (2013, equation 6.4), it is shown that   Pr∂h∙(PrPrhPr)Pr=Pr∂h∙(PrPrhPr−Pr)Pr=−Pr∂h∙(Pr nhnhTPr)Pr. Since the projections are bounded, we need only the bounds   |Pr nh|=|nh−(n⋅nh)n|≤ chk,|∂h∙(Pr nh)|=|∂h∙(nh−(n⋅nh)n)|≤ chk, where the first inequality has been shown above, while the second follows by the chain rule and by $$|\partial^{\bullet}_h(\mathrm{n}-\mathrm{n}_h^k)| = {\mathscr {O}}(h^k)$$ proved above, together with the boundedness of $$\partial^{\bullet}_h \mathrm{n}$$ and $$\partial^{\bullet}_h \mathrm{n}_h$$ (which follows from velocity approximation and the smoothness of $$\mathscr{G}_T$$ for the first term, while for the latter by the same arguments and an additional triangle inequality). □ References Antonietti P., Dedner A., Madhavan P., Stangalino S., Stinner B. & Verani M. ( 2015) High order discontinuous Galerkin methods for elliptic problems on surfaces. SIAM J. Num. Anal , 53, 1145– 1171. Google Scholar CrossRef Search ADS   Aubin T. ( 1998) Some Nonlinear Problems in Riemannian Geometry,  1st edn. Berlin: Springer. Google Scholar CrossRef Search ADS   Brenner S. & Scott R. ( 2007) The Mathematical Theory of Finite Element Methods , vol. 15. New York: Springer. Burman E., Hansbo P. & Larson M. ( 2015) A stabilized cut finite element method for partial differential equations on surfaces: The Laplace–Beltrami operator. Comput. Methods Appl. Mech. Eng. , 285, 188– 207. Google Scholar CrossRef Search ADS   Dedner A., Madhavan P. & Stinner B. ( 2013) Analysis of the discontinuous Galerkin method for elliptic problems on surfaces. IMA J. Numer. Anal. , 33, 952– 973. Google Scholar CrossRef Search ADS   Demlow A. ( 2009) Higher-order finite elment methods and pointwise error estimates for elliptic problems on surfaces. SIAM J. Numer. Anal. , 47, 805– 807. Google Scholar CrossRef Search ADS   Dunavant D. ( 1985) High degree efficient symmetrical Gaussian quadrature rules for the triangle. Int. J. Numer. Meth. Eng. , 21, 1129– 1148. Google Scholar CrossRef Search ADS   Dziuk G. ( 1988) Finite elements for the Beltrami operator on arbitrary surfaces. Partial Differential Equations and Calculus of Variations . Lecture Notes in Mathematics, 1357. Berlin: Springer. pp. 142– 155. Google Scholar CrossRef Search ADS   Dziuk G. & Elliott C. ( 2007) Finite elements on evolving surfaces. IMA J. Numer. Anal. , 27, 262– 292. Google Scholar CrossRef Search ADS   Dziuk G. & Elliott C. ( 2013a) Finite element methods for surface PDEs. Acta Numerica , 22, 289– 396. Google Scholar CrossRef Search ADS   Dziuk G. & Elliott C. ( 2013b) $$L^2$$-estimates for the evolving surface finite element method. Math. Comp. , 82, 1– 24. Google Scholar CrossRef Search ADS   Dziuk G., Lubich C. & Mansour D. ( 2012) Runge–Kutta time discretization of parabolic differential equations on evolving surfaces. IMA J. Numer. Anal. , 32, 394– 416. Google Scholar CrossRef Search ADS   Elliott C. & Fritz H. ( 2016) On approximations of the curve shortening flow and of the mean curvature flow based on the DeTurck trick. To appear i. IMA J. Numer. Anal. , doi: 10.1093/imanum/drw020. Elliott C. M. & Ranner T. ( 2013) Finite element analysis for a coupled bulk–surface partial differential equation. IMA J. Numer. Anal. , 33, 377– 402. Google Scholar CrossRef Search ADS   Elliott C. M. & Ranner T. ( 2015) Evolving surface finite element method for the Cahn–Hilliard equation. Numer. Math. , 129, 483– 534. Google Scholar CrossRef Search ADS   Elliott C. M. & Venkataraman C. ( 2014) Error analysis for an ALE evolving surface finite element method. Numer. Methods Partial Differential Equations , 31, 459– 499. Google Scholar CrossRef Search ADS   Grande J. & Reusken A. ( 2014) A higher order finite element method for partial differential equations on surfaces. SIAM J. Numer. Anal. , 54, 388– 414. Google Scholar CrossRef Search ADS   Hairer E. & Wanner G. ( 1996) Solving Ordinary Differential Equations II: Stiff and Differential–Algebraic Problems , 2nd edn. Berlin: Springer. Kovács B., Li B., Lubich Ch. & Power Guerra C. A. ( 2017) Convergence of finite elements on an evolving surfacedriven by diffusion on the surface, arXiv:1607.07170v2. Kovács B. & Power Guerra C. A. ( 2014) Higher–order time discretizations with ALE finite elements for parabolic problems on evolving surfaces. to appear i. IMA J. Numer. Anal. , arXiv:1410.0486. Kovács B. & Power Guerra C. A. ( 2016) Error analysis for full discretizations of quasilinear parabolic problems on evolving surfaces. Numer. Methods Partial Differential Equations , 32, 1200– 1231. Google Scholar CrossRef Search ADS   Lubich C. & Mansour D. ( 2015) Variational discretization of wave equations on evolving surfaces. Math. Comp. , 84, 513– 542. Google Scholar CrossRef Search ADS   Lubich C., Mansour D. & Venkataraman C. ( 2013) Backward difference time discretization of parabolic differential equations on evolving surfaces. IMA J. Numer. Anal. , 33, 1365– 1385. Google Scholar CrossRef Search ADS   Mansour D. ( 2013) Numerical analysis of partial differential equations on evolving surfaces. Ph.D. Thesis , Universität Tübingen. http://hdl.handle.net/10900/49925. Accessed 1 March 2017. Olshanskii M., Reusken A. & Grande J. ( 2009) A finite element method for elliptic equations on surfaces. SIAM J. Numer. Anal. , 47, 3339– 3358. Google Scholar CrossRef Search ADS   Reusken A. ( 2015) Analysis of trace finite element methods for surface partial differential equations. IMA J. Numer. Anal. , 35, 1568– 1590. Google Scholar CrossRef Search ADS   Thomée V. ( 2006) Galerkin Finite Element Methods for Parabolic Problems , 2nd edn. Berlin: Springer. © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

# High-order evolving surface finite element method for parabolic problems on evolving surfaces

, Volume 38 (1) – Jan 1, 2018
30 pages

/lp/ou_press/high-order-evolving-surface-finite-element-method-for-parabolic-189tcaeo9x
Publisher
Oxford University Press
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drx013
Publisher site
See Article on Publisher Site

### Abstract

Abstract High-order spatial discretizations and full discretizations of parabolic partial differential equations on evolving surfaces are studied. We prove convergence of the high-order evolving surface finite element method by showing high-order versions of geometric approximation errors and perturbation error estimates and by the careful error analysis of a modified Ritz map. Furthermore, convergence of full discretizations using backward difference formulae and implicit Runge–Kutta methods are also shown. 1. Introduction Numerical methods for partial differential equations (PDEs) on stationary and evolving surfaces and for coupled bulk–surface PDEs have been under intensive research in recent years. Surface finite element methods are all based on the fundamental article of Dziuk (1988), further developed for evolving surface parabolic problems by Dziuk & Elliott (2007, 2013b). High-order versions of various finite element methods for problems on a stationary surface have received attention in a number of publications previously. We give a brief overview of this literature here: The surface finite element method of Dziuk (1988) was extended to higher-order finite elements on stationary surfaces by Demlow (2009) Some further important results for higher-order surface finite elements were shown by Elliott & Ranner (2013). Discontinuous Galerkin methods for elliptic surface problems were analysed by Dedner et al. (2013) and then extended to high-order discontinuous Galerkin methods in Antonietti et al. (2015). Recently, unfitted (also called trace or cut) finite element methods have been investigated intensively (see, e.g., Olshanskii et al., 2009; Burman et al., 2015; Reusken, 2015). A higher-order version of the trace finite element method was analysed by Grande & Reusken (2014). However, to our knowledge, there are no articles showing convergence of the high-order evolving surface finite element method for parabolic partial differential equations on evolving surfaces. In this article, we extend the $$H^1$$- and $$L^2$$-norm convergence results of Dziuk & Elliott (2007, 2013b) to high-order evolving surface finite elements applied to parabolic problems on evolving surfaces with prescribed velocity. Furthermore, convergence results for full discretizations using Runge–Kutta methods, based on Dziuk et al. (2012), and using backward difference formulae (BDF), based on Lubich et al. (2013), are also shown. To prove high-order convergence of the spatial discretization, three main groups of errors have to be analysed: Geometric errors, resulting from the appropriate approximation of the smooth surface. Many of these results carry over from Demlow (2009) by careful investigation of time dependencies, while others are extended from Dziuk et al. (2012), Mansour (2013) and Lubich et al. (2013). Perturbation errors of the bilinear forms, whose higher-order version can be shown by carefully using the core ideas of Dziuk & Elliott (2013b). High-order estimates for the errors of a modified Ritz map, which was defined in Lubich & Mansour (2015). These projection error bounds rely on the nontrivial combination of the mentioned geometric error bounds and on the well-known Aubin–Nitsche trick. We further show convergence results for full discretizations using Runge–Kutta and BDF methods. The error estimates, based on energy estimates, for Runge–Kutta methods shown in Dziuk et al. (2012) and for BDF methods in Lubich et al. (2013) are applicable without any modifications, because the semidiscrete problem can be written in a matrix–vector formulation (cf. Dziuk et al., 2012), where the matrices have exactly the same properties as in the linear finite element case. Therefore, the fully discrete convergence results transfer to high-order evolving surface finite elements using the mentioned error estimates of the Ritz map. The implementation of the high-order method is also a nontrivial task. The matrix assembly of the time-dependent mass and stiffness matrices is based on the usual reference element technique. Similar to isoparametric finite elements, the approximating surface is parameterized over the reference element, and therefore all the computations are done there. It was pointed out by Grande & Reusken (2014) that the approach of Demlow (2009) requires explicit knowledge of the exact signed distance function to the surface$${\it{\Gamma}}$$. However, the signed distance function is used only in the analysis, but it is not required for the numerical computations away from the initial time level. This is possible since the high-order element is uniquely determined by its elements, and by using the usual reference element technique. It is used only for generating the initial surface approximation. Here, we consider only linear parabolic PDEs on evolving surfaces; however, we believe our techniques and results carry over to other cases, such as to the Cahn–Hilliard equation (Elliott & Ranner, 2015), to wave equations (Lubich & Mansour, 2015), to ALE methods (Elliott & Venkataraman, 2014; Kovács & Power Guerra, 2014), nonlinear problems (Kovács & Power Guerra, 2016) and to evolving versions of bulk–surface problems (Elliott & Ranner, 2013). For more details, see a remark later on. Furthermore, while, in this article, we consider only evolving surfaces with prescribed velocity, many of the high-order geometric estimates of this article are essential for the numerical analysis of parabolic problems where the surface velocity depends on the solution (cf. Kovács et al., 2016). The article is organized in the following way. In Section 2, we recall the basics of linear parabolic problems on evolving surfaces together with some notation based on Dziuk & Elliott (2007). Section 3 deals with the description of higher-order evolving surface finite elements based on Demlow (2009)Section 4 contains the time discretizations and states the main results of this article: semidiscrete and fully discrete convergence estimates. In Section 5, we turn to the estimates of geometric errors and geometric perturbation estimates. In Section 6, the errors in the generalized Ritz map and its material derivatives are estimated. Section 7 contains the proof of the main results. Section 8 briefly describes the implementation of the high-order evolving surface finite elements. In Section 9, we present some numerical experiments illustrating our theoretical results. 2. The problem Let us consider a smooth-evolving compact surface, given by a smooth signed distance function $$d$$, $${\it{\Gamma}}(t) = \{ x\in \mathbb{R}^{m+1} \mid d(x,t) = 0 \} \subset \mathbb{R}^{m+1}$$ ($$m\leq 3$$), $$0 \leq t \leq T$$, which moves with a given smooth velocity $$v$$. Let $$\partial^{\bullet} u = \partial_{t} u + v \cdot \nabla u$$ denote the material derivative of the function $$u$$, where $$\nabla_{{\it{\Gamma}}}$$ is the tangential gradient given by $$\nabla_{{\it{\Gamma}}} u = \nabla u -\nabla u \cdot \mathrm{n} \mathrm{n}$$, with unit normal $$\mathrm{n}$$. We are sharing the setting of Dziuk & Elliott (2007) and Dziuk & Elliott (2013b). We consider the following linear problem on the above surface, for $$u=u(x,t)$$:   ∂∙u+u∇Γ(t)⋅v−ΔΓ(t)u =f on Γ(t),u(⋅,0) =u0 on Γ(0), (2.1) where for simplicity we set $$f=0$$, but all our results hold with a nonvanishing inhomogeneity as well. Let us briefly recall some important concepts used later on, whereas, in general, for basic formulae we refer to Dziuk & Elliott (2013a). An important tool is the Green’s formula on closed surfaces,   ∫Γ(t)∇Γ(t)z⋅∇Γ(t)ϕ=−∫Γ(t)(ΔΓ(t)z)ϕ. We use Sobolev spaces on surfaces: for a smooth surface $${\it{\Gamma}}(t)$$, for fixed $$t\in[0,T]$$ or for the space–time manifold given by $$\mathscr{G}_T= \cup_{t\in[0,T]} {\it{\Gamma}}(t)\times\{t\}$$, we define   H1(Γ(t))= {η∈L2(Γ(t))∣∇Γ(t)η∈L2(Γ(t))m+1},H1(GT)= {η∈L2(GT)∣∇Γ(t)η∈L2(Γ(t))m+1,∂∙η∈L2(Γ(t))} and analogously for higher-order versions $$H^{k}({\it{\Gamma}}(t))$$ and $$H^{k}(\mathscr{G}_T)$$ for $$k\in \mathbb{N}$$ (cf. Dziuk & Elliott, 2007, Section 2.1). The variational formulation of this problem reads as follows: find $$u\in H^1(\mathscr{G}_T)$$ such that   ddt∫Γ(t)uφ+∫Γ(t)∇Γ(t)u⋅∇Γ(t)φ=∫Γ(t)u∂∙φ (2.2) holds for almost every $$t\in(0,T)$$ for every $$\varphi(\cdot,t) \in H^1({\it{\Gamma}}(t))$$ with $$\partial^{\bullet}\varphi(\cdot,t) \in L^2({\it{\Gamma}}(t))$$ and $$u(\cdot,0)=u_0$$ holds. For suitable $$u_{0}$$, existence and uniqueness results for (2.2) were obtained in Dziuk & Elliott (2007, Theorem 4.4). 3. High-order evolving surface finite elements We define the high-order evolving surface finite element method (ESFEM) applied to our problem following Dziuk (1988), Dziuk & Elliott (2007) and Demlow (2009) We use simplicial elements and continuous piecewise polynomial basis functions of degree $$k$$. 3.1 Basic notions The smooth initial surface $${\it{\Gamma}}(0)$$ is approximated by a $$k$$-order interpolating discrete surface constructed in (a) below. This high-order approximation surface is then evolved in time by the a priori known surface velocity $$v$$, detailed in (b). The construction presented here is from Demlow (2009, Section 2). (a) Let $${\it{\Gamma}}_h^1(0)$$ be a polyhedron having triangular elements (denoted by $$E$$) with vertices lying on the initial surface $${\it{\Gamma}}(0)$$, forming a quasiuniform triangulation $$\mathscr{T}_h^1(0)$$, with mesh parameter $$h$$, and unit outward normal $$\mathrm{n}_h^1$$. Next, for $$k \geq 2$$, we define $${\it{\Gamma}}_h^k(0)$$, the $$k$$-order polynomial approximations to $${\it{\Gamma}}(0)$$. For a given fixed element $$E \in \mathscr{T}_h^1(0)$$, with nodes (numbered locally) $$x_1, x_2, \dotsc, x_{n_k}$$, the corresponding Lagrange basis functions of degree $$k$$ on $$E$$ are denoted by $$\tilde\chi^k_1, \tilde\chi^k_2, \dotsc, \tilde\chi^k_{n_k}$$. For arbitrary $$x \in E$$, we define the discrete projection   pk(x,0)=∑j=1nkp(xj,0)χ~jk(x,0), with Γ(0)∋p(xj,0)=xj−d(xj,0)n(p(xj,0)). This definition yields a continuous piecewise polynomial map on $${\it{\Gamma}}_h^1(0)$$. Then the $$k$$-order approximation surface is defined by   Γhk(0)={pk(x,0)∣x∈Γh1(0)}. The vertices of $${\it{\Gamma}}_h^k(0)$$ are denoted by $$a_i(0)$$, $$i=1,2,\dotsc,N$$ (here numbered globally). Note that these vertices are sitting on the exact surface. The high-order triangulation is denoted by $$\mathscr{T}_h^k(0)$$. Further small details can be found in Demlow (2009). (b) The surface approximation $${\it{\Gamma}}_h^k(t)$$ at time $$t\in[0,T]$$ is given by evolving the nodes of the initial triangulation by the velocity $$v$$ along the space-time manifold. For instance, the nodes $$a_i(t)$$ are determined by the ordinary differential equation (ODE) system   ddtai(t)=v(ai(t),t). Then the base triangulation and the basis functions move along as well. Hence, for $$0 \leq t \leq T$$, the nodes $$(a_i(t))_{i=1}^N$$ define an approximation of $${\it{\Gamma}}(t)$$ of degree $$k$$,   Γhk(t)={pk(x,t)∣x∈Γh1(t)}, defined analogously as for $$t=0$$. Therefore, the discrete surface $${\it{\Gamma}}_h^k(t)$$ remains an interpolation of $${\it{\Gamma}}(t)$$ for all times. The high-order triangulation at $$t$$ is denoted by $$\mathscr{T}_h^k(t)$$. Remark 3.1 It is highly important to note here that computing $$p(\cdot,t)$$ or $$p^k(\cdot,t)$$, for $$t>0$$, is never used during the numerical computations, hence explicit knowledge of the distance function is avoided (except for $$t=0$$). Instead, the positions of the nodes appearing in an element of $${\it{\Gamma}}_h^k(t)$$ are used to construct a reference mapping over the reference triangle $$E_0$$. In fact, the nodes (locally numbered here) $$a_1(t), a_2(t), \dotsc, a_{n_k}(t)$$ of the high-order element uniquely determine an approximating surface of degree $$k$$ (being an image of a polynomial of $$m$$ variables over $$E_0$$). The ESFEM matrices are then assembled using the reference element and the reference mapping. The discrete tangential gradient on the discrete surface $${\it{\Gamma}}_h^k(t)$$ is given by   ∇Γhk(t)ϕ:=∇ϕ−∇ϕ⋅nhknhk, understood in an elementwise sense, with $$\mathrm{n}_h^k$$ denoting the normal to $${\it{\Gamma}}_h^k(t)$$. For every $$t\in[0,T]$$ and for the $$k$$-order approximation surface $${\it{\Gamma}}_h^k(t)$$, we define the finite element subspace of order $$k$$, denoted by $$S_h^k(t)$$, spanned by the continuous evolving basis functions $$\chi_j$$ of piecewise degree $$k$$, satisfying $$\chi_j(a_i(t),t) = \delta_{ij}$$ for all $$i,j = 1, 2, \dotsc, N$$:   Shk(t)= {ϕ∈C(Γhk(t))∣ϕ=ϕ~∘pk(⋅,t)−1 where ϕ~∈S~hk(t)}= span{χ1(⋅,t),χ2(⋅,t),…,χN(⋅,t)}, where   χj(⋅,t)=χ~j(pk(⋅,t)−1,t)(j=1,2,…,N), with $$\tilde\chi_j( \cdot ,t)$$ being the degree $$k$$ polynomial basis function over the base triangulation $${\it{\Gamma}}_h^1$$, spanning the space $$\tilde S_h^k(t) = \{ \tilde\phi \in C({\it{\Gamma}}_h^1(t)) \mid \tilde\phi \text{ is a piecewise polynomial of degree } k \}$$. By construction $$S_h^k(t)$$ is an isoparametric finite element space. We interpolate the surface velocity on the discrete surface using the basis functions and denote it by $$V_h$$. Then the discrete material derivative is given by   ∂h∙ϕh=∂tϕh+Vh⋅∇ϕh(ϕh∈Shk(t)). The key transport property derived by Dziuk & Elliott (2007, Proposition 5.4) is   ∂h∙χj=0forj=1,2,…,N. (3.1) This result translates from the linear ESFEM case by the original proof of Dziuk & Elliott (2007, Section 5.2) using barycentric coordinates and the chain rule. The spatially discrete problem (for a fixed degree $$k$$) then reads as follows: find $$U_h\in S_h^k(t)$$, with $$\partial^{\bullet}_h U_h\in S_h^k(t)$$ and temporally continuous, such that   ddt∫Γhk(t)Uhϕh+∫Γhk(t)∇ΓhUh⋅∇Γhϕh=∫Γhk(t)Uh∂h∙ϕh(∀ϕh∈Shk(t) with ∂h∙ϕh∈Shk(t)), (3.2) with the initial condition $$U_h^0\in S_h^k(0)$$ being a suitable approximation to $$u_0$$. Later on we will always work with a high-order approximation surface $${\it{\Gamma}}_h^k(t)$$ and with the corresponding high-order evolving surface finite element space $$S_h^k(t)$$ (with $$2\leq k\in\mathbb{N}$$); therefore from now on, we drop the upper index $$k$$, unless we would like to emphasize it or it is not clear from the context. 3.2 The ODE system Similarly to Dziuk et al. (2012), the ODE form of the above problem (3.2) can be derived by setting   Uh(⋅,t)=∑j=1Nαj(t)χj(⋅,t) in the semidiscrete problem, and by testing with $$\phi_h=\chi_j$$ ($$j=1,2,\dotsc,N$$), and using the transport property (3.1). The spatially semidiscrete problem (3.2) is equivalent to the following ODE system for the vector $$\alpha(t)=(\alpha_j(t))_{j=1}^N \in\mathbb{R}^N$$, collecting the nodal values of $$U_h(\cdot,t)$$:   ddt(M(t)α(t))+A(t)α(t) =0,α(0) =α0, (3.3) where the evolving mass matrix $$M(t)$$ and the stiffness matrix $$A(t)$$ are defined as   M(t)|kj=∫Γh(t)χjχkandA(t)|kj=∫Γh(t)∇Γhχj⋅∇Γhχk, for $$j,k = 1,2, \dotsc, N$$. 3.3 Lifts For the error analysis we need to compare functions on different surfaces. This is conveniently done by the lift operator, which was introduced in Dziuk (1988) and further investigated in Dziuk & Elliott (2007, 2013b). The lift operator maps a function on the discrete surface onto a function on the exact surface. Let $${\it{\Gamma}}_h(t)$$ be a $$k$$-order approximation to the exact surface $${\it{\Gamma}}(t)$$. Using the oriented distance function $$d$$ (cf. Dziuk & Elliott, 2007, Section 2.1), the lift of a continuous function $$\eta_h \colon {\it{\Gamma}}_h(t) \rightarrow \mathbb{R}$$ is defined as   ηhl(p,t):=ηh(x,t),x∈Γ(t), where for every $$x\in {\it{\Gamma}}_h(t)$$ the point $$p=p(x,t)\in{\it{\Gamma}}(t)$$ is uniquely defined via   p=x−n(p,t)d(x,t). (3.4) By $$\eta^{-l}$$ we denote the function on $${\it{\Gamma}}_h(t)$$ whose lift is $$\eta$$. In particular, we will often use the space of lifted basis functions   (Sh(t))l=(Shk(t))l={ϕhl | ϕh∈Shk(t)}. 4. Convergence estimates 4.1 Convergence of the semidiscretization We now formulate the convergence theorem for the semidiscretization using high-order evolving surface finite elements. This result is the higher-order extension of Dziuk & Elliott (2013b, Theorem 4.4). Theorem 4.1 Consider the ESFEM of order $$k$$ as a space discretization of the parabolic problem (2.1). Let $$u$$ be a sufficiently smooth solution of the problem, and assume that the initial value satisfies   ‖uh0−u(⋅,0)‖L2(Γ(0))≤C0hk+1. Then there exists $$h_0>0$$, such that for mesh size $$h \leq h_0$$, the following error estimate holds, for $$t \leq T$$:   ‖uh(⋅,t)−u(⋅,t)‖L2(Γ(t))+h(∫0t‖∇Γ(s)(uh(⋅,s)−u(⋅,s))‖L2(Γ(s))2ds)1/2≤Chk+1. The constant $$C$$ is independent of $$h$$ and $$t$$ but depends on $$T$$. The proof of this result is postponed to a later section, after we have shown some preparatory results. 4.2 Time discretization: BDF We apply a $$p$$-step BDF for $$p\leq5$$ as a discretization of the ODE system (3.3), coming from the ESFEM space discretization of the parabolic evolving surface PDE. We briefly recall the $$p$$-step BDF method applied to system (3.3) with step size $$\tau>0$$:   1τ∑j=0pδjM(tn−j)αn−j+A(tn)αn=0(n≥p), (4.1) where the coefficients of the method are given by $$\delta(\zeta)=\sum_{j=0}^p \delta_j \zeta^j=\sum_{\ell=1}^p \frac{1}{\ell}(1-\zeta)^\ell$$, while the starting values are $$\alpha_0, \alpha_1, \dotsc, \alpha_{p-1}$$. The method is known to be $$0$$-stable for $$p\leq6$$ and have order $$p$$ (for more details, see Hairer & Wanner, 1996, Chapter V). In the following result, we compare the fully discrete solution   Uhn=∑j=1Nαjnχj(⋅,tn), obtained by solving (4.1) and the Ritz map $$\widetilde{\mathscr{P}}_h : H^{1}({\it{\Gamma}}(t)) \rightarrow S_h^k(t)$$ ($$t\in[0,T]$$) of the sufficiently smooth solution $$u$$. The precise definition of the Ritz map is given later. The $$H_h^{-1}({\it{\Gamma}}_h(t))$$-norm of the ESFEM function $$R_h$$ is defined as   ‖Rh(⋅,t)‖Hh−1(Γh(t))=sup0≠ϕh∈Shk(t)mh(Rh(⋅,t),ϕh)‖ϕh‖H1(Γh(t)) . The following error bound was shown in Lubich et al. (2013) for BDF methods up to order 5 (see also Mansour, 2013). Theorem 4.2 (Lubich et al., 2013, Theorem 5.1) Consider the parabolic problem (2.1), having a sufficiently smooth solution for $$0\leq t \leq T$$. Couple the $$k$$-order ESFEM as space discretization with time discretization by a $$p$$-step backward difference formula with $$p\leq5$$. Assume that the Ritz map of the solution has continuous discrete material derivatives up to order $$p+1$$. Then there exists $$\tau_0>0$$, independent of $$h$$, such that for $$\tau \leq \tau_0$$, for the error $$E_h^n=U_h^n-\widetilde{\mathscr{P}}_h u(\cdot,t_n)$$ the following estimate holds for $$t_n=n\tau \leq T$$:   ‖Ehn‖L2(Γh(tn))+ (τ∑j=1n‖∇Γh(tj)Ehj‖L2(Γh(tj))2)1/2≤Cβ~h,pτp+ (τ∑j=1n‖Rh(⋅,tj)‖Hh−1(Γh(tj))2)1/2+Cmax0≤i≤p−1‖Ehi‖L2(Γh(ti)), where $$R_h$$ is the high-order ESFEM residual. The constant $$C$$ is independent of $$h, n$$ and $$\tau$$ but depends on $$T$$. Furthermore,   β~h,p2=∫0T∑ℓ=1p+1‖(∂h∙)(ℓ)(P~hu)(⋅,t)‖L2(Γh(t))dt. Remark 4.3 The most important technical tools in the proof of Theorem 4.2, and also in the proof of the BDF stability result in Lubich et al. (2013, Lemma 4.1), are the ODE formulation (3.3) and the key estimates first shown in Dziuk et al. (2012, Lemma 4.1), where the following estimates are shown: there exist $$\mu,\kappa>0$$ such that, for $$w,z \in \mathbb{R}^N$$,   wT(M(s)−M(t))z≤(eμ(s−t)−1)‖w‖M(t)‖z‖M(t),wT(A(s)−A(t))z≤(eκ(s−t)−1)‖w‖A(t)‖z‖A(t). The proof of these bounds involves only basic properties of the mass and stiffness matrices $$M(t)$$ and $$A(t)$$; hence it is independent of the order of the basis functions. Hence, the high-order ESFEM versions of these two inequalities also hold. Therefore, the original results of Lubich et al. (2013) hold here as well. Later on, when suitable results are at hand, we give some remarks on the smoothness assumption of the Ritz map. 4.3 Convergence of the full discretization We are now in a position to formulate one of the main results of this article, which yields optimal-order error bounds for high-order finite element semidiscretization coupled to BDF methods up to order 5 applied to an evolving surface PDE. Theorem 4.4 ($$k$$-order ESFEM and BDF-$$p$$) Consider the ESFEM of order $$k$$ as space discretization of the parabolic problem (2.1), coupled to the time discretization by a $$p$$-step backward difference formula with $$p\leq5$$. Let $$u$$ be a sufficiently smooth solution of the problem and assume that the starting values satisfy (with $$\mathscr{P}_h u=(\widetilde{\mathscr{P}}_h u)^l)$$)   max0≤i≤p−1‖uhi−(Phu)(⋅,ti)‖L2(Γ(ti))≤C0hk+1. Then there exist $$h_0>0$$ and $$\tau_0>0$$, such that for $$h\leq h_0$$ and $$\tau \leq \tau_0$$, the following error estimate holds for $$t_n=n\tau \leq T$$:   ‖uhn−u(⋅,tn)‖L2(Γ(tn))+h(τ∑j=pn‖∇Γ(tj)(uhj−u(⋅,tj))‖L2(Γ(tj))2)1/2≤C(τp+hk+1). The constant $$C$$ is independent of $$h, \ \tau$$ and $$n$$ but depends on $$T$$. The proof of this result is also postponed to a later section, after we have shown some preparatory results. Remark 4.5 We remark here that an analogous fully discrete convergence result is readily available for algebraically stable implicit Runge–Kutta methods (such as the Radau IIA methods), since the Runge–Kutta analogue of Theorem 4.2 has been proved in Dziuk et al. (2012). The combination of this result with our high-order semidiscrete error bounds (Theorem 7.1) proves the Runge–Kutta analogue of Theorem 4.4. 5. Geometric estimates In this section, we present further notation and some technical lemmas that will be used later on in the proofs leading to the convergence result. These estimates are high-order analogues of some previous results proved in Dziuk (1988), Dziuk & Elliott (2007), Demlow (2009), Mansour (2013) and Dziuk & Elliott (2013b). 5.1 Geometric approximation results In the following, we state and prove estimates for the errors resulting from the geometric surface approximation. Most of these estimates hold for a sufficiently small $$h$$, which here means for $$h\leq h_0$$ with a sufficiently small $$h_0>0$$. Lemma 5.1 (Equivalence of norms, Dziuk, 1988; Demlow, 2009) Let $$\eta_h : {\it{\Gamma}}_h(t) \rightarrow \mathbb{R}$$ with lift $$\eta_h^l : {\it{\Gamma}}(t) \rightarrow \mathbb{R}$$. Then, for a sufficiently small $$h$$, the discrete and continuous $$L^p$$ and Sobolev norms are equivalent, independently of the mesh size $$h$$. For instance, there is a constant $$c>0$$ such that for all $$h$$ sufficiently small,   c−1‖ηh‖L2(Γh(t))≤ ‖ηhl‖L2(Γ(t))≤c‖ηh‖L2(Γh(t)),c−1‖ηh‖H1(Γh(t))≤ ‖ηhl‖H1(Γ(t))≤c‖ηh‖H1(Γh(t)). We now turn to the study of some geometric concepts and their errors. By $$\delta_h$$ we denote the quotient between the continuous and discrete surface measures, $$\textrm{d} A$$ and $$\textrm{d} A_h$$, defined as $$\delta_h \textrm{d} A_h = \textrm{d} A$$. Further, we recall that $${\text{Pr}}$$ and $${\text{Pr}}_h$$ are the projections onto the tangent spaces of $${\it{\Gamma}}(t)$$ and $${\it{\Gamma}}_h(t)$$, respectively. We further set, from Dziuk & Elliott (2013b),   Qh=1δh(Id−dH)PrPrhPr(Id−dH), (5.1) where $$\mathscr{H}$$ ($$\mathscr{H}_{ij} = \partial_{x_j}\mathrm{n}_i$$) is the (extended) Weingarten map. Using this notation and (3.4), in the proof of Dziuk & Elliott (2013b, Lemma 5.5), it is shown that   ∇Γhϕh(x)⋅∇Γhϕh(x)=δhQh∇Γϕhl(p)⋅∇Γϕhl(p). (5.2) For these quantities we show some results analogous to their linear ESFEM version showed in Dziuk & Elliott (2013b, Lemma 5.4), Demlow (2009, Proposition 2.3) or Mansour (2013, Lemma 6.1). Later on the following estimates will play a key technical role. Lemma 5.2 For $${\it{\Gamma}}_h^k(t)$$ and $${\it{\Gamma}}(t)$$ as above, for a sufficiently small $$h$$, we have the following geometric approximation estimates:    ‖d‖L∞(Γh(t))≤chk+1, ‖1−δh‖L∞(Γh(t))≤chk+1, ‖n−nh‖L∞(Γh(t))≤chk, and ‖Id−δhQh‖L∞(Γh(t))≤chk+1, ‖(∂h∙)(ℓ)d‖L∞(Γh(t))≤chk+1, ‖(∂h∙)(ℓ)δh‖L∞(Γh(t))≤chk+1, ‖Pr((∂h∙)(ℓ)Qh)Pr‖L∞(Γh)≤chk+1, with constants depending only on $$\mathscr{G}_T$$ but not on $$h$$ or $$t$$. The first three bounds were shown in Demlow (2009, Proposition 2.3) for the stationary case. Noting that the constants depend only on $$\mathscr{G}_T$$ these inequalities are shown. The last four bounds are simply the higher-order extensions of the corresponding estimates of Mansour (2013, Lemma 6.1). They can be proved in the exact same way using the bounds of the first three estimates. These proofs are included in the Appendix. 5.2 Interpolation estimates for evolving surface finite elements The following result gives estimates for the error in the interpolation. Our setting follows that of Demlow (2009, Section 2.5). Let us assume that the surface $${\it{\Gamma}}(t)$$ is approximated by the interpolation surface $${\it{\Gamma}}_h^k(t)$$. Then for any $$w\in H^{k+1}({\it{\Gamma}}(t))$$, there is a unique $$k$$-order surface finite element interpolation $$\widetilde I_h^k w \in S_h^k(t)$$; furthermore we set $$(\widetilde I_h^k w)^l = I_h^k w$$. Lemma 5.3 (Demlow, 2009, Proposition 2.7) Let $$w:\mathscr{G}_T \rightarrow \mathbb{R}$$ such that $$w\in H^{k+1}({\it{\Gamma}}(t))$$ for $$0 \leq t \leq T$$. There exists a constant $$c>0$$ depending on $$\mathscr{G}_T$$, but independent of $$h$$ and $$t$$, such that for $$0 \leq t \leq T$$ and for a sufficiently small $$h$$,   ‖w−Ihkw‖L2(Γ(t))+h‖∇Γ(w−Ihkw)‖L2(Γ(t))≤ chk+1‖w‖Hk+1(Γ(t)). We distinguish the special case of a linear surface finite element interpolation on $${\it{\Gamma}}_h^k(t)$$. For $$w \in H^2({\it{\Gamma}}(t))$$, the linear surface finite element interpolant is denoted by $$I_h^{(1)} w$$, and it satisfies, with $$c>0$$,   ‖w−Ih(1)w‖L2(Γ(t))+h‖∇Γ(w−Ih(1)w)‖L2(Γ(t))≤ ch2‖w‖H2(Γ(t)). Note that for $$I_h^{(1)}$$ the underlying approximating surface $${\it{\Gamma}}_h^k(t)$$ is still of high order. For $$k=1$$, $$I_h^1$$ and $$I_h^{(1)}$$ simply coincide. The upper $$k$$ indices are again dropped later on. 5.3 Velocity of lifted material points and material derivatives Following Dziuk & Elliott (2013b) and Lubich & Mansour (2015) we define the velocity of the lifted material points, denoted by $$v_h$$, and the corresponding discrete material derivative $$\partial^{\bullet}_h$$. For arbitrary $$y(t) = p(x(t),t) \in {\it{\Gamma}}(t)$$, with $$x(t) \in {\it{\Gamma}}_h(t)$$, cf. (3.4), we have   ddty(t)=vh(y(t),t)=∂tp(x(t),t)+Vh(x(t),t)⋅∇p(x(t),t); (5.3) hence for $$y=p(x,t)$$ (see Dziuk & Elliott, 2013b),   vh(y,t)=(Pr−dH)(x,t)Vh(x,t)−∂td(x,t)n(x,t)−d(x,t)∂tn(x,t). Following Lubich & Mansour (2015, Section 7.3), we note that $$- \partial_t d(x,t)\mathrm{n}(x,t)$$ is the normal component of $$v(p,t)$$ and that the other terms are tangent to $${\it{\Gamma}}(t)$$; hence   v−vhis a tangent vector. (5.4) It is also important to note that $$v_h \neq V_h^l$$ (cf. Dziuk & Elliott, 2013b). The discrete material derivative of the lifted points on $${\it{\Gamma}}(t)$$ reads   ∂h∙φh=∂tφh+vh⋅φh(φh∈(Sh(t))l). The lifted basis functions also satisfy the transport property   ∂h∙χjl=0(j=1,2,…,N). To prove error estimates for higher-order material derivatives of the Ritz map, we need high-order bounds for the error between the continuous velocity $$v$$ and the discrete velocity $$v_h$$. We generalize here Dziuk & Elliott (2013b, Lemma 5.6) and Lubich & Mansour (2015, Lemma 7.3) to the high-order case. Lemma 5.4 For $$\ell \geq 0$$, there exists a constant $$c_\ell>0$$ depending on $$\mathscr{G}_T$$, but independent of $$t$$ and $$h$$, such that, for a sufficiently small $$h$$,   ‖(∂h∙)(ℓ)(v−vh)‖L∞(Γ(t))+h‖∇Γ(∂h∙)(ℓ)(v−vh)‖L∞(Γ(t))≤cℓhk+1. Proof. We follow the steps of the original proofs from Dziuk & Elliott (2013b) and Lubich & Mansour (2015). (a) For $$\ell=0$$. Using the definition (5.3) and the fact $$V_h = \widetilde I_h v$$, we have   |v(p,t)−vh(p,t)|=|Pr(v−Ihv)(p,t)+d(HIhv(p,t)+∂tn)|≤chk+1, where we have used the interpolation estimates, Lemma 5.2 and the boundedness of the other terms. For the gradient estimate we use the fact that $$\nabla_{{\it{\Gamma}}} d = \nabla_{{\it{\Gamma}}} \partial^{\bullet}_h d = 0$$ and the geometric bounds of Lemma 5.2:   |∇Γ(v−vh)| ≤c|v−Ihv|+c|∇Γ(v−Ihv)| +|(∇Γd)(HIhv+∂tn)|+|d(∇Γ(HIhv+∂tn))| ≤chk. (b) For $$\ell=1$$, we use the transport property and again Lemma 5.2:   |∂h∙(v−vh)| ≤|(∂h∙Pr)(v−Ihv)|+|Pr(∂h∙v−Ih∂h∙v)| +|(∂h∙d)(HIhv+∂tn)|+|d(∂h∙(HIhv+∂tn))| ≤chk+1. Again using $$\nabla_{{\it{\Gamma}}} d = \nabla_{{\it{\Gamma}}} \partial^{\bullet}_h d = 0$$ and the geometric bounds of Lemma 5.2,   |∇Γ∂h∙(v−vh)| ≤c|v−Ihv|+c|∇Γ(v−Ihv)| +c|∂h∙(v−Ihv)|+c|∇Γ(∂h∙v−Ih∂h∙v)| +|(∇Γ∂h∙d)(HIhv+∂tn)|+|d(∇Γ∂h∙(HIhv+∂tn))| ≤chk. (c) For $$\ell >1$$, the proof uses similar arguments. □ 5.4 Bilinear forms and their estimates We use the time-dependent bilinear forms defined in Dziuk & Elliott (2013b): for arbitrary $$z,\varphi \in H^1({\it{\Gamma}}(t))$$ and for their discrete analogues for $$Z_h, \phi_h \in S_h(t)$$,   m(t;z,φ) =∫Γ(t)zφ,a(t;z,φ) =∫Γ(t)∇Γz⋅∇Γφ,g(t;v;z,φ) =∫Γ(t)(∇Γ⋅v)zφ,b(t;v;z,φ) =∫Γ(t)B(v)∇Γz⋅∇Γφ,mh(t;Zh,ϕh) =∫Γh(t)Zhϕh,ah(t;Zh,ϕh) =∫Γh(t)∇ΓhZh⋅∇Γhϕh,gh(t;Vh;Zh,ϕh) =∫Γh(t)(∇Γh⋅Vh)Zhϕh,bh(t;Vh;Zh,ϕh) =∫Γh(t)Bh(Vh)∇ΓhZh⋅∇Γhϕh, where the discrete tangential gradients are understood in a piecewise sense, and with the tensors given by   B(v)|ij =δij(∇Γ⋅v)−((∇Γ)ivj+(∇Γ)jvi),Bh(Vh)|ij =δij(∇Γh⋅Vh)−((∇Γh)i(Vh)j+(∇Γh)j(Vh)i), for $$i,j=1,2,\dotsc,m+1$$. For more details, see Dziuk & Elliott (2013b, Lemma 2.1) (and the references in the proof) or Dziuk & Elliott (2013a, Lemma 5.2). We will omit the time dependency of the bilinear forms if it is clear from the context. 5.4.1 Discrete material derivative and transport properties The following transport equations, especially the one with discrete material derivatives on the continuous surface, are of great importance during the defect estimates later on. Lemma 5.5 (Dziuk & Elliott, 2013b, Lemma 4.2) Consider $${\it{\Gamma}}(t)$$ as the lift of the discrete surface $${\it{\Gamma}}_h^k(t)$$ (i.e., $${\it{\Gamma}}(t)$$ can be decomposed into curved elements that are lifts of the elements of $${\it{\Gamma}}_h^k(t)$$), moving with the velocity $$v_h$$ from (5.3). Then for any $$z,\varphi, \partial^{\bullet}_h z, \partial^{\bullet}_h \varphi \in H^1({\it{\Gamma}}(t))$$ we have   ddtm(z,φ) =m(∂h∙z,φ)+m(z,∂h∙φ)+g(vh;z,φ),ddta(z,φ) =a(∂h∙z,φ)+a(z,∂h∙φ)+b(vh;z,φ). The same formulae hold for $${\it{\Gamma}}(t)$$ when $$\partial^{\bullet}_h$$ and $$v_h$$ are replaced with $$\partial^{\bullet}$$ and $$v$$, respectively. Similarly, in the discrete case, for arbitrary $$z_h, \phi_h, \partial^{\bullet}_hz_h, \partial^{\bullet}_h\phi_h \in S_h(t)$$, we have   ddtmh(zh,ϕh) =mh(∂h∙zh,ϕh)+mh(zh,∂h∙ϕh)+gh(Vh;zh,ϕh),ddtah(zh,ϕh) =ah(∂h∙zh,ϕh)+ah(zh,∂h∙ϕh)+bh(Vh;zh,ϕh). 5.4.2 Geometric perturbation errors The following estimates are the most important technical results of this paper. Later on, they will play a crucial role in the defect estimates. We note here that these results extend the first-order ESFEM theory of Dziuk & Elliott (2013b) (for the first three inequalities) and Lubich & Mansour (2015) (for the last inequality) to the higher-order ESFEM case. The high-order version of the inequalities for time-independent bilinear forms $$a(\cdot,\cdot)$$ and $$m(\cdot,\cdot)$$ were shown in Elliott & Ranner (2013, Lemma 6.2). The following results generalizes these for the time-dependent case. Lemma 5.6 For any $$Z_h,\phi_h \in S_h^k(t)$$, and for their lifts $$Z_h^l,\phi_h^l \in H^1({\it{\Gamma}}(t))$$, we have the following bounds, with a sufficiently small $$h$$:   |m(Zhl,ϕhl)−mh(Zh,ϕh)| ≤chk+1‖Zhl‖L2(Γ(t))‖ϕhl‖L2(Γ(t)),|a(Zhl,ϕhl)−ah(Zh,ϕh)| ≤chk+1‖∇ΓZhl‖L2(Γ(t))‖∇Γϕhl‖L2(Γ(t)),|g(vh;Zhl,ϕhl)−gh(Vh;Zh,ϕh)| ≤chk+1‖Zhl‖L2(Γ(t))‖ϕhl‖L2(Γ(t)),|b(vh;Zhl,ϕhl)−bh(Vh;Zh,ϕh)| ≤chk+1‖∇ΓZhl‖L2(Γ(t))‖∇Γϕhl‖L2(Γ(t)), where the constants $$c>0$$ are independent of $$h$$ and $$t$$, but depend on $$\mathscr{G}_T$$. Proof. The proof of the first two estimates is a high-order generalization of Dziuk & Elliott (2013b, Lemma 5.5), while the proof of the last two estimates is a high-order extension of Lubich & Mansour (2015, Lemma 7.5). Their proofs follow these references. Both parts use the geometric estimates shown in Lemma 5.2. To show the first inequality, we estimate   |m(Zhl,ϕhl)−mh(Zh,ϕh)| =|∫Γ(t)ZhlϕhldA−∫Γh(t)ZhϕhdAh| =|∫Γ(t)(1−δh−1)ZhlϕhldA| ≤chk+1‖Zhl‖L2(Γ(t))‖ϕhl‖L2(Γ(t)). For the second inequality, similarly we have   |a(Zhl,ϕhl)−ah(Zh,ϕh)| =|∫Γ(t)∇ΓZhl⋅∇ΓϕhldA−∫Γh(t)∇ΓhZh⋅∇ΓhϕhdAh| =|∫Γ(t)(Id−δhQh)∇ΓZhl⋅∇ΓϕhldA| ≤chk+1‖∇ΓZhl‖L2(Γ(t))‖∇Γϕhl‖L2(Γ(t)). For the third estimate we start by taking the time derivative of the equality $$m(Z_h^l,\phi_h^l) = m_h(Z_h,\phi_h \delta_h)$$, using the first transport property from Lemma 5.5 to obtain   ddtm(Zhl,ϕhl) =m(∂h∙Zhl,ϕhl)+m(Zhl,∂h∙ϕhl)+g(vh;Zhl,ϕhl) =ddtmh(Zh,ϕhδh) =mh(∂h∙Zh,ϕhδh)+mh(Zh,(∂h∙ϕh)δh) +gh(Vh;Zh,ϕhδh)+mh(Zh,(∂h∙δh)ϕh). Hence, using $$\partial^{\bullet}_h w_h^l = (\partial^{\bullet}_h w_h)^l$$,   g(vh;Zhl,ϕhl)−gh(Vh;Zh,ϕhδh) = m((∂h∙Zh)l,ϕhl)−mh(∂h∙Zh,ϕhδh) +m(Zhl,(∂h∙ϕh)l)−mh(Zh,(∂h∙ϕh)δh)+mh(Zh,(∂h∙δh)ϕh) =mh(Zh,(∂h∙δh)ϕh). Hence, together with Lemma 5.2, the bound   |g(vh;Zhl,ϕhl)−gh(Vh;Zh,ϕh)| ≤|gh(Vh;Zh,ϕh(δh−1))|+|mh(Zh,(∂h∙δh)ϕh)| ≤c(‖∂h∙δh‖L∞(Γh(t))+‖1−δh‖L∞(Γh(t)))‖Zhl‖L2(Γ(t))‖ϕhl‖L2(Γ(t)) finishes the proof of the third estimate. For the last inequality we take a similar approach, using the second transport property from Lemma 5.5 and the relation (5.2) to obtain   ddtah(Zh,ϕh) =ah(∂h∙Zh,ϕh)+ah(Zh,∂h∙ϕh)+bh(Vh;Zh,ϕh) =ddt∫Γ(t)Qhl∇ΓZhl⋅∇Γϕhl =∫Γ(t)Qhl∇Γ∂h∙Zhl⋅∇Γϕhl+∫Γ(t)Qhl∇ΓZhl⋅∇Γ∂h∙ϕhl +∫Γ(t)(∂h∙Qhl)∇ΓZhl⋅∇Γϕhl+∫Γ(t)B(vh)Qhl∇ΓZhl⋅∇Γϕhl. Again, using $$\partial^{\bullet}_h w_h^l = (\partial^{\bullet}_h w_h)^l$$, together with (5.2) and the geometric estimates of Lemma 5.2 we obtain   |bh(Vh;Zh,ϕh)−b(vh;Zhl,ϕhl)| =|∫Γ(t)(∂h∙Qhl)∇ΓZhl⋅∇Γϕhl+∫Γ(t)B(vh)(Qhl−Id)∇ΓZhl⋅∇Γϕhl| ≤chk+1‖∇ΓZhl‖L2(Γ(t))‖∇Γϕhl‖L2(Γ(t)), completing the proof. □ 6. Generalized Ritz map and higher-order error bounds We recall the generalized Ritz map for evolving surface PDEs from Lubich & Mansour (2015). Definition 6.1 (Ritz map) For any given $$z\in H^{1}({\it{\Gamma}}(t))$$, there is a unique $$\widetilde{\mathscr{P}}_h z\in S_h^k(t)$$ such that for all $$\phi_h\in S_h^k(t)$$, with the corresponding lift $$\varphi_h=\phi_h^l$$, we have   ah∗(P~hz,ϕh)=a∗(z,φh), (6.1) where we let $$a^{\ast}=a+m$$ and $$a_h^{\ast}=a_h+m_h$$, so that the forms $$a$$ and $$a_h$$ are positive definite. Then $$\mathscr{P}_h z \in (S_h^k(t))^l$$ is defined as the lift of $$\widetilde{\mathscr{P}}_h z$$, i.e., $$\mathscr{P}_h z = (\widetilde{\mathscr{P}}_h z)^l$$. We note here that originally in Lubich & Mansour (2015, Definition 8.1) an extra term appeared involving $$\partial^{\bullet} z$$ and the surface velocity, which is not needed for the parabolic case. The Ritz map above is still well defined. Galerkin orthogonality does not hold in this case, just up to a small defect. Lemma 6.2 (Galerkin orthogonality up to a small defect) For any $$z\in H^{1}({\it{\Gamma}}(t))$$ and for all $$\varphi_{h}\in (S_h^k(t))^l$$ and for a sufficiently small $$h$$,   |a∗(z−Phz,φh)|≤chk+1‖Phz‖H1(Γ(t))‖φh‖H1(Γ(t)), (6.2) where $$c$$ is independent of $$\xi$$, $$h$$ and $$t$$. Proof. Using the definition of the Ritz map and Lemma 5.6, we estimate   |a∗(z−Phz,φh)| =|ah∗(P~hz,ϕh)−a∗(Phz,φh)|≤chk+1‖Phz‖H1(Γ(t))‖φh‖H1(Γ(t)). □ 6.1 Errors in the Ritz map Now we prove higher-order error estimates for the Ritz map (6.1) and also for its material derivatives; the analogous results for the first-order ESFEM case can be found in Dziuk & Elliott (2013b, Section 6) or in Mansour (2013, Section 7). 6.1.1 Error bounds for the Ritz map Theorem 6.3 Let $$z:\mathscr{G}_T \rightarrow \mathbb{R}$$ with $$z\in H^{k+1}({\it{\Gamma}}(t))$$ for every $$0 \leq t \leq T$$. Then the error in the Ritz map satisfies the bound, for $$0 \leq t \leq T$$ and for $$h \leq h_0$$ with a sufficiently small $$h_0$$,   ‖z−Phz‖L2(Γ(t))+h‖z−Phz‖H1(Γ(t))≤chk+1‖z‖Hk+1(Γ(t)), where the constant $$c>0$$ is independent of $$h$$ and $$t$$. Proof. To ease the presentation, we suppress all time arguments $$t$$ appearing in the norms within the proof (except for some special occasions). (a) We first prove the gradient estimate. Starting with the definition of the $$H^1({\it{\Gamma}}(t))$$-norm, then using the estimate (6.2), we have   ‖z−Phz‖H1(Γ)2 =a∗(z−Phz,z−Phz) =a∗(z−Phz,z−Ihz)+a∗(z−Phz,Ihz−Phz) ≤‖z−Phz‖H1(Γ)‖z−Ihz‖H1(Γ)+chk+1‖Phz‖H1(Γ)‖Ihz−Phz‖H1(Γ) ≤chk‖z−Phz‖H1(Γ)‖z‖Hk+1(Γ) +chk+1(2‖z−Phz‖H1(Γ)2+‖z‖H1(Γ)2+ch2k‖z‖Hk+1(Γ)2), using interpolation error estimate, and for the second term, we used the estimate   ‖Phz‖H1(Γ)‖Ihz−Phz‖H1(Γ) ≤(‖Phz−z‖H1(Γ)+‖z‖H1(Γ))(‖Ihz−z‖H1(Γ)+‖z−Phz‖H1(Γ)) ≤2‖z−Phz‖H1(Γ)2+‖z‖Hk+1(Γ)2+ch2k‖z‖Hk+1(Γ)2. Now using Young’s inequality, and for a sufficiently small (but $$t$$ independent) $$h\leq h_0$$, we have the gradient estimate   ‖z−Phz‖H1(Γ(t))≤chk‖z‖Hk+1(Γ(t)). (b) The $$L^2$$-estimate follows from the Aubin–Nitsche trick. Let us consider the problem   −ΔΓ(t)w+w=z−PhzonΓ(t); then the usual elliptic theory (see, e.g., Dziuk & Elliott, 2013a, Section 3.1; Aubin, 1998) yields the following: the solution $$w\in H^2({\it{\Gamma}}(t))$$ satisfies the bound, with $$c>0$$ independent of $$t$$,   ‖w‖H2(Γ(t))≤c‖z−Phz‖L2(Γ(t)). By testing the elliptic weak problem with $$z-\mathscr{P}_h z$$, using (6.2) again, and using the linear finite element interpolation $$I_h^{(1)}$$ on $${\it{\Gamma}}_h^k(t)$$, we obtain   ‖z−Phz‖L2(Γ)2 =a∗(z−Phz,w) =a∗(z−Phz,w−Ih(1)w)+a∗(z−Phz,Ih(1)w) ≤‖z−Phz‖H1(Γ)‖w−Ih(1)w‖H1(Γ)+chk+1‖Phz‖H1(Γ)‖Ih(1)w‖H1(Γ) ≤chk‖z‖Hk+1(Γ)ch‖w‖H2(Γ)+chk+1‖Phz‖H1(Γ)‖Ih(1)w‖H1(Γ), where for the second term we have now used   ‖Phz‖H1(Γ)‖Ih(1)w‖H1(Γ) ≤(‖z−Phz‖H1(Γ)+‖z‖H1(Γ))(‖w−Ih(1)w‖H1(Γ)+‖w‖H1(Γ)) ≤(1+chk)‖z‖Hk+1(Γ)(1+ch)‖w‖H2(Γ). Then the combination of the gradient estimate for the Ritz map and the interpolation error yields   ‖z−Phz‖L2(Γ(t))1c‖w‖H2(Γ(t))≤‖z−Phz‖L2(Γ(t))2≤chk+1‖z‖Hk+1(Γ(t))‖w‖H2(Γ(t)), which completes the proof. □ 6.1.2 Error bounds for the material derivatives of the Ritz map Since, in general, $$\partial^{\bullet}_h \mathscr{P}_h z = \mathscr{P}_h \partial^{\bullet}_h z$$ does not hold, we need the following result. Theorem 6.4 The error in the material derivatives of the Ritz map, for any $$\ell\in\mathbb{N}$$, satisfies the following bounds, for $$0 \leq t \leq T$$ and for $$h \leq h_0$$ with a sufficiently small $$h_0$$:   ‖(∂h∙)(ℓ)(z−Phz)‖L2(Γ(t))+h ‖∇Γ(∂h∙)(ℓ)(z−Phz)‖L2(Γ(t))≤cℓhk+1∑j=0ℓ‖(∂∙)(j)z‖Hk+1(Γ(t)), where the constant $$c_\ell>0$$ is independent of $$h$$ and $$t$$. Proof. The proof is a modification of Mansour (2013, Theorem 7.3). Again, to ease the presentation, we suppress the $$t$$ argument of the surfaces norms (except for some special occasions). For $$\ell=1$$: (a) We start by taking the time derivative of the definition of the Ritz map (6.1), use the transport properties Lemma 5.5 and apply the definition of the Ritz map once more; we arrive at   a∗(∂h∙z,φh) =−b(vh;z,φh)−g(vh;z,φh) +ah∗(∂h∙P~hz,ϕh)+bh(Vh;P~hz,ϕh)+gh(Vh;P~hz,ϕh). Then we obtain   a∗(∂h∙z−∂h∙Phz,φh) =−b(vh;z−Phz,φh)−g(vh;z−Phz,φh) +F1(φh), (6.3) where   F1(φh) =(ah∗(∂h∙P~hz,ϕh)−a∗(∂h∙Phz,φh)) +(bh(Vh;P~hz,ϕh)−b(vh;Phz,φh)) +(gh(Vh;P~hz,ϕh)−g(vh;Phz,φh)). Using the geometric estimates of Lemma 5.6, $$F_1$$ can be estimated as   |F1(φh)|≤chk+1(‖∂h∙Phz‖H1(Γ(t))+‖Phz‖H1(Γ(t)))‖φh‖H1(Γ(t)). (6.4) The velocity error estimate Lemma 5.4 yields   ‖∂h∙z‖H1(Γ)≤‖∂∙z‖H1(Γ)+chk‖z‖H2(Γ). Then using $$\partial^{\bullet}_h \mathscr{P}_h z$$ as a test function in (6.3), and using the error estimates of the Ritz map, together with the estimates above, with $$h\leq h_0$$, we have   ‖∂h∙Phz‖H1(Γ)2 = a∗(∂h∙Phz,∂h∙Phz) =b(vh;z−Phz,∂h∙Phz)+g(vh;z−Phz,∂h∙Phz)+a∗(∂h∙z,∂h∙Phz)−F1(∂h∙Phz) ≤ chk‖z‖Hk+1(Γ)‖∂h∙Phz‖H1(Γ)+‖∂h∙z‖H1(Γ)‖∂h∙Phz‖H1(Γ) +chk+1(‖∂h∙Phz‖H1(Γ)+‖Phz‖H1(Γ))‖∂h∙Phz‖H1(Γ) ≤chk‖z‖Hk+1(Γ)‖∂h∙Phz‖H1(Γ)+(‖∂∙z‖H1(Γ)+chk‖z‖H2(Γ))‖∂h∙Phz‖H1(Γ) +chk+1(‖∂h∙Phz‖H1(Γ)+‖z−Phz‖H1(Γ)+‖z‖H1(Γ))‖∂h∙Phz‖H1(Γ) ≤(chk‖z‖Hk+1(Γ)+‖∂∙z‖H1(Γ))‖∂h∙Phz‖H1(Γ)+chk+1‖∂h∙Phz‖H1(Γ)2, absorption using an $$h\leq h_0$$, with a sufficiently small $$h_0>0$$, and dividing through yields   ‖∂h∙Phz‖H1(Γ)≤c‖∂∙z‖H1(Γ)+chk‖z‖Hk+1(Γ). Combining all the previous estimates and using Young’s inequality, the Cauchy–Schwarz inequality and Theorem 6.3, for a sufficiently small $$h\leq h_0$$, we obtain   a∗(∂h∙z−∂h∙Phz,φh) ≤ c‖z−Phz‖H1(Γ)‖φh‖H1(Γ) +chk+1(‖∂h∙Phz‖H1(Γ)+‖Phz‖H1(Γ))‖φh‖H1(Γ) ≤ chk‖z‖Hk+1(Γ)‖φh‖H1(Γ) +chk+1(‖∂∙z‖H1(Γ)+(1+chk)‖z‖Hk+1(Γ))‖φh‖H1(Γ) ≤ chk(‖z‖Hk+1(Γ)+h‖∂∙z‖H1(Γ))‖φh‖H1(Γ). Then, similarly to the previous proof, we have   ‖∂h∙z−∂h∙Phz‖H1(Γ)2 ≤a∗(∂h∙z−∂h∙Phz,∂h∙z−∂h∙Phz) =a∗(∂h∙z−∂h∙Phz,∂h∙z−Ih∂∙z)+a∗(∂h∙z−∂h∙Phz,Ih∂∙z−∂h∙Phz) ≤‖∂h∙z−∂h∙Phz‖H1(Γ)‖∂h∙z−Ih∂∙z‖H1(Γ) +chk(‖z‖Hk+1(Γ)+h‖∂∙z‖H1(Γ))‖Ih∂∙z−∂h∙Phz‖H1(Γ). Finally, the interpolation estimates, Young’s inequality and absorption using a sufficiently small $$h\leq h_0$$, yields the gradient estimate. (b) The $$L^2$$ estimate again follows from the Aubin–Nitsche trick. Let us now consider the problem   −ΔΓ(t)w+w=∂h∙z−∂h∙PhzonΓ(t), together with the usual elliptic estimate, for the solution $$w\in H^2({\it{\Gamma}}(t))$$,   ‖w‖H2(Γ(t))≤c‖∂h∙z−∂h∙Phz‖L2(Γ(t)); again, $$c$$ is independent of $$t$$ and $$h$$. Following the proof of Dziuk & Elliott (2013b, Theorem 6.2), let us first bound   −b(vh;z−Phz,Ih(1)w) =b(vh;z−Phz,w−Ih(1)w)−b(vh;z−Phz,w) ≤chk‖z‖Hk+1(Γ) ch‖w‖H2(Γ)−b(vh;z−Phz,w) =chk+1‖z‖Hk+1(Γ)‖w‖H2(Γ)+b(v;z−Phz,w) +b(v;z−Phz,w)−b(vh;z−Phz,w), where again $$I_h^{(1)}$$ denotes the linear finite element interpolation operator on $${\it{\Gamma}}_h^k(t)$$. The pair in the last line can be estimated, using Lemma 5.4, by   b(v;z−Phz,w)−b(vh;z−Phz,w)≤ ∫Γ(t)|B(v)−B(vh)||∇Γ(z−Phz)||∇Γw|≤ ch2k‖z‖Hk+1(Γ)‖w‖H1(Γ). Finally, for the remaining term, the proof of Dziuk & Elliott (2013b, Theorem 6.2) yields   b(v;z−Phz,w)≥−c‖z−Phz‖L2(Γ)‖w‖H2(Γ)≥−chk+1‖z‖Hk+1(Γ)‖w‖H2(Γ). For the other bilinear form in (6.3), we have   g(vh;z−Phz,φh)≤chk+1‖z‖Hk+1(Γ)‖w‖H1(Γ). The combination of all these estimates with (6.4) yields   a∗(∂h∙z−∂h∙Phz,Ihw)≤chk+1‖z‖Hk+1(Γ)‖w‖H2(Γ). By testing the above elliptic weak problem with $$z-\mathscr{P}_h z$$, and using the above bound and the gradient estimate from (a), we obtain   ‖∂h∙z−∂h∙Phz‖L2(Γ)2= a∗(∂h∙z−∂h∙Phz,w)= a∗(∂h∙z−∂h∙Phz,w−Ih(1)w)+a∗(∂h∙z−∂h∙Phz,Ih(1)w)≤ chk(‖z‖Hk+1(Γ)+h‖∂∙z‖H1(Γ)) ch‖w‖H2(Γ(t))+chk+1‖z‖Hk+1(Γ)‖w‖H2(Γ)≤ chk+1(‖z‖Hk+1(Γ)+h‖∂∙z‖H1(Γ))‖w‖H2(Γ). For $$\ell>1,$$ the proof is analogous. □ 7. Error bounds for the semidiscretization and full discretization 7.1 Convergence proof for the semidiscretization By combining the error estimates in the Ritz map and in its material derivatives and the geometric results of Section 5, we prove convergence of the high-order ESFEM semidiscretization. Proof of Theorem 4.1. The result is simply shown by repeating the arguments of Dziuk & Elliott (2013b, Section 7) for our setting but using the high-order versions for all results: geometric estimates Lemma 5.2, perturbation estimates of bilinear forms Lemma 5.6 and Ritz map error estimates Theorems 6.3 and 6.4. □ 7.2 Convergence proof for the full discretization 7.2.1 Bound of the semidiscrete residual We follow the approach of Mansour (2013, Section 8.1) and Lubich et al. (2013, Section 5) by defining the ESFEM residual $$R_h(\cdot,t) = \sum_{j=1}^N r_j(t) \chi_j(\cdot,t)\in S_h^k(t)$$ as   ∫Γhk(t)Rhϕh=ddt∫Γhk(t)P~huϕh+∫Γhk(t)∇Γh(P~hu)⋅∇Γhϕh−∫Γhk(t)(P~hu)∂h∙ϕh, (7.1) where $$\phi_h\in S_h^k(t)$$, and $$\widetilde{\mathscr{P}}_h u(\cdot,t)$$ is the Ritz map of the smooth solution $$u$$. We now show the optimal-order $$H_h^{-1}$$-norm estimate of the residual $$R_h$$. Theorem 7.1 Let the solution $$u$$ of the parabolic problem be sufficiently smooth. Then there exist $$C>0$$ and $$h_0>0$$, such that for all $$h\leq h_0$$ and $$t\in[0,T]$$, the finite element residual $$R_h$$ of the Ritz map is bounded as   ‖Rh‖Hh−1(Γh(t))≤Chk+1. Proof. (a) We start by applying the discrete transport property to the residual equation (7.1):   mh(Rh,ϕh)= ddtmh(P~hu,ϕh)+ah(P~hu,ϕh)−mh(P~hu,∂h∙ϕh)= mh(∂h∙P~hu,ϕh)+ah(P~hu,ϕh)+gh(Vh;P~hu,ϕh). (b) We continue by the transport property with discrete material derivatives from Lemma 5.5 but for the weak form, with $$\varphi:=\varphi_h=(\phi_h)^l$$:   0= ddtm(u,φh)+a(u,φh)−m(u,∂∙φh)= m(∂h∙u,φh)+a(u,φh)+g(vh;u,φh)+m(u,∂h∙φh−∂∙φh). (c) Subtraction of the two equations, using the definition of the Ritz map (6.1) and using that   ∂h∙φh−∂∙φh=(vh−v)⋅∇Γφh holds, we obtain   mh(Rh,ϕh) =mh(∂h∙P~hu,ϕh)−m(∂h∙u,φh) +gh(Vh;P~hu,ϕh)−g(vh;u,φh) +m(u,φh)−mh(P~hu,ϕh) +m(u,(vh−v)⋅∇Γφh). All the pairs can be easily estimated separately as $$c h^{k+1} \|\varphi_h\|_{H^1({\it{\Gamma}}(t))}$$ by combining the geometric perturbation estimates of Lemma 5.6, the velocity estimate of Lemma 5.4 and the error estimates of the Ritz map from Theorems 6.3 and 6.4. The proof is finished using the definition of the $$H_h^{-1}$$-norm and the equivalence of norms Lemma 5.1. □ 7.2.2 Proof of Theorem 4.4. Using the error estimate for the BDF methods Theorem 4.2 and using the bounds for the semidiscrete residual Theorem 7.1, we give here a proof for the fully discrete error estimates of Theorem 4.4. Proof of Theorem 4.4. The global error is decomposed into two parts,   uhn−u(⋅,tn)=(uhn−(Phu)(⋅,tn))+((Phu)(⋅,tn)−u(⋅,tn)), and then the terms are estimated separately by results from above. The first term is estimated, analogously to Thomée (2006) or exactly as in Lubich et al. (2013) and Mansour (2013) as follows. The vectors collecting the nodal values of the error $$u_h^n - (\mathscr{P}_h u)(\cdot,t_n)$$ satisfy the fully discrete problem (4.1) perturbed by the semidiscrete residual (7.1) (cf. Mansour, 2013). Then applying results for BDF methods Theorem 4.2, together with the residual bound Theorem 7.1 (and by the assumption on the initial value approximation), gives the desired bound $${\mathscr {O}}(h^{k+1} + \tau^p)$$. The second term is directly estimated by the error estimates for the Ritz map and for its material derivatives, Theorems 6.3 and 6.4. □ Remark 7.2 In Remark 4.5, we noted that the analogous fully discrete results can be shown for algebraically stable implicit Runge–Kutta methods. To be more precise, for the Runge–Kutta analogue, instead of Theorem 4.2 one has to use the error bounds of Dziuk et al. (2012, Theorem 5.1), but otherwise the proof above remains the same. Remark 7.3 The statement in the introduction on various extensions is supported by the following facts. The first two points in the introduction, geometric errors and perturbation errors of the bilinear forms (Sections 5.1 and 5.4), and the basic high-order evolving surface finite element setting (Section 3.1), are independent of the considered problem. The velocity estimates (Section 5.3) depend only on the surface evolution as well. Furthermore, these general high-order results could be used in analogous ways, as their linear counterparts in the cited articles; or missing ones can be easily shown based on the ideas presented here (e.g., involving the additional term for ALE maps (Elliott & Venkataraman, 2014; Kovács & Power Guerra, 2014). The proof of the error bounds of a modified Ritz map (although perhaps defined slightly differently; see, e.g., wave equations Lubich & Mansour, 2015, or quasilinear problems Kovács & Power Guerra, 2016), rely on these geometric approximation results and some general techniques, such as the almost Galerkin orthogonality and the Aubin–Nitsche trick. Naturally, the estimates for the semidiscrete residual greatly depend on the problem itself, but in the above-mentioned articles, similar (if not the same) ideas and techniques were used. Finally, convergence of time discretizations, proved by energy estimates, carry over to high-order discretizations of various problems, from Dziuk et al. (2012), Lubich et al. (2013), Kovács & Power Guerra (2014), Lubich & Mansour (2015) and Kovács & Power Guerra (2016). Remark 7.4 Theorem 4.2 requires sufficient temporal regularity of the Ritz map. By having sufficient regularity of the solution and the surface evolution on $$[0,T]$$, using Theorem 6.4 and equivalence of norms, we obtain   ‖(∂h∙)(ℓ)(P~hu)(⋅,t)‖L2(Γh(t))≤ c‖(∂∙)(ℓ)(u−Phu)(⋅,t)‖L2(Γ(t))+c‖(∂∙)(ℓ)u‖L2(Γ(t))≤ c(cℓhk+1+1)∑j=0ℓ‖(∂∙)(j)u‖Hk+1(Γ(t)). Here, the $$H^{k+1}({\it{\Gamma}}(t))$$-norm on the right-hand side could be replaced by $$H^{2}({\it{\Gamma}}(t))$$ (also the power in $$h^{k+1}$$ would be reduced to $$2$$) by modifiying the proofs of Theorems 6.3 and 6.4, using the linear interpolation $$I_h^{(1)}$$ on $${\it{\Gamma}}_h^k(t)$$ instead of $$I_h$$. Alternatively, a weaker condition can be obtained in the following way. By having sufficient regularity of the solution at $$t=0$$ and having smooth evolution of the surface, by repeating the proof of Dziuk et al. (2012, Theorem 9.1) for the Ritz map instead of the semidiscrete solution, the following estimate is obtained:   supt∈[0,T]‖(∂h∙)(ℓ)P~hu(⋅,t)‖L2(Γhk(t))+∫0T‖∇Γh((∂h∙)(ℓ)P~hu(⋅,s))‖L2(Γhk(s))ds≤c∑j=1ℓ‖(∂h∙)(j)u(⋅,0)‖L2(Γhk(0)). The result could even be obtained directly by using Dziuk et al. (2012, Theorem 9.1) and modifying our results using the semidiscrete solution instead of the Ritz map. 8. Implementation The implementation of the high-order ESFEM code follows the typical method of finite element mass matrix, stiffness matrix and load vector assembly, mixed with techniques from isoparametric FEM theory. This was also used for linear ESFEM (see Dziuk & Elliott, 2007, Section 7.2). Similarly to the linear case, a curved element $$E_h$$ of the $$k$$-order interpolation surface $${\it{\Gamma}}_h^k(t)$$ is parameterized over the reference triangle $$E_0$$, chosen to be the unit simplex. Then the polynomial map of degree $$k$$ between $$E_h$$ and $$E_0$$ is used to compute the local matrices. All the computations are done on the reference element, using the Dunavant quadrature rule (see Dunavant, 1985). Then the local values are summed to their correct places in the global matrices. In a typical case, the surface is evolved by solving a series of ODEs, hence only the initial mesh is created based on $${\it{\Gamma}}(0)$$. Naturally, the problem of velocity-based grid distortion is still present. Possible ways to overcome this are methods using the DeTurck trick (see Elliott & Fritz, 2016) or using ALE finite elements (see Elliott & Venkataraman, 2014; Kovács & Power Guerra, 2014). 9. Numerical experiments We performed various numerical experiments with quadratic approximation of the surface $${\it{\Gamma}}(t)$$ and using quadratic ESFEM to illustrate our theoretical results. 9.1 Example 1: Parabolic problem on a stationary surface Let us briefly report on numerical experiments for parabolic problems on a stationary surface, as a benchmark problem. Let $${\it{\Gamma}}\subset\mathbb{R}^3$$ be the unit sphere, and let us consider the parabolic surface PDE   ∂tu−ΔΓu=f, with given initial value and inhomogeneity $$f$$ chosen such that the solution is $$u(x,t)=e^{-6t}x_1x_2$$. Let $$(\mathscr{T}_k)_{k=1,2,\dotsc,n}$$ and $$(\tau_k)_{k=1,2,\dotsc,n}$$ be series of meshes and time steps, respectively, such that $$2 h_k = h_{k-1}$$ and $$2 \tau_k = \tau_{k-1}$$, with $$h_1=\sqrt2$$ and $$\tau_1=0.2$$. For each mesh $$\mathscr{T}_k$$ with corresponding step size $$\tau_k$$, we numerically solve the surface PDE using second-order ESFEM combined with third-order BDF methods. Then, by $$e_k$$ we denote the error corresponding to the mesh $$\mathscr{T}_k$$ and step size $$\tau_k$$. We then compute the errors between the lifted numerical solution and the exact solution using the following norm and seminorm:   L∞(L2): max1≤n≤N‖uhn−u(⋅,tn)‖L2(Γ(tn)),L2(H1): (τ∑n=1N‖∇Γ(tn)(uhn−u(⋅,tn))‖L2(Γ(tn))2)1/2. Using the above norms, the experimental order of convergence rates (EOCs) are computed by   EOCk=ln⁡(ek/ek−1)ln⁡(2)(k=2,3,…,n). In Table 1, we report on the EOCs for the second-order ESFEM coupled with the BDF3 method; theoretically, we expect EOC$$\approx3$$ in the $$L^\infty(L^2)$$-norm, and EOC$$\approx2$$ in the $$L^2(H^1)$$-seminorm. Table 1. Errors and EOCs in the$$L^\infty(L^2)$$- and$$L^2(H^1)$$-norms for the stationary problem Level  dof  $$L^\infty(L^2)$$  EOC  $$L^2(H^1)$$  EOC  $$1$$  $$6$$  $$5.3113\cdot 10^{-3}$$  —  $$8.0694\cdot 10^{-3}$$  —  $$2$$  $$18$$  $$2.9257\cdot 10^{-3}$$  $$0.9511$$  $$4.3162\cdot 10^{-3}$$  $$0.99805$$  $$3$$  $$66$$  $$9.2303\cdot 10^{-4}$$  $$1.7122$$  $$1.9050\cdot 10^{-3}$$  $$1.2139$$  $$4$$  $$258$$  $$1.7285\cdot 10^{-4}$$  $$2.4338$$  $$5.4403\cdot 10^{-4}$$  $$1.8207$$  $$5$$  $$1026$$  $$2.6463\cdot 10^{-5}$$  $$2.7124$$  $$1.2265\cdot 10^{-4}$$  $$2.1529$$  $$6$$  $$4098$$  $$3.6845\cdot 10^{-6}$$  $$2.8457$$  $$2.4772\cdot 10^{-5}$$  $$2.3088$$  Level  dof  $$L^\infty(L^2)$$  EOC  $$L^2(H^1)$$  EOC  $$1$$  $$6$$  $$5.3113\cdot 10^{-3}$$  —  $$8.0694\cdot 10^{-3}$$  —  $$2$$  $$18$$  $$2.9257\cdot 10^{-3}$$  $$0.9511$$  $$4.3162\cdot 10^{-3}$$  $$0.99805$$  $$3$$  $$66$$  $$9.2303\cdot 10^{-4}$$  $$1.7122$$  $$1.9050\cdot 10^{-3}$$  $$1.2139$$  $$4$$  $$258$$  $$1.7285\cdot 10^{-4}$$  $$2.4338$$  $$5.4403\cdot 10^{-4}$$  $$1.8207$$  $$5$$  $$1026$$  $$2.6463\cdot 10^{-5}$$  $$2.7124$$  $$1.2265\cdot 10^{-4}$$  $$2.1529$$  $$6$$  $$4098$$  $$3.6845\cdot 10^{-6}$$  $$2.8457$$  $$2.4772\cdot 10^{-5}$$  $$2.3088$$  In Figs 1 and 2, we report on the errors   ‖u(⋅,Nτ)−uhN‖L2(Γ)and‖∇Γ(u(⋅,Nτ)−uhN)‖L2(Γ), at time $$N\tau=1$$. The logarithmic plots show the errors against the mesh width $$h$$ (in Fig. 1) and against time-step size $$\tau$$ (in Fig. 2). Fig. 1. View largeDownload slide Spatial convergence of the BDF3 / quadratic SFEM discretization for the stationary surface PDE. Fig. 1. View largeDownload slide Spatial convergence of the BDF3 / quadratic SFEM discretization for the stationary surface PDE. Fig. 2. View largeDownload slide Temporal convergence of the BDF3/quadratic SFEM discretization for the stationary surface PDE. Fig. 2. View largeDownload slide Temporal convergence of the BDF3/quadratic SFEM discretization for the stationary surface PDE. The different lines correspond to different time-step sizes and to different mesh refinements, respectively in Figs 1 and 2. In both figures we can observe two regions. In Fig. 1, a region where the spatial discretization error dominates, matching the convergence rates of our theoretical results, and a region, with fine meshes, where the time discretization error dominates (the error curves flatten out). In Fig. 2, the same description applies but with reversed roles. First, the time discretization error dominates, while for smaller step sizes the spatial error dominates. The convergence in space (Fig. 1) and in time (Fig. 2) can both be observed to be nicely in agreement with the theoretical results (note the reference lines). 9.2 Example 2: Evolving surface parabolic problem In the following experiment, we consider the parabolic problem (2.1) on the evolving surface given by   Γ(t)={x∈R3 | a(t)−1x12+x22+x32−1=0}, where $$a(t)=1+\frac{1}{4}\sin(2\pi t)$$ (see, e.g., Dziuk & Elliott, 2007; Dziuk et al., 2012; Mansour, 2013), with given initial value and inhomogeneity $$f$$ chosen such that the solution is $$u(x,t)=e^{-6t}x_1x_2$$. Similarly to the stationary surface case, we again report on the experimental orders of convergence and similar spatial and temporal convergence plots. They are all produced exactly as described above. The EOCs for the evolving surface problem solved with BDF method of order $$3$$ and evolving surface finite elements of second order can be seen in Table 2. Table 2. Errors and EOCs in the$$L^\infty(L^2)$$- and$$L^2(H^1)$$-norms for the evolving surface problem Level  dof  $$L^\infty(L^2)$$  EOC  $$L^2(H^1)$$  EOC  1  6  $$5.7898\cdot 10^{-3}$$  —  $$8.4446\cdot 10^{-3}$$  —  2  18  $$9.5840\cdot 10^{-4}$$  2.8688  $$1.8173\cdot 10^{-3}$$  2.4504  3  66  $$4.0725\cdot 10^{-4}$$  1.2704  $$1.6549\cdot 10^{-3}$$  0.13895  4  258  $$9.1096\cdot 10^{-5}$$  2.1755  $$5.4513\cdot 10^{-4}$$  1.6133  5  1026  $$1.3847\cdot 10^{-5}$$  2.7226  $$1.2774\cdot 10^{-4}$$  2.0971  6  4098  $$1.9534\cdot 10^{-6}$$  2.8267  $$2.6135\cdot 10^{-5}$$  2.2901  Level  dof  $$L^\infty(L^2)$$  EOC  $$L^2(H^1)$$  EOC  1  6  $$5.7898\cdot 10^{-3}$$  —  $$8.4446\cdot 10^{-3}$$  —  2  18  $$9.5840\cdot 10^{-4}$$  2.8688  $$1.8173\cdot 10^{-3}$$  2.4504  3  66  $$4.0725\cdot 10^{-4}$$  1.2704  $$1.6549\cdot 10^{-3}$$  0.13895  4  258  $$9.1096\cdot 10^{-5}$$  2.1755  $$5.4513\cdot 10^{-4}$$  1.6133  5  1026  $$1.3847\cdot 10^{-5}$$  2.7226  $$1.2774\cdot 10^{-4}$$  2.0971  6  4098  $$1.9534\cdot 10^{-6}$$  2.8267  $$2.6135\cdot 10^{-5}$$  2.2901  The errors at time $$N\tau=1$$ in different norms can be seen in the following plots: the different lines again correspond to different time-step sizes and to different mesh refinements in Figs 3 and 4, respectively. Fig. 3. View largeDownload slide Spatial convergence of the BDF3/quadratic ESFEM discretization for the evolving surface PDE. Fig. 3. View largeDownload slide Spatial convergence of the BDF3/quadratic ESFEM discretization for the evolving surface PDE. Fig. 4. View largeDownload slide Temporal convergence of the BDF3/quadratic ESFEM discretization for the evolving surface PDE. Fig. 4. View largeDownload slide Temporal convergence of the BDF3/quadratic ESFEM discretization for the evolving surface PDE. Acknowledgements The author is grateful for the valuable discussions with Prof. Christian Lubich during the preparation of the manuscript. We would also like to thank the anonymous referees for their extremely constructive suggestions and comments. Funding Deutsche Forschungsgemeinschaft, SFB 1173 to B.K. Appendix. Proof of the geometric approximation results of Lemma 5.2 For clarity we recall Lemma 5.2. Lemma A.1 For $${\it{\Gamma}}_h^k(t)$$ and $${\it{\Gamma}}(t)$$ as above, for $$h\leq h_0$$ with a sufficiently small $$h_0>0$$, we have the geometric approximation estimates:    ‖d‖L∞(Γh(t))≤chk+1, ‖1−δh‖L∞(Γh(t))≤chk+1, ‖n−nh‖L∞(Γh(t))≤chk, and ‖Id−δhQh‖L∞(Γh(t))≤chk+1, ‖(∂h∙)(ℓ)d‖L∞(Γh(t))≤chk+1, ‖(∂h∙)(ℓ)δh‖L∞(Γh(t))≤chk+1, ‖Pr((∂h∙)(ℓ)Qh)Pr‖L∞(Γh)≤chk+1, with constants depending only on $$\mathscr{G}_T$$, but not on $$h$$ or $$t$$. Proof of Lemma 5.2. The proofs follow the ideas of Dziuk & Elliott (2007, Lemma 5.1) and Mansour (2013, Lemma 6.1), in combination with the ideas and techniques of the proof of Demlow (2009, Proposition 2.3). Let $$E=E(t) \subset {\it{\Gamma}}_h^k(t)$$ be an element of the discrete surface. By $$I_h$$ we denote the $$k$$-order interpolation operator of Section 5.2. (i) Since the nodes of $$E$$ lie on the exact surface $${\it{\Gamma}}(t)$$ we have that the interpolate $$\widetilde I_h d$$ vanishes on $$E$$. Then by using standard interpolation estimates (from Lemma 5.3 or from Brenner & Scott, 2007) we obtain   ‖d‖L∞(E)=‖d−I~hd‖L∞(E)≤chk+1‖d‖Wk+1,∞(E)≤chk+1‖d‖Ck,1(GT). Higher-order norm estimates are shown analogously:   ‖p−pk‖Wi,∞≤chk+1−i. (A.1) (iii) For the normal vector estimate,   |n(x^)−nhk(x^)|≤ |n(pk(x~))−n(p(x~))|+|n(p(x~))−nhk(pk(x~))|≤ c(GT)hk+1+c(GT)hk, (A.2) where for the last estimate we used the smoothness of $$\mathscr{G}_T$$ and the above bound on $$d$$, and the bounds (A.1) and the Gram–Schmidt orthonormalization algorithm (cf. Demlow, 2009). (ii) The second estimate is shown by recalling, from Demlow (2009, (2.10)), that (for a fixed $$t\in[0,T]$$)   δhk(x)=n(x)⋅nhk(x)∏j=1m(1−d(x,t)Kj(x))(x∈Γhk(t)), (A.3) where $$K_j(x) = {\kappa_j(p(x,t))}/{1+d(x,t)\kappa_j(p(x,t))}$$ with $$\kappa_j$$ being the principle curvatures; cf. Demlow (2009) Then, following the proof of Demlow (2009, Proposition 4.1), using $$\|d\|_{L^\infty} = {\mathscr {O}}(h^{k+1})$$ and (A.2), we obtain   |1−δhk|≤ c(GT)hk+1+c(GT)|1−n⋅nhk|≤ c(GT)hk+1+c(GT)|n−nhk|2≤c(GT)hk+1. (A.4) (iv) To show the fourth estimate, we use the idea of the linear ESFEM case. Using the previous estimates and (5.1) the definition of $$Q_h$$,   |Id−δhQh|≤|Pr−PrPrhPr|+chk+1. Then, for an arbitrary unit vector $$z$$,   |(Pr−PrPrhPr)z|=|z⋅(nh−(nh⋅n)n)(nh−(nh⋅n)n)|≤ch2k, where the estimate follows, using (A.2), from   |nh−(nh⋅n)n|≤ |(n⋅n)nh−(nh⋅n)nh|+|(nh⋅n)nh−(nh⋅n)n|≤ |((n−nh)⋅n)nh|+|(nh⋅n)(nh−n)|≤chk. See also the proof of Demlow (2009, Proposition 4.1). The proofs of the estimates with material derivatives are similar to their nondifferentiated versions; we follow the ideas of Mansour (2013). (v) Again, since $${\left( {\partial _h^ \bullet } \right)^{(l)}} d$$ vanishes at the nodes of $$E$$, hence the interpolant $$\widetilde I_h {\left( {\partial _h^ \bullet } \right)^{(l)}} d$$ vanishes on $$E$$ completely. Again by interpolation estimates we obtain   ‖(∂h∙)(l)d‖L∞(E)=‖(∂h∙)(l)d−I~h(∂h∙)(l)d‖L∞(E)≤chk+1‖(∂h∙)(l)d‖Wk+1,∞(E)≤chk+1‖(∂h∙)(l)d‖Ck,1(GT). (vi) The sixth estimate is shown by taking the material derivative of (A.3) and by a similar argument to (A.4),   ‖(∂h∙)(ℓ)δh‖L∞(Γh(t))≤ c(GT)hk+1+c(GT)|∂h∙(n−nhk)|2≤c(GT)hk+1, where the last inequality follows by the chain rule and since the combination of (A.1) and (A.2), together with $$\mathrm{n}=\nabla d$$, yields $$|\partial^{\bullet}_h(\mathrm{n}-\mathrm{n}_h^k)| = {\mathscr {O}}(h^k)$$. (vii) Let us show the last estimate for $$\ell=1$$. The higher-order version follows by a similar argument recursively. It is clear from the definition of $$Q_h$$ (see (5.1)) that with a smooth remainder function $$R$$ we can write   Qh=1δhPrPrhPr+dR(δh,Pr,Prh,H). Now using the facts $$\|d\| , \|\partial^{\bullet}_h d\| ={\mathscr {O}}(h^{k+1})$$, $$\delta_h=1+{\mathscr {O}}(h^{k+1})$$ and $$\|\partial^{\bullet}_h \delta_h\|={\mathscr {O}}(h^{k+1})$$ we bound   Pr(∂h∙Qh)Pr=Pr∂h∙(PrPrhPr)Pr+O(hk+1). The first term here is estimated separately. Using $$\partial^{\bullet}_h \mathrm{n} \cdot \mathrm{n} = 0$$, in Mansour (2013, equation 6.4), it is shown that   Pr∂h∙(PrPrhPr)Pr=Pr∂h∙(PrPrhPr−Pr)Pr=−Pr∂h∙(Pr nhnhTPr)Pr. Since the projections are bounded, we need only the bounds   |Pr nh|=|nh−(n⋅nh)n|≤ chk,|∂h∙(Pr nh)|=|∂h∙(nh−(n⋅nh)n)|≤ chk, where the first inequality has been shown above, while the second follows by the chain rule and by $$|\partial^{\bullet}_h(\mathrm{n}-\mathrm{n}_h^k)| = {\mathscr {O}}(h^k)$$ proved above, together with the boundedness of $$\partial^{\bullet}_h \mathrm{n}$$ and $$\partial^{\bullet}_h \mathrm{n}_h$$ (which follows from velocity approximation and the smoothness of $$\mathscr{G}_T$$ for the first term, while for the latter by the same arguments and an additional triangle inequality). □ References Antonietti P., Dedner A., Madhavan P., Stangalino S., Stinner B. & Verani M. ( 2015) High order discontinuous Galerkin methods for elliptic problems on surfaces. SIAM J. Num. Anal , 53, 1145– 1171. Google Scholar CrossRef Search ADS   Aubin T. ( 1998) Some Nonlinear Problems in Riemannian Geometry,  1st edn. Berlin: Springer. Google Scholar CrossRef Search ADS   Brenner S. & Scott R. ( 2007) The Mathematical Theory of Finite Element Methods , vol. 15. New York: Springer. Burman E., Hansbo P. & Larson M. ( 2015) A stabilized cut finite element method for partial differential equations on surfaces: The Laplace–Beltrami operator. Comput. Methods Appl. Mech. Eng. , 285, 188– 207. Google Scholar CrossRef Search ADS   Dedner A., Madhavan P. & Stinner B. ( 2013) Analysis of the discontinuous Galerkin method for elliptic problems on surfaces. IMA J. Numer. Anal. , 33, 952– 973. Google Scholar CrossRef Search ADS   Demlow A. ( 2009) Higher-order finite elment methods and pointwise error estimates for elliptic problems on surfaces. SIAM J. Numer. Anal. , 47, 805– 807. Google Scholar CrossRef Search ADS   Dunavant D. ( 1985) High degree efficient symmetrical Gaussian quadrature rules for the triangle. Int. J. Numer. Meth. Eng. , 21, 1129– 1148. Google Scholar CrossRef Search ADS   Dziuk G. ( 1988) Finite elements for the Beltrami operator on arbitrary surfaces. Partial Differential Equations and Calculus of Variations . Lecture Notes in Mathematics, 1357. Berlin: Springer. pp. 142– 155. Google Scholar CrossRef Search ADS   Dziuk G. & Elliott C. ( 2007) Finite elements on evolving surfaces. IMA J. Numer. Anal. , 27, 262– 292. Google Scholar CrossRef Search ADS   Dziuk G. & Elliott C. ( 2013a) Finite element methods for surface PDEs. Acta Numerica , 22, 289– 396. Google Scholar CrossRef Search ADS   Dziuk G. & Elliott C. ( 2013b) $$L^2$$-estimates for the evolving surface finite element method. Math. Comp. , 82, 1– 24. Google Scholar CrossRef Search ADS   Dziuk G., Lubich C. & Mansour D. ( 2012) Runge–Kutta time discretization of parabolic differential equations on evolving surfaces. IMA J. Numer. Anal. , 32, 394– 416. Google Scholar CrossRef Search ADS   Elliott C. & Fritz H. ( 2016) On approximations of the curve shortening flow and of the mean curvature flow based on the DeTurck trick. To appear i. IMA J. Numer. Anal. , doi: 10.1093/imanum/drw020. Elliott C. M. & Ranner T. ( 2013) Finite element analysis for a coupled bulk–surface partial differential equation. IMA J. Numer. Anal. , 33, 377– 402. Google Scholar CrossRef Search ADS   Elliott C. M. & Ranner T. ( 2015) Evolving surface finite element method for the Cahn–Hilliard equation. Numer. Math. , 129, 483– 534. Google Scholar CrossRef Search ADS   Elliott C. M. & Venkataraman C. ( 2014) Error analysis for an ALE evolving surface finite element method. Numer. Methods Partial Differential Equations , 31, 459– 499. Google Scholar CrossRef Search ADS   Grande J. & Reusken A. ( 2014) A higher order finite element method for partial differential equations on surfaces. SIAM J. Numer. Anal. , 54, 388– 414. Google Scholar CrossRef Search ADS   Hairer E. & Wanner G. ( 1996) Solving Ordinary Differential Equations II: Stiff and Differential–Algebraic Problems , 2nd edn. Berlin: Springer. Kovács B., Li B., Lubich Ch. & Power Guerra C. A. ( 2017) Convergence of finite elements on an evolving surfacedriven by diffusion on the surface, arXiv:1607.07170v2. Kovács B. & Power Guerra C. A. ( 2014) Higher–order time discretizations with ALE finite elements for parabolic problems on evolving surfaces. to appear i. IMA J. Numer. Anal. , arXiv:1410.0486. Kovács B. & Power Guerra C. A. ( 2016) Error analysis for full discretizations of quasilinear parabolic problems on evolving surfaces. Numer. Methods Partial Differential Equations , 32, 1200– 1231. Google Scholar CrossRef Search ADS   Lubich C. & Mansour D. ( 2015) Variational discretization of wave equations on evolving surfaces. Math. Comp. , 84, 513– 542. Google Scholar CrossRef Search ADS   Lubich C., Mansour D. & Venkataraman C. ( 2013) Backward difference time discretization of parabolic differential equations on evolving surfaces. IMA J. Numer. Anal. , 33, 1365– 1385. Google Scholar CrossRef Search ADS   Mansour D. ( 2013) Numerical analysis of partial differential equations on evolving surfaces. Ph.D. Thesis , Universität Tübingen. http://hdl.handle.net/10900/49925. Accessed 1 March 2017. Olshanskii M., Reusken A. & Grande J. ( 2009) A finite element method for elliptic equations on surfaces. SIAM J. Numer. Anal. , 47, 3339– 3358. Google Scholar CrossRef Search ADS   Reusken A. ( 2015) Analysis of trace finite element methods for surface partial differential equations. IMA J. Numer. Anal. , 35, 1568– 1590. Google Scholar CrossRef Search ADS   Thomée V. ( 2006) Galerkin Finite Element Methods for Parabolic Problems , 2nd edn. Berlin: Springer. © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

### Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Jan 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### Monthly Plan • Read unlimited articles • Personalized recommendations • No expiration • Print 20 pages per month • 20% off on PDF purchases • Organize your research • Get updates on your journals and topic searches$49/month

14-day Free Trial

Best Deal — 39% off

### Annual Plan

• All the features of the Professional Plan, but for 39% off!
• Billed annually
• No expiration
• For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588$360/year

billed annually

14-day Free Trial