An implicitly extended Crank-Nicolson scheme for the heat equation on a time-dependent domain

We consider a time-stepping scheme of Crank-Nicolson type for the heat equation on a moving domain in Eulerian coordinates. As the spatial domain varies between subsequent time steps, an extension of the solution from the previous time step is required. Following Lehrenfeld & Olskanskii [ESAIM: M2AN, 53(2): 585-614, 2019], we apply an implicit extension based on so-called ghost-penalty terms. For spatial discretisation, a cut ﬁnite element method is used. We derive a complete a priori error analysis in space and time, which shows in particular second-order convergence in time under a parabolic CFL condition. Finally, we present numerical results in two and three space dimensions that conﬁrm the analytical estimates, even for much larger time steps.


Introduction
Partial differential equations (PDEs) posed on moving domains are significant in many areas of science and engineering.They arise for example in flow problems around moving structures, such as pumps [4], wind or water turbines [54], within moving objects [15], or as sub-problems in fluidstructure interactions or multiphase flows.Fluid-structure interactions arise in aerodynamical applications like flow around airplanes or parachutes [58], in biomedical problems such as blood flow through the cardiovascular system [53,60,25] or the airflow within the respiratory system [62] and even in tribological applications [46].Multiphase problems include for instance gas-liquid and particle-laden gas flows [36,19,45], rising bubbles [43], droplets in microfluidic devices [17] or the simulation of tumor growth [34].For further details and applications we refer to the textbooks [55,3] and [35], respectively.
In this article we consider the time discretisation of a parabolic model problem (namely the heat equation) which is posed on a moving domain Ω(t) ⊂ R d (d = 2, 3) that evolves smoothly in time for t ∈ I = [0, t max ]: In literature, two major numerical approaches can be found for the simulation of partial differential equations on moving domains: the Arbitrary Lagrangian Eulerian (ALE) [21,22] approach, where the equations are transformed to an arbitrary reference domain which is independent of time, and Eulerian approaches, where the equations are solved in the time-dependent Eulerian framework [23,27,48].
The ALE approach is a popular technique for the numerical simulation of PDEs on moving domains, in particular for flow problems [42,41].For details, we refer to the textbooks [3,55] and reference cited therein.Convection-diffusion problems on moving domains were, for example, solved in [33,57] using a stabilised ALE method.The ALE approach is very attractive in the case of moderate domain movements, but shows problems when the shape of the domain changes significantly in time.In particular, topology changes of Ω(t), as occurring for example in contact problems or considering the separation or union of bubbles can not be modelled by means of an ALE approach [55,17,8,9].Another example of extreme variations of Ω(t) are so-called fingering phenomena, which can be frequently observed in multi-phase flows or even for tumor growth [34].
In such cases, a numerical approach that discretises the equations directly in the moving Eulerian coordinate framework is preferable.The Eulerian framework is also the coordinate framework, which is typically used to model flow problems and consequently, in multi-phase flows [35] and fluid-structure interactions with large displacements [23,10,27].However, as the domains Ω(t) to be discretised vary with time t, additional difficulties arise concerning a proper and accurate discretisation, both in space and in time.
In recent decades, a great amount of works have been contributed concerning the spatial discretisation of curved or moving boundaries by means of finite elements.The techniques can be categorised in fitted and unfitted finite element methods.In fitted methods, the boundary ∂Ω(t) is resolved in each time step by the finite element mesh [24,5,30].If the domain is time-dependent, this means that new meshes need to be created in each time step.Several approaches have been proposed to alleviate this issue, such as the locally fitted finite element method [30,29], which is based on a fixed coarse and a variable fine mesh.However, different issues might arise, such as anisotropic fine cells that complicate the numerical discretisation [28].
The idea in unfitted finite element methods, on the other hand, is to use the same finite element mesh for all times t, independently of the position of the boundary ∂Ω(t).A popular approach is the cut finite element method (CutFEM) [7,13,38,51,37,64], where cells of the finite element mesh are cut into parts that lie inside Ω(t) and parts outside for numerical integration.Boundary values are then incorporated weakly by means of Nitsche's method [52].The method shows similarities to the extended finite element method [20,16,32] and the generalised finite element method [2], where the finite element spaces are enriched by suitable functions to account for the position of the boundary.
Much less works can be found in literature concerning a proper time discretisation on moving domains.In the case of moving domains, standard time discretisation based on the method of lines is not directly applicable.The reason is that the domain of definition of the variables changes from time step to time step.As an example consider the finite difference discretisation of the time derivative within a variational formulation (∆t = t n − t n−1 ) The function u h (t n−1 ) is only well-defined on Ω(t n−1 ), but is needed on Ω(t n ).A possible remedy is to use characteristic-based approaches based on trajectories that follow the motion of the domain, see e.g.[40].Similar time-stepping schemes result when applying the ALE method only locally within one time step and projecting back to the original reference frame after each step [18], or based on Galerkin time discretisations with modified Galerkin spaces [31].The disadvantage of these approaches is the necessity for a projection that needs to be computed within each or after a certain number of steps.
A further alternative are space-time approaches [39,47], where a d + 1-dimensional domain is discretised.These are, however, computationally demanding, in particular within complex three-dimensional applications.The implementation of higher-dimensional discretisations and accurate quadrature formulas pose additional challenges.If a discontinuous Galerkin approach is applied in time for the test functions, the formulation decouples in certain time intervals and can be seen as an Eulerian time-stepping scheme [39,64,26].
In this work, we follow a slightly different approach first used by Schott [56] and later analysed by Lehrenfeld & Olshanskii [48].Here, the idea is to define extensions of the solution u(t n−1 ) from previous time steps to a domain Ω δ (t n−1 ) that spans at least Ω(t n ).On the finite element level these extensions can be incorporated implicitly in the time-stepping scheme by so-called ghost penalty stabilisations [6] to a sufficiently large domain.These techniques have originally been proposed to extend the coercivity of elliptic bilinear forms from the physical to the computational domain in the context of CutFEM or fictitious domain approaches [6].
Lehrenfeld & Olshanskii [48] analysed the so extended Backward Euler method in detail for a convection-diffusion problem and gave hints on how to transfer the argumentation to the secondorder backward difference scheme (BDF2).Recently, the analysis has been extended to higher order in space and time using an isoparametric finite element approach [49].In [11,61], extended BDF time-stepping schemes were applied and analysed for the non-stationary Stokes equations on moving domains.
The reason why only BDF-type time-stepping schemes have been considered in previous works, is that in these schemes spatial derivatives appear only on the "new" time step, i.e. ∇u(t n ).We will see below that the appearance of additional derivatives on u(t n−1 ) will complicate the error analysis severely.This paper gives a first step towards the analysis of time-stepping schemes that require derivatives at different time instants, such as the Crank-Nicolson method, the Fractional-step-θ-, implicit Runge-Kutta-or Adams-Bashforth schemes.
As a first step, we focus in this work on the popular Crank-Nicolson time-stepping scheme.Up to now, it has been largely open, if and under what conditions a Crank-Nicolson-type scheme can be used within an Eulerian time discretisation on moving domains.We give a detailed stability and convergence analysis.While the analysis requires a strong parabolic CFL condition of type ∆t ≤ ch 2 , our numerical results indicate that the scheme is stable also for much larger time steps.
The article is organised as follows: In Section 2, we introduce the discretisation of the model problem (1.1) in time and space.Section 3 presents a stability analysis for the fully discrete scheme using a CFL condition.In Section 4, we show a detailed a priori convergence analysis.Numerical experiments in two and three space dimensions are presented in Section 5. Section 6 summarises this article with some concluding remarks.

Discretisation
In this section, we present the numerical approximation of the model problem (1.1).We start with discretisation in time, and continue with the spatial discretisation of the resulting time-discrete formulation.

Temporal discretisation
For time discretisation, we divide the time interval of interest I = [0, t max ] in intervals I n = (t n−1 , t n ].For simplicity, we take a uniform time step ∆t = t max N and define t n = n∆t.We define the domain Ω n := Ω(t n ) with boundary Γ n := Γ(t n ) and write u n = u(t n ) for the exact solution of the continuous problem (1.1) at time t n .A δ-neighborhood of Ω(t) at time step n is chosen large enough such that (Ω n ∪ Ω n+1 ) ⊂ Ω n δ , see Figure 1.Therefore we choose δ ≥ w max ∆t, The required regularity of the domain mapping T will be ensured in Assumption 1 below.For the error analysis we will also assume the upper bound with a constant c δ > 1.Finally, we introduce the following notations for some space-time domains Now, the Crank-Nicolson method applied to (1.1) writes formally The main issue of this formulation is that u n−1 is needed on Ω n , while it is defined on Ω n−1 .Thus, we will add implicit extension operators below to define u n on Ω n δ ⊃ Ω n+1 , where it is needed in the following time step.Similarly, f n−1 might be undefined on Ω n \ Ω n−1 .If f n−1 is given analytically, it can typically be extended in a canonical way to Ω n .To cover different scenarios, we do not want to restrict the analysis in this work to a particular extension, but assume only that f n−1 is smoothly extended to Ω n .
In this article, we will use the abbreviation c to refer to a generic positive constant, which is independent of discretisation parameters (∆t, h) and the relative positions of the boundary with respect to the mesh.

Extension operator
In this part, we introduce an extension operator of the exact solution u to extend variables to larger domains, as the spatial domain evolves.We make the following assumption (see also [11,Assumption 3.2]) for the analysis of this article.
Assumption 1 The boundary of the initial domain Ω(0) is assumed to be piecewise smooth and Lipschitz, and the domain motion T (t) is a W 1,∞ -diffeomorphism for each t, that fulfills T ∈ W r,∞ ( Q), where r = max{3, m + 1} and m is the polynomial degree of the finite element space defined in the following subsection.
By using assumption 1, there exist W r,∞ -stable extension operators E n from Ω n to Ω n δ that satisfy the following analytical properties: The properties (2.3) and (2.4) are discussed in [11].In an analogous way, one can derive the estimate for the third-order time derivative in (2.5).

Spatial discretisation
For spatial discretisation, we introduce a polygonal domain D, which is chosen large enough, such that Ω n δ ⊂ D for all n.We introduce a quasi-uniform family of triangulations (T h ) h>0 of D with maximum cell size h, which will serve as background meshes.
In each time step, we extract from T h all cells of non-empty intersection with Ω n δ and define We write Ω n h,δ for the domain spanned by all cells K ∈ T n h,δ and define the following finite element space: V n,m h The set of elements that lie (at least partially) outside of Ω n−1 , but in Ω n , will be of particular interest in the analysis.The domain spanned by them will be denoted by Moreover, we introduce the following notations for the facets of T n h,δ , see Figure 2 for an illustration: • F n h,δ : the set of interior facets of T n h,δ .
• F n,int h,δ : the set of facets that belong exclusively to elements K ∈ T n h,δ that lie completely in the interior of Ω n .
• F n,cut h,δ : the set of facets that belong to some element K ∈ T n h,δ with K ∩ ∂Ω n = ∅.
Assumption 2 (CFL condition) We assume the parabolic CFL condition ∆t ≤ c CF L h 2 , where c CF L is an arbitrary constant for m = 1, while we assume c CF L sufficiently small for m > 1.
The inequality (2.1) and the CFL condition (Assumption 2) lead to Remark 2.1 The inequality (2.6) implies that the distance between ∂Ω n and ∂Ω n−1 is bounded by O(h 2 ).This implies the following property, which will be needed in the analysis below: For each cell K ∈ S n,n−1 h , there exists a path of cells and the final cell K M lies fully in the interior of Ω n .Furthermore, the number of cases, in which an element K M ⊂ Ω n is utilised as a final element among all paths, can be bounded independently of h and ∆t.

Discrete variational formulation
In the numerical approximation, the boundary condition of the discrete problem (2.7) is implemented weakly by means of Nitsche's method.Moreover, the function u n h is extended by means of a ghost penalty term g n h (•, •).In each time step n = 1, 2, . . ., N , we consider the following discrete variational formulation: Find where and We assume that all integrals in (2.8) are evaluated exactly.For a consideration of additional quadrature errors, that result when cut cells are approximated linearly in the computation of the integrals, we refer to [48].
The ghost penalty stabilization is defined by where [[•]] is the jump operator and ∂ n the exterior normal derivative.For further possibilities for the extension g n h , we refer to [48].The variant chosen here based on the jump of derivatives over edges has the advantage that it is fully consistent, in the sense that g n h (u, v) vanishes for u ∈ H m+1 (Ω n δ ).The purpose of the ghost penalty is twofold: First, it serves to extend the solution u n h implicitly to Ω n h,δ .Secondly, it ensures the discrete coercivity of the formulation (2.7) on T n h,δ .To incorporate the initial condition, we set u 0 , where E 1 u 0 is a smooth, e.g. a canonical extension, of the initial value u 0 .This corresponds to the following Ritz projection of the initial value u 0 The following lemma is the key to extend the discrete coercivity to Ω n δ : In addition, for v, w ∈ H m+1 (Ω n δ ), m ≥ 1, it holds A proof of this lemma is given in [48].
At the end of this section, we briefly show that the variational formulation (2.7) is well-posed for each n.We define the discrete energy as (2.10) and the energy norm as We will show the coercivity relation for sufficiently large γ g , γ D .The well-posedness of (2.7) follows then by standard arguments.
From the definition of the bilinear form A(u n h , 0; u n h ), we have The term − 1 2 (∂ n u n h , u n h ) ∂Ω n is estimated by means of Young's inequality, an inverse inequality and Lemma 2.2 as follows for sufficiently small > 0: For sufficiently large parameters γ D , γ g , we have This proves the coercivity (2.11).

Stability analysis
In this section, a detailed stability analysis of the discrete problem (2.7) is developed.One of the main issues in the analysis is that the discrete functions u n−1 h and ∇u n−1 h appear on Ω n in the n-th time step, whereas bounds are only available for u n−1 h Ω n−1 and ∇u n−1 h Ω n−1 from the previous time step.We start with some technical lemmas that will enable us to deal with this issue.

Auxiliary estimates
h satisfy the following inequality: ). (3.1) In the general case (m > 1), we have still . By Remark 2.1, there is a set of neighbouring cells K 2 , ..., K M , such that ) and K M lies fully in the interior of Ω n .Let e 1,2 be the edge that separates the cells K 1 and K 2 , see Figure 2.Then, using arguments from [50, Lemma 5.1] and [48, Lemma 5.2], we can deduce that (3.4) We follow this process from K 2 to K M by crossing edges e 2,3 to e M −1,M to obtain For the first term on the right-hand side of (3.5), we use an inverse inequality and the CFL condition with a sufficiently small constant c CF L As all edges e j−1,j belong to both F g,n h,δ and F g,n−1 h,δ , we obtain after summation over all cells in Using (3.2) instead of (3.3) for m = 1, the same result follows under the CFL condition ∆t ≤ c CF L h 2 for an arbitrary constant c CF L .
We note that the CFL condition is required to estimate the "mismatch" Ω n (see the following lemma).Next, we discuss how the term ∇u n−1 h Ω n can be bounded by ∇u n−1 h Ω n−1 plus further terms that can be controlled in the following stability analysis.

Lemma 3.2 Under the assumptions of Lemma 3.1 it holds for
). (3.7) Proof.By the triangle inequality, we have By means of the inequality (a + b) The statement follows by using Lemma 3.1 for the second term in (3.8) ). (3.9) Next, we provide the following Poincaré-type estimate: Proof.The proof follows the lines of [55, Lemma 4.34] and uses the fact that dist(Ω n , Ω n−1 ) ≤ c∆t.
Using this, we can derive bounds for Proof.By means of Lemma 3.3 for p = 2 we have for The statement follows by applying Lemma 3.2 to the second term in (3.12).

Stability result
Before discussing the stability result, we introduce some abbreviations for the space-time Bochner norms to simplify the mathematical expressions where m ∈ N ∪ {0} and H 0 (Ω(t)) := L 2 (Ω(t)).Now we are ready to prove the following stability result.
Theorem 3.5 (Stability) Let Assumptions 2 be valid and let f ∈ L ∞ (I, L 2 (Ω(t))), u 0 ∈ H 1 (Ω 0 ) and let the mapping T be a W 1,∞ -diffeomorphism.For sufficiently large γ g , γ D the solution We estimate the fourth term in (3.14) by means of Young's inequality with a sufficiently small > 0 followed by an inverse inequality where γ D ≥ 8/ .By using the relation 2(a + b, a) = (a + b) 2 + a 2 − b 2 for the first two terms in (3.14), we obtain that (3.15) For n > 1, we bring the terms u n−1 h Ω n and ∇u n−1 h Ω n to the domain Ω n−1 .By employing Lemmas 3.2 and 3.4, we have Inserting (3.16) into (3.15) and using 2∆t( (3.17) For n = 1, we obtain from (3.15) and the stability of the extension E 1 (3.18) Taking the sum over k = 1, 2, . . ., n, this yields for sufficiently large γ g and γ D The statement follows by means of the discrete Gronwall lemma.
Remark 3.6 (CFL condition) For m = 1, the stability result in Theorem 3.5 could be shown under the weaker CFL condition ∆t ≤ c CF L h 3/2 for sufficiently small c CF L .This is due to the possibility to use the estimate ) is constant in each cell T ).However, in the following section the stronger CFL condition ∆t ≤ c CF L h 2 with arbitrary c CF L (see Assumption 2) will be used in order to show optimal-order convergence estimates.The CFL condition is needed to estimate the term ) and Lemma 3.2.

A priori error analysis
In this section, we show an a priori error estimate for the discrete problem (2.7).We define the discretisation error as where u n := u(t n ) is assumed to be at least in H 2 (Ω n ).Further regularity assumptions on u will be made below.The error is decomposed into an interpolation error η n and a discrete error ξ n h terms defined by where I n h u n denotes the standard Lagrangian nodal interpolation of u n on T n h,δ .For n = 0 we have, by definition of u 0 h , that e 0 = u 0 − u 0 h = 0 in Ω 0 and thus, we also set η 0 = ξ 0 h = 0. We will make use of the following standard interpolation estimates for n ≥ 1

Consistency error
The exact solution u ∈ H 1 (Ω(t)) of the continuous problem (1.1) satisfies the following weak formulation: for v ∈ H 1 (Ω(t)) and the bilinear form At time t n−1 , we have where To estimate the consistency error, we will need an analogous equality for u n−1 on Ω n instead of Ω n−1 .Therefore, we define ũ as a smooth extension of the exact solution u to Q n δ .Moreover, we define a smooth extension of the source term f as follows: However, as we have allowed an arbitrary smooth extension of f to Ω n−1 δ in Section 2, this does not necessarily hold in By adding the equations (4.5) and (4.8), we obtain We note that the right-hand side in (4.9) differs from the discrete formulation (2.7) by ( Hence, (4.9) can be rewritten as where given by In the next lines, we will discuss a bound for the term E n−1 f .Lemma 4.1 Under the assumptions of Lemma 3.1, the error term E n−1 f defined in (4.11) satisfies the following estimate for v n h ∈ V n,m h and sufficiently regular u: Proof.We apply Lemma 3. .
(4.15) Inserting these estimates into (4.14) and using Young's inequality and Assumption 2, we obtain ) . (4.16) Now we are ready to estimate the consistency error related to the discrete problem (2.7).By subtracting (2.7) from (4.10), the global error term e n satisfies the equality .17) where the consistency error E n c (v h ) is given by .
(4.18)The terms I 3 and I 4 vanish due to the homogeneous boundary condition and the regularity assumption on the exact solution u n ∈ H 2 (Ω n ) and its extension E n u n ∈ H 2 (Ω n δ ).The remaining terms are estimated in the following lemma.
Proof.First, we will show a bound for the term I 1 .By following the argumentation in [59, Chapter 1], we have Using the stability of the extension operator E n given in (2.5) and the Cauchy-Schwarz inequality, we have (4.20) A bound for the second term I 2 follows in a similar way, see also [59,Chapter 1].The statement follows using Young's inequality.

Interpolation error
To derive an interpolation error estimate, we devise a discrete problem associated with the discrete error ξ n h .By definition of ξ n h (4.2) and using (4.17), we have for where the interpolation error E n I (v h ) is given by We remark that we are using smooth extensions of, e.g., η n−1 , in the following without further notice, whenever variables would be undefined on parts of Ω n .
. Based on the assumptions 1 and 2, the interpolation error is bounded by Proof.The first term in the interpolation error E n I (v h ) from (4.22) is estimated as follows: Using the stability of the extension (2.4), we obtain Next, we use interpolation estimates (4.3) and (4.4) in combination with (2.3) to deduce for k Using that, by an inverse inequality and the CFL condition ∆t ≤ c CF L h 2 , (4.27)For the Nitsche penalty term, we have (4.28) Finally, the ghost-penalty term is approximated by using The statement (4.23) follows by combining estimates (4.25)-(4.29)and using Young's inequality.

Convergence estimate
Q) be the solution of (1.1) and {u k h } n k=1 the discrete solution of (2.7), respectively.Under Assumptions 1 and 2, the discrete error term ξ n satisfies for γ g , γ D sufficiently large where , with R C and R I specified in Lemma 4.1 and 4.3, respectively.
Proof.By taking v h = 2∆tξ n h in (4.21) and using the argumentation from the stability proof (Theorem 3.5), see (3.17), we obtain that By combining results from Lemmas 4.1, 4.2 and 4.3 we have Inserting (4.32) into (4.31) and absorbing terms into the left-hand side yields for n > 1 and γ D , γ g sufficiently large As ξ 0 h = 0, we obtain for n = 1 Now, summing over k = 1, 2, . . ., n, we deduce that (4.34) Finally, we use the discrete Gronwall lemma to conclude the result.This completes the proof.
where R(u) is defined in Lemma 4.4.

Numerical results
In this section, we show numerical results in two and three space dimensions to verify the theoretical findings and the practical behaviour of the numerical method.All the numerical experiments have been obtained using the CutFEM library [7], which is based on FEniCS [1].
To verify the theoretical results, we will analyse the error terms e k = u k − u k h , k = 1, ..., n in the following norms Given the respective CFL condition (Assumption 2), Theorem 4.5 guarantees second-order convergence in time and convergence of order m in space in the L 2 -norm at the end time e n L 2 (Ω n ) and in the averaged L 2 (H 1 )-norm e L 2 (H 1 av ) .

2d example
Example 5.1 We consider a circle traveling with constant velocity w = (1, 0) towards the right.The domain is given by The data of the model example is chosen in such a way that the exact solution is u(x, y, t) = exp(−4π 2 t) cos(2πx) cos(2πy).
An illustration of the numerical solution is given in Figure 3.As background domain D, we use the unit square, i.e.D = [0, 1] 2 .The background triangulations T h are created by a uniform subdivision of the unit square into triangles and successive refinement.For each time step n, the active triangulation T n h,δ is then extracted from T h , as described in Section 2. We use δ = 4∆t.Estimated orders of convergence We will show results for different time-step sizes ∆t i = 1 /50 • 2 −i , i = 0, ..., 4 and mesh sizes h j = 1 /32 • 2 −j , j = 0, ..., 3. From the computed errors, we will estimate the temporal and spatial order of convergence.Therefore, we assume that the total error can be decomposed into a temporal and a spatial component as follows with constants c ∆t , c h and estimated orders of convergence eoc ∆t , eoc h .For a fixed mesh size h j , this relation becomes with a fixed spatial error part g h j .We will use (5.1) to estimate the order of convergence in time eoc ∆t by means of a least-squares fit of all available error values for a fixed h j to find the three parameters g h j , c ∆t and eoc ∆t .Analogously, we estimate the spatial order of convergence eoc h by a least-squares fit of the function against all available error values for a fixed time-step size ∆t i .Finally, we will also compute estimated orders of convergence for the "diagonal values", which correspond to fixing ∆t = ch for c ∈ { 32 /50, 32 /100}.Here we fit the two parameters c ∆t,h and eoc ∆t,h of the function g(ch, h) = c ∆t,h h eoc ∆t,h .
(5.3) against the computed error values.

P 1 finite elements
Firstly, we consider P 1 finite elements.We choose the numerical parameters as γ D = 1 and γ g = 10 −3 .The errors in the L 2 -norm at the end time, the L 2 (L 2 )-norm and the L 2 (H 1 av )-norm are shown in Table 1 for different ∆t and h.The estimated orders of convergence are shown, if the asymptotic standard error (computed by gnuplot [63]) was below 20%; otherwise we draw a '-'.
We observe estimated spatial convergence orders close to two in the L 2 -norms and close to one in the L 2 (H 1 av )-norm.Note that Theorem 4.5 guarantees only first-order convergence in space.We expect, however, that using a duality argument second-order convergence in space could be shown in the L 2 (L 2 )-norm, as in [11].
End-time error e n L 2 (Ω n ) h ↓ /∆t → 1 /50 The estimated temporal orders of convergence are close to two or larger in the L 2 -norms.This is in agreement with Theorem 4.5.In the L 2 (H 1 av )-norm the estimated eoc ∆t seems to be even larger than three.This has to be read carefully, however, as the spatial error part clearly dominates the overall error in this case.
Finally, the diagonal orders are around two in the L 2 (L 2 )-norm and even slightly higher in the L 2 -norm at the end time.In the L 2 (H av )-norm the spatial part is dominant and we obtain eoc ∆t,h close the one, in agreement with Theorem 4.5.

P 2 finite elements
Next, we consider Example 5.1 with P 2 finite elements.We increase the Nitsche parameter to γ D = 10, as for higher polynomial degree a larger Nitsche parameter is required, see e.g.[44].The ghost-penalty parameter γ g is still chosen as 10 −3 , but now, according to (2.9), second derivatives are included in the ghost-penalty term.The L 2 -norm errors at the end time, the L 2 (L 2 )-and the L 2 (H P 1 finite elements.The spatial orders of convergence are close to two in all norms for smaller ∆t.We note that in Theorem 4.5 second order in space has been shown for the end-time L 2 -norm and the L 2 (H 1 av )-norm.Using a duality argument, one could even hope for convergence order three in the L 2 (L 2 )-norm.We need to consider, however, that for these results the quadrature error related to a curved boundary has not been taken into account.In the CutFEM library used here, the geometry is approximated linearly in the set-up of the quadrature rule, see [7].This can lead to a reduced order of convergence, namely order 1.5 in the H 1 -norms and 2 in the L 2 -norms, see [12] for results for a CutFEM approach applied to an elliptic problem on curved domains.This reduction can also be observed in the L 2 -norm errors in the example considered here.
The estimated temporal orders of convergence are again close to two in the L 2 -norms, which confirm the estimates in Theorem 4.5.In the L 2 (H 1 av )-norm the error is still clearly dominated by the spatial part for h ≥ 1  64 .This changes, however, on the finer levels, where the temporal part gets dominant and the eoc ∆t is very close to two, in agreement with Theorem 4.5.
The spatial and temporal convergence orders are confirmed by the "diagonal" orders, which are close to two in all cases.  .The background triangulations T h are created by uniform subdivisions of D into tetrahedra and successive refinement.We use again δ = 4∆t and choose γ g = 0.1 for P 1 and γ g = 1 for P 2 finite elements, respectively and γ D = 10 in both cases.We note that in this example the quadrature error is zero, as the boundary ∂Ω k consists of plane surfaces for all k.An illustration of the numerical solution at times t = 0 and t = 1 is given in Figure 4.

P 1 finite elements
As the numerical experiments in three space dimensions are much more time-consuming compared to two dimensions, we focus on simultaneous refinement in space and time by choosing h i = ∆t i 10 = 2 −i−1 for i = 0, ..., 3. The resulting errors in the four norms introduced above are plotted in Figure 5 over the mesh size (blue curves) and compared to linear (red) and quadratic convergence (pink).We observe second-order convergence in the L 2 -norms and first-order convergence in the L 2 (H 1 av )norm.We note that for m = 1 first-order convergence in space has been shown in Theorem 4.5 in the L 2 -norm at the end time and in the L 2 (H 1 av )-norm.The numerical results indicate again that second-order convergence in space could be shown in the L 2 (L 2 )-norm using a duality argument.The first-order convergence in the L 2 (H 1 av )-norm is optimal, as the spatial error part dominates the overall error.

P 2 finite elements
In Figure 6, we illustrate the errors under simultaneous refinement (h i = ∆t i 10 ) for P 2 finite elements.We observe again (at least) second-order convergence in the L 2 -norms, in agreement with Theorem 4.5.The convergence in the L 2 (H 1 av )-norm lies between linear and quadratic convergence and decreases slightly for finer mesh sizes.This reduction in the convergence rate could be related to a violation of the CFL condition, which has been used in the error analysis.On the other hand, we did not find any stability issues in our computations.

Concluding remarks
We have analysed a Crank-Nicolson variant of the implicitly extended Eulerian time-stepping scheme for the heat equation on time-dependent domains.Theoretically, stability and optimalorder convergence estimates were derived in the energy norm under the assumption of a parabolic CFL condition.In the numerical results, on the other hand, we did not observe any stability issue related to a violated CFL condition.The three-dimensional results for second-order polynomials indicate that a violated CFL condition could result in a slightly reduced convergence order in the L 2 (H 1 av )-norm.To our knowledge this is the first work, in which an implicitly extended Eulerian time-stepping scheme is applied with a scheme that requires derivative information at different time steps.As mentioned in the introduction, this could be the basis for an analysis of a whole zoo of time-stepping schemes, such as the Fractional-step-θ method, implicit Runge-Kutta-or Adams-Bashforth schemes.Moreover, we plan to apply the developed time-stepping scheme to flow problems on time-dependent domains and to fluid-structure interactions with large displacements, see e.g.[27,14].

Figure 2 :
Figure 2: Left: Illustration of the triangulations T n h and T n h,δ and the sets of facetsF n h,δ = F n,int h,δ ∪ F n,cut and let the CFL condition (Assumption 2) be valid.It holds for l ∈ {n − 1, n}

Theorem 4 . 5 (
Global error) Based on the assumption made in Lemma 4.4, the global error e

Figure 4 :
Figure 4: Illustration of the numerical solution of Example 5.2 at time t = 0 (left) and t = 1 (right).
the first term in (4.13) vanishes.By definition of f in (4.7) we can estimate further ∆t

Table 1 :
L 2 (T ), L 2 (L 2 ) and L 2 (H 1 av )-norm errors for P 1 finite elements applied to Example 5.1.The estimated orders of convergence are computed according to (5.1)-(5.3).The diagonal orders are computed from the underlined error values.

Table 2 :
1 av )-norm errors are shown in Table2.Firstly, we observe that the absolute values of the errors are significantly smaller compared to L 2 (T ), L 2 (L 2 ) and L 2 (H 1 av )-norm errors for P 2 finite elements applied to Example 5.1.