{'fraction_non_alphanumeric': 0.04544747888272788, 'fraction_numerical': 0.011020711336822649, 'mean_word_length': 4.13500088699663, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 8, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Geometrically Higher Order Unfitted Space-Time Methods for PDEs on Moving Domains

In this paper, we propose new geometrically unfitted space-time Finite Element methods for partial differential equations posed on moving domains of higher order accuracy in space and time. As a model problem, the convection-diffusion problem on a moving domain is studied. For geometrically higher order accuracy, we apply a parametric mapping on a background space-time tensor-product mesh. Concerning discretisation in time, we consider discontinuous Galerkin, as well as related continuous (Petrov-)Galerkin and Galerkin collocation methods. For stabilisation with respect to bad cut configurations and as an extension mechanism that is required for the latter two schemes, a ghost penalty stabilisation is employed. The article puts an emphasis on the techniques that allow to achieve a robust but higher order geometry handling for smooth domains. We investigate the computational properties of the respective methods in a series of numerical experiments. These include studies in different dimensions for different polynomial degrees in space and time, validating the higher order accuracy in both variables.

1. Introduction. Many problems in physics, engineering, chemistry and biology can be described in terms of Partial Differential Equations (PDEs) posed on moving domains. Examples include multi-phase flows [17], blood flow in the human heart [31] or evolving biological cells. Finite Element methods (FEM) provide a powerful framework to solve these problems numerically. Traditional fitted mesh approaches, where the computational mesh also prescribes a parametrisation of the geometry, requirewhen applied to moving domain problems -mesh adaptations or remeshings in every time step. Such approaches are typically employed in Arbitrary Lagrangian-Eulerian (ALE) methods [21], and are computationally attractive for problems involving small deformations or without topology changes, but become difficult to realize when strong deformations or topology changes occur.
In recent years, geometrically unfitted Finite Element methods have been introduced and studied for a variety of problems under different names such as CutFEM [10], Finite Cell Method [33], XFEM [16], fictitious domain methods [12,13] and TraceFEM [32]. In these methods, the computational mesh is not aligned with the geometry of the problem, avoiding the necessity for a mesh generation or remeshing procedure. However, unfitted FEM have to deal with different issues such as robustness w.r.t. the position of the geometry relative to the computational mesh, the imposition of interface and boundary conditions, accurate geometry handling and time integration on moving domains. In this manuscript we present a class of unfitted Finite Element methods aimed at addressing the last two aspects: Robust and higher order 1 space-time integration for moving domain problems.
1.1. Approaches in the literature. Most fitted FEM discretisations use a version of the Method-of-Lines approach to deal with the time derivative exploiting that every unknown has a well-defined history as it is associated to an entity of the mesh. For unfitted FEM, due to the decoupling of the computational mesh and the geometry, the situation is entirely different. A naive Method-of-Lines approach may not even be well-defined when unknowns associated to entities of the computational mesh enter or leave the time-dependent domain, cf. Figure 1 for a sketch.
In characteristic-based methods also known as semi-Lagrangian methods such as those discussed in [30,5] instead of the partial derivative ∂ ∂t the material derivative d dt = ∂ ∂t + w · ∇ is approximated. The material derivative is the rate of change along the trajectories of the flow field w that also advects the geometry. Provided that these trajectories can be resolved numerically in an efficient and accurate way it can be ensured that a finite difference approximation, as in the Method-of-Lines approach, is a reasonable approximation of the material derivative.
Extension-based methods use a finite difference approximation of the partial time derivative, but apply an additional extension to the solution within every time step. The extension makes sure that the discrete solution is well-defined in a small neighbourhood of the domain, providing a reasonable history for each unknown involved in the finite difference approximation at the next time step. Such a method has been introduced in [25] for a scalar convection-diffusion equation. Stokes and Navier-Stokes problems with this approach have been considered in [35,11,37].
Space-time methods do not apply finite difference-type approximations to the time derivative, but opt for a variational time discretisation which together with an FEbased space discretisation results in a variational formulation in a space-time domain. Often (but not always) the space-time domain is decomposed into time slabs that correspond to small time intervals and a tensor-product discretisation (w.r.t. meshes and FE spaces) is chosen. This allows to keep the structure of a time stepping scheme and to build on components of related spatial problems (basis functions, quadrature rules, etc.). The former aspect implies that the PDE solution on the space-time domain can still be solved for with a series of problems with an essential complexity of a spatial PDE problem.
In this manuscript we restrict to the latter class of space-time methods in the context of geometries that are described by level set functions. Among this class of methods, [26,22] consider Discontinuous Galerkin approaches in time of second order accuracy for two-phase flows. Moreover, in [18,38] space-time methods for moving surface/ coupled surface-bulk problems are presented for second order accuracy as far as the level set description is concerned, and third order accuracy in the context of a geometry description involving splines. In [29] unfitted space-time methods based on a geometry description with overlapping meshes with piecewise linear-in-time motion are considered. Recently, in [1] discontinuous-in-time methods with higher order accuracy have been investigated numerically for a Navier-Stokes problem.

Main contributions.
In this manuscript, we aim to address three inherent challenges of higher order unfitted space-time methods: 1. The derivation of discretisations that are robust, i.e. stable and accurate inde-pendent of the essentially arbitrary cut configurations that can occur in space-time, is a crucial component to obtain a reasonable method. 2. One inherent challenge of unfitted methods where the geometry is described by level set functions -especially for the space-time setting -is that of finding proper means to handle the implicit space-time geometry robustly and accurately for the potentially complicated domains of numerical integration.
3. Another challenge is the computational complexity. Although we can maintain a time stepping structure, in the higher order (in time) case we will still have considerably more unknowns to solve for in each time step than for a standard time stepping scheme. This renders space-time methods typically much more expensive compared to other time stepping schemes.
To obtain stable variational formulations, we will build on the ideas of the previous works [26,18] in the context of DG-in-time methods of lower order. One major contribution in this work is the extension of these ideas to higher order in space and time and to Petrov-Galerkin methods with discrete solutions that are (higher order) continuous in time. The transition from discontinuous solutions to continuous solutions or even solutions with higher regularity in time reduces the number of unknowns that are involved in each time step which compensates for the higher complexity of space-time methods. The second major contribution in this work is the extension to the space-time setting of a higher order accurate geometry handling based on isoparametric mappings for level set domains as introduced in [23]. Thereby we can ensure to obtain arbitrary high order accuracy in space and time, including considerations of geometry approximation. Parts of the ideas in this manuscript have also been discussed in the Master's theses [34,19].
1.3. Structure of this paper. The remainder of this paper is organised as follows: In section 2, the convection-diffusion problem is introduced as a model problem. Furthermore, the different time discretisations are illustrated by semi-discrete (discrete in time) variants of the space-time methods coming up in the subsequent sections. Section 3 is devoted to the higher order approximation of moving space-time domains described by level sets. Afterwards, in section 4, the fully discrete methods are defined in their discrete variational form. In section 5, a specific strategy for higher order integration in time is discussed. Then, section 6 presents a series of numerical experiments, which demonstrate the higher order accuracy of each method, and investigate aspects of accuracy and computational complexity. Finally, section 7 concludes the paper with a summary and an outlook on open problems.

Model problem and time discretisation.
2.1. The model problem. Before we come to the specific methods, we introduce the model problem used in this study. The model problem is given as a convection-diffusion equation of a species of concentration u on a moving domain Ω(t). We denote the corresponding space-time domain as Q := t Ω(t) × {t}. As initial condition, we assume a function u 0 : Ω(0) → R to be given, and for the boundary condition we assume that no transport of the species across the domain boundary occurs. We arrive at the following problem: Find u : Q → R, s.t.
with initial data u(·, 0) = u 0 in Ω(0). Here, w represents a divergence-free convection field with w L ∞ (Q) ≤ c < ∞ and f sources or sinks for the species. We further assume that the convection field w coincides with the domain motion on the spatial boundary, so that on Γ * := t∈(0,T ] ∂Ω(t) × {t} there is (w, 1) ⊥ n * where n * is the space-time normal to Γ * . In order to formulate a well-posed variational formulation in a space-time setting, we define the Sobolev space of functions with weak derivatives in spatial direction: and denote its dual as H −1,0 (Q) := (H 1,0 (Q)) * . The weak formulation then reads: where ·, · denotes the duality pairing between H 1,0 (Q) and H −1,0 (Q) and (·, ·) G denotes the L 2 inner product on a domain G. This weak formulation is well-posed, cf. also [17,Section 10.3] for the more involved case of a moving interface problem.
Remark 2.1. Let us stress that although this model is comparably simple, as it is a scalar parabolic equation with very simple boundary conditions and the domain motion is assumed to be known, the techniques that are to be discussed in the remainder of this manuscript can also be applied to more complex problems, e.g. vectorial free boundary problems with possibly different boundary conditions. This especially holds for the geometry handling, which is a main feature of this study.
Next, we introduce three different suitable time-discrete variational formulations of (2.3). With these, several features of moving domain problems in an Eulerian setting and the considered time discretisation approaches become visible. In this sense this section also serves as a preparation for the fully discrete formulations in section 4.

Discontinuous Galerkin in time formulation.
We assume that the time interval [0, T ] is subdivided into time intervals, I n = (t n−1 , t n ], n = 1, . . . , N where 0 = t 0 < t 1 < · · · < t n−1 < t N = T , and, for ease of presentation, we further assume that all slices have the same length, ∆t = t n − t n−1 = T /N . Let Q n := t∈In Ω(t) × {t} be the space-time domain corresponding to one time interval I n . Further, set Ω n = t∈In Ω(t) and let Q n = Ω n × I n be the smallest tensorproduct domain containing Q n , see Figure 2 for an illustration. We define the ansatz and trial space so that on each time interval the functions are of tensor-product form and polynomial in time with k t ∈ N 0 Here and in the following, we denote by P k (E) the set of polynomials on the geometric entity E of order smaller or equal k. The space-time formulation is obtained from Fig. 2. Schematic illustration of (tensor-product) space-time domains on the time slabs. reducing (2.3) to one time slab and adding an upwind stabilisation in time. For one time interval n ∈ {1, .., N } it takes the form: Find u ∈ W n , s.t. for all v ∈ W n there holds where u n−1 + (·) = lim t t n−1 u(·, t), and similarly for v. Here, u n−1 − takes the role of the initial data for the current time slab n. For n = 1 it is simply u 0 , the given initial data, for n > 1 it is lim t t n−1 w(·, t) for w ∈ W n−1 the discrete solution from the previous time step. Subsequent solution of (2.5) for n = 1, .., N yields the global (discontinuous-in-time) space-time solution.
Remark 2.2. With similar techniques as in [26] 2 it is fairly easy to show unique solvability and some basic error estimates in the H 1,0 norm. We are however not aware of any analysis providing superconvergence results also for the L 2 (Ω(T ))-norm as known from the stationary domain case (see e.g. [36]).

A Continuous Galerkin in time formulation.
From the case of a stationary domain Ω it is known that a Continuous Galerkin (CG) (in time) trial space can be used while keeping a time stepping structure if a discontinuous (in time) test space (of one degree less) is used, see e.g. [4]. In this case within a time slab one dof (in time) corresponds to an H 1 (Ω) function known from the previous time slab. Hence, only k t (instead of k t + 1) unknown H 1 (Ω) functions remain.
Difficulties arise if one wants to apply this strategy to the case of a moving domain. To illustrate this, let us consider a time slab Q n , n > 1. To solve on Q n we want to use the tensor-product space H 1 (Ω n ) ⊗ P kt (I n ) again. However, in the previous time slab we only solved for a solution in H 1 (Ω n−1 ) ⊗ P kt (I n−1 ) so that the initial data for the current time slab n is only defined on Ω n−1 ⊃ Ω(t n−1 ), but we can have Ω n Ω n−1 , i.e. the initial data is insufficient to remove the first dof (in time). The geometrical situation is sketched in Figure 2. To overcome this issue we combine the solution step on every time slab n with a Sobolev extension E of the solution u(t n ) ∈ H 1 (Ω n−1 ) on H 1 (Ω n + ) with Ω n + s.t. Ω n , Ω n−1 ⊂ Ω n + . Thereby we can match up the degrees of freedom (dofs) (in time) for trial and test space again, which (on each time slab n) we choose as where u n−1 − ∈ H 1 (Ω n−1 ) is the solution from the previous time slab (or given initial data), φ 0 is a function in P kt (I n ) with φ 0 (t n−1 ) = 1 and P kt 0 (I n ) is the space of polynomials up to degree k t with value zero at t n−1 . The space-time formulation on time slab n ∈ {1, .., N } then takes the form: Find u ∈ U n , s.t.
Again, subsequent solution for n = 1, .., N yields the global space-time solution.
2.4. Galerkin Collocation formulations with higher regularity. In the recent decade space-time finite element methods of Petrov-Galerkin-type with higher regularity became popular [7,3,6,2]. The idea is to impose a higher regularity m r ≥ 1 in time on the trial space which allows to reduce the number of temporal unknowns per time slab. Straightforward generalizations of the previously introduced scheme are however not necessarily stable which leads to the following modification: Instead of using test functions of degree k t − m r − 1 to match the dofs (in time), a smaller degree m t is used and the remaining dofs are used to impose the PDE (in weak form) at m c fixed time instances (the collocation points) yielding a Galerkin-Collocation (GCC) method.
With U n 0 := H 1 (Ω n ) ⊗ P kt (mr),0 (I n ) we choose the trial and test spaces where u n (l),− ∈ H 1 (Ω n−1 ), l = 0, .., m r are given from the previous time slab (or given initial data), φ l is a function in P kt (I n ) with φ where t n l ∈ I n , l = 1, .., m c are the collocation points. Note that for m r = m c = 0 we recover the CG method from Subsection 2.3.
3. Geometry handling. Before we consider temporally and spatially discrete versions of the formulations in the previous section, we discuss the discrete geometry handling in the case of an unfitted geometry. In this setting the geometry description is separated from the computational mesh. More specifically, we assume that we are given a level set function on a background domainΩ. The level set function φ :Ω × [0, T ] → R implicitly describes the space-time geometry by Due to the implicit nature of the geometry description it has no explicit parametrisation. It is hence not immediately clear how to realize numerical integration on Q or Q n or suitable approximations. In this section we introduce an approach that allows for higher order accuracy in space-time while still guaranteeing positive quadrature weights. This approach relies on the idea of isoparametric unfitted FEM that has been established for the stationary case in [23]. After a brief repetition of these techniques we extend the approach to the space-time setting exploiting tensor-product structure.
3.1. Isoparametric mappings for stationary level set domains. In the following, we summarise the method introduced in [23] which prepares the techniques and notations used later on in the space-time setting. We assume a domain Ω with smooth boundary to be given as well as a background domainΩ ⊃ Ω with a shape regular triangulation T h . The domain Ω is described by a (within the scope of this subsection only spatial) level set function φ : Ω → R, such that Ω = {x ∈Ω | φ(x) < 0}. The geometry description in terms of φ is only implicit and hence does not offer a parametrisation that can be directly used for setting up numerical integration routines. To enable the use of robust numerical integration routines, we define a piecewise linear approximation of φ in two steps: Fig. 3. Sketch of a simple subdivision that allows for setting up quadrature on T ∩ Ω lin (left), the discrete geometry Ω lin (middle) and the mapped geometry Ω h obtained from the mapping Θ h (right); stationary case.
(1) We introduce a polynomial approximation of higher order This will lead to an approxi- As with φ, this description of the approximated domain is still implicit so that it cannot be exploited for numerical integration directly. However, in the following we will construct an approximation which is as good as Ω h and has an explicit parametrisation.
(2) Next, we define by φ lin ∈ V 1 h a continuous elementwise (multi-)linear approximation of φ h leading to the geometry approximation Ω lin = {x ∈Ω | φ lin (x) < 0}. On simplex meshes Ω lin is polygonal and each cut element can easily be decomposed into uncut elements on which numerical integration can easily be applied, cf. the sketch on the left of Figure 3. Also on non-simplex meshes setting up numerical integration is simplified for Ω lin compared to Ω (or Ω h ), cf. [20].
This simplification in the geometry handling comes at the price of accuracy. By construction Ω lin is only a second order approximation to Ω. To improve the accuracy of the approximation we introduce a vectorial finite element function that serves as a mesh transformation Θ h : The transformation Θ h is constructed and applied on T h such that Ω h := Θ h (Ω lin ) is a higher order approximation to Ω, cf. Figure 3 for a sketch. Note that for technical reasons discussed in [23], Ω h is not exactly equal to Ω h stemming from φ h . However, their difference is asymptotically small, dist(∂Ω h , ∂Ω h ) h qs+1 , i.e. small enough to yield dist(∂Ω h , ∂Ω) h qs+1 .
Problems of numerical integration on Ω h can then be reformulated as numerical integration problems on Ω lin with transformation factors involving Θ h . Corresponding to the deformation of the mesh the involved finite element spaces are also mapped as usual in isoparametric finite element methods.
For the construction, initially a preliminary mapping is set up on a subset of elements in T h which is typically chosen as the domain of cut elements In a blending step, a transition to the identity mapping is applied on all neighbours to cut elements. On all remaining elements the mapping is set to the identity. For the construction of the mapping we introduce the map where we refer to the power set of T h by Pow(T h ). We refer the reader to [23] and [27] for a detailed account on the construction and analysis of the mapping for the stationary case.
3.2. Isoparametric mappings for space-time level set domains. We want to transfer the previous ideas to the space-time setting exploiting a tensor-product structure. We consider a time interval I n = (t n−1 , t n ] and define the space-time tensorproduct mesh of prisms by extrusion {T × I n | T ∈ T h }. In this setting, we assume that our space-time domain Q n is described by a level set function φ :Ω × I n → R which is approximated by a space-time level set function in a tensor-product finite element space, φ h ∈ V qs h ⊗ P qt , where q t is the order of the temporal accuracy. In the following we also use q = (q s , q t ) to summarise the orders of spatial and temporal geometrical accuracy and write q = r if q s = q t = r. Additionally we make the following assumption.
Assumption 3.1. The approximated level set function φ h :Ω × I n → R is continuous across all time slabs and hence ensures a connected approximated space-time geometry, cf. Figure 6 for a sketch.
We are now going to construct a space-time mapping Θ n h :Ω × I n →Ω based on the tensor-product structure of φ h and hence look for a space-time transformation Similarly to the stationary case discussed in the previous section, we will use a level set function that is elementwise (multi-)linear in space to set up a reference configuration. Note that we keep the order in time fixed, i.e. higher order. Let I 1 h : C(Ω) → V 1 h be the nodal (spatial) interpolation operator. We define We develop the transformation Θ n h in the same tensor-product structure, and construct each spatial transformation Θ h,i (x) based on φ i h (x) and I 1 h φ i h (x) as in the stationary case using the map C. It is crucial to make sure that the set of active elements is the same within one time slab, cf. Remark 3.2 below. We hence define The overall procedure is applied to each time slice I n , leading to a mapping Θ n h :Ω × I n →Ω for each I n respectively. A schematic illustration is provided in Figure 4. Moreover, by Θ n,st  denote the function mapping each ( The constructed mapping will in general not be continuous in time due to different active elements on different time slabs. However, Ω lin (t) and Ω h (t) are continuous.
In numerical experiments (cf. Figure 5, section 6), we confirm that Q h,n satisfies the following approximation error bound Remark 3.2. A subtle, but very important detail is the decision on the third argument of C, the subset of actively deformed elements. A seemingly straightforward approach would associate to every i ∈ {0, . . . , q t } a set of active elements T i h according to the cut elements associated to I 1 h φ i h . This would easily allow to make sure that the constructed transformation is also continuous in time 4 . However, we observed in numerical experiments that this method would not yield the optimal approximation quality (3.8) for higher orders. In the remainder we consider the previously presented construction with a discontinuous-in-time mapping. The question if a similar construction yielding a continuous deformation can be found which preserves the desired accuracy is left open for future research.
A proper treatment of the discontinuity in the mesh deformation is discussed next.

Discontinuous-in-time deformations and mesh transfer operations.
With the previously presented construction of the deformation, for successive time intervals I n = (t n−1 , t n ], I n+1 = (t n , t n+1 ], we will obtain deformations Θ n h (x, t) and Θ n+1 h (x, t), respectively, where in general. One specific T ∈ T h can be associated to different regions of the isoparametric mapping construction, e.g. T ∈ T Γ h,n+1 but T ∈ T Γ h,n , cf. Figure 6 (left) for an illustration of that situation.
The challenge within this setting stems from the parametrically mapped discrete function spaces. To take the deformation into account, a discrete space-time function defined on I n , restricted to t n , will have the structure for a discrete representation of the initial data for I n+1 for the CG and GCC method. A standard L 2 projection from one deformed mesh to the other would be a valid option. Within the context of Eulerian time stepping methods for higher order methods in space and time, a similar transfer operation is needed. In [28, Section 5] a transfer operation between differently deformed meshes is presented (and thoroughly analysed). It is explicit and requires only essentially element-local operations leading to an operation that is more efficient to evaluate than a standard L 2 projection and has negligible costs compared to the remaining operations required to solve the discretisations discussed in the remainder. We can fortunately exploit this operation also in our (in comparison simpler 5 ) setting as well which we especially did in the numerical examples. For further details of the construction and analysis of the transfer operator we refer to [28].
Concerning the DG method, there is no need for a discrete projection of u − ∈ V qs h • (Θ − h (t n )) −1 , as the initial data is only imposed weakly. It is however computationally convenient to have a discrete approximation living on the mesh deformed by Θ + h (t n ). This does not need to be a continuous discrete function, so that we recommend to apply the previously discussed projection only element-wise and skip the final averaging step, called P h in [28, Subsection 5.1.3]. The so-constructed transferred discrete projection is denoted as Π * u (u − ) from now on, so that Π u (u − ) = P h Π * u (u − ). Remark 3.3 (Accumulation of projection errors from the transfer operation). In general the application of a transfer operator in each time step of a time stepping schemes can lead to unbounded accumulated projection errors for increasing number of time steps, cf. for instance [8]. However, one important insight in [28,Section 3.4] is that the special situation of only very slightly different mesh deformations considered here leads to the following important property: For a fixed element the number of times that the change of the mesh transformation is not Lipschitz-continuous in time is bounded independent of the number of time steps. Furthermore, the projection error accumulation of transfer operations corresponding to Lipschitz-continuous changes in the mesh deformation as well as the projection error accumulation due to the small number of discontinuous changes in the mesh deformation stay bounded independent of the time step size. Hence, overall the accumulation of projection errors stays bounded. 4. Fully discrete weak forms. In this section, we introduce the discrete weak forms of the three methods discussed in the semi-discrete setting in section 2. The fully discrete version of the DG method is presented in Subsection 4.1, the CG method is introduced in Subsection 4.2 and the GCC method in Subsection 4.3. In the following when we distinguish the methods as DG, CG or GCC we refer to the discretisation in time only. In space we always only consider the same continuous finite element space.

Discontinuous Galerkin method.
To state the discrete variational formulation of the DG method, we have to introduce several discrete regions. Some of them rely on the spatially linear interpolation of the level set function φ lin . First, we define an extended space-time domain, which contains all points in Q lin,n and has tensor-product structure per time slice: In addition, we define an interior space-time domain, which is the largest domain with time-slice tensor-product structure which is completely contained in Q lin,n : I(Ω lin,n ) : = {x ∈ T for some T ∈ T h such that (T × I n ) ⊆ Q lin,n }.
For both spatial regions, we define a space-time counterpart: E(Q lin,n ) := E(Ω lin,n ) × I n , I(Q lin,n ) = I(Ω lin,n ) × I n .
As usual in unfitted methods, we use a ghost penalty stabilisation to handle ill-posed cut configurations. It is defined almost on the facets between interior and exterior, For the case when E(Ω lin,n ) \ I(Ω lin,n ) has a width that is significantly larger than h we include a slightly larger set of nearby interior facets, F n,ext R ⊇ F n R , so that the number of facets in the interior is proportional to that in the exterior. We refer to [19, Assumption 3.1] for a precise explanation. These regions are illustrated for a setting of a one-dimensional domain in Figure 7. Now, we are in a position to introduce the discrete function spaces. To this end, first a finite element space for each time slab is defined as where k = (k s , k t ) is the tuple describing the combination of the polynomial orders in space and time and again for k s = k t = r we write k = r. Now we introduce a general notation for cut discrete spaces depending on a geometrical domain E ⊆Ω × [0, T ]: This means that only those dofs that are associated to E stay active whereas all others are set to zero. Note that with this definition of the space only its restriction to E is uniquely defined 6 , which however suffices as functions in W n,n h,cut (E) will not be involved outside of E. This allows to define the finite element space for the DG method as follows: Using the discrete regions and function spaces introduced above, we are now able to define the first components of the DG method for one time slab: where u n ± := lim s tn u(·, s) and u n−1 − is the initial value for the current time slab. For n = 1 this is the initial value of the global problem whereas for n > 1 this stems from the solution of the previous time step. The transfer operation Π * u is applied here as Θ st h may be discontinuous across t n−1 .
In addition, we define a ghost penalty stabilisation. We use a variant which exploits facet patches, which we call direct version. 7 Let F = T 1 ∩ T 2 be a facet of the undeformed triangulation, T 1 , T 2 ∈ T h . Then, the facet patch is defined as ω F := T 1 ∪ T 2 . Taking into account the mesh deformation, we define also ω h F (t) = Θ h (ω F , t). In addition, we need a volumetric jump operation on these curved elements. Let us assume u =û • (Θ st h ) −1 , and we are interested in the jump for some point x in the mapped element T 1 . Then, , where the latter term extends the (mapped) polynomial on T 2 to T 1 based on the canonical extension of polynomials from their element T ∈ T h to R d (denoted by E p ) combined with the involved mesh transformations. Then, defines the ghost penalty bilinear form which -roughly speaking -allows to carry over control from one element to its neighbour. Here,γ J = 1 + ∆t h γ J , where we refer to [25,Section 4.4] for an explanation of the scaling with ∆t and h that is relevant in the anisotropic case h ∼ ∆t. Note that the analysis guarantees that the chosen scaling is sufficient for stability and bounded condition-numbers. The question if this scaling is also necessary is left open for future research. γ J is a stabilisation constant to be chosen sufficiently large, c.f. [34,19] for an analysis and Subsection 6.1.3 for a computational study. In total, the discrete time slab problem reads: Find u ∈ W n h such that The time slab solutions can be computed one after another yielding a global spacetime solution u h . The same applies to all following discrete problem statements. 6 outside of E the decay to zero depends on the chosen basis for W n,k h . 7 Other versions appear in the literature, such as the original version [9] or the normal-derivative jump version [13]. See [34,Remark 6] for a comparison in terms of analysis.
Ω lin (t) Fig. 8. Illustration of the domains of definition of the initial data at tn as needed and given from the previous time step if the DG-regions were used. The region where an extension of u from In is needed is depicted in purple.

Continuous Galerkin method.
In general, the structure of the discrete regions for the CG method resembles that of the DG method. For the dof associated to the end point of a time slice, we introduce E + (Ω lin,n ), an extension of E(Ω lin,n ), E(Ω lin,n ) ⊆ E + (Ω lin,n ), which additionally should satisfy the constraint This is the discrete version of the domain Ω n + introduced in Subsection 2.3, with the only difference that the discrete extended domain aligns with the mesh, as illustrated in Figure 8. We start with a variant of Ω lin (t), Ω lin (t) := {x ∈Ω | φ lin (x, t) < }, which is used to define a strip around the discrete boundary S n δ , and E + (Ω lin,n ): Note that for a time slice I n , the regions E + (Ω lin,n ) and S(Ω lin,n ) depend only on the discrete function φ lin at t n , the end time point of the interval. Moreover, δ is assumed to be given as follows: Assumption 4.1. We assume δ to be some small number such that E(Ω lin,n ) ⊆ E + (Ω lin,n ) and the constraint (4.6) is satisfied.
If φ lin would be a signed distance function, we could set δ ≥ ∆t w L ∞ (0,T,L ∞ (Ω)) . For practical computations, we recommend to multiply this expression by a factor f > 1 and to check (4.6) in each time step to optimize f when necessary. For a sufficiently large f , our framework reduces to a global extension on the whole mesh, such as suggested in [1]. We comment on this option in Subsection 6.2.2.
The discrete extension should be realised by a merely spatial ghost penalty stabilisation. The according facet set is hence given as F n,+ R = {F ∈ F s.t. ∃T 1 = T 2 , T 1 ⊆ E + (Ω lin,n ), T 2 ⊆ S(Ω lin,n ) with F = T 1 ∩ T 2 }. Again, we give a sketch of this construction in Figure 9. In order to introduce the discrete function spaces next, we pose an assumption on the basis of the function space P kt ([t n−1 , t n ]) contained in W n,k h : Assumption 4.2. The basis {p 0 , . . . , p kt } of the space P kt ([t n−1 , t n ]) contains two functions p i lo (t) and p iup (t) which can be associated to the lower and upper domain boundary t n−1 and t n , s.t. p i (t n−1 ) = δ i lo i and p i (t n ) = δ iupi ∀i = 0, . . . , k t .
The discrete trial function space is then defined as follows in terms of the nomenclature known from the DG method: where Π u (u n−1 − ) serves as initial data for the solution as for the DG method, but with a strong imposition. The relevant discrete region for each time slice is given as E(Q lin,n )∪(E + (Ω lin,n )×{t n }), where the expression E + (Ω lin,n ) × {t n } will allow for the discrete extension motivated above at the end of each time slice. The dofs from W k h related to the region E + (Ω lin,n ) × {t n } will stem from p iup , so the construction is well-posed by Assumption 4.2. The structure of the domain is illustrated in Figure 9 by black circles. For the space of test functions, continuity beyond time slice boundaries is not required, so we define: Note that the polynomial order in time is one less as usual with CG methods. Now we have all tools in place to state the variational formulation of the CG method within one time slab. It reads: Here, both ghost penalty parts only act on the dofs of u 0 that needs to be solved for, but not on u init . We further note that the decomposition u = u 0 + u init and hence the solution u depends on definition of the basis function p i lo .
Remark 4.3. There exists a variant of the suggested CG method which we denote as CG , where only the ghost penalty stabilisation J n (u, v) = tn tn−1 j h (F n,+ R , t; u, v) dt is applied, which is defined on the smallest tensor-product domain which contains both domains used before. This can be regarded as an implementational simplification, although more dofs are needed as elaborated on in Subsection 6.2.2. Further details about this variant are given in the supplementary material.

Galerkin-Collocation method.
The GCC method builds on the same idea as the CG method, but extends the strategy by one step: In addition to the continuity of the discrete solution function, we ask for continuity also in higher order time derivatives up to order m r . Correspondingly, we strengthen the assumption of Assumption 4.2 w.r.t. to the basis of P kt (I n ): and u init (x, t) = mr =0 Π u (∂ t u n−1 − (x)) · p i (t) where Π u (∂ t u n−1 − ) serve as initial data for the solution as for the CG and the DG method. We note however that we have to prescribe m r + 1 functions as initial data here. For the space of test functions, we now pick only order m t ≤ k t − m r − 1: For the ghost penalty we use the tensor-product variant CG from Remark 4.3 8 . The formulation of the GCC method then reads: Find u = u 0 + u init with u 0 ∈ U n h,0 s.t.
with collocation points {t n 1 , . . . , t n mc }. Here, similar to the CG method, we made sure that the ghost penalty only acts on the dofs of u 0 . Again, the solution u will depend on the concrete choice of the basis functions p i , = 0, .., m r .
A possible choice for the triple (k t , m r , m c ) is (3, 1, 1) which we also consider in the numerical experiments below. The one collocation point in this setting is chosen as the end point of each time interval, t n 1 = t n mc = t n . 5. Implementation of numerical integration in time. The discrete (bi)linear forms presented in the previous section involve summands defined on the spacetime domains Q h,n , c.f. e.g. (4.4). Hence, a numerical integration procedure is needed for those domains. More specifically, as Q h,n was defined as the image of Q lin,n , cf. (3.7), and finite element computations are done element-wise, we apply the following calculation: For the other summands, similar equations hold. Hence, we only need to implement an algorithm to approximate for each T ∈ T h for some integrandf . For a given timet, an integration rule on Ω lin (t ) ∩ T is assumed to be available, e.g. on simplices by tessellation. Hence, the inner integral in (5.3) is well-suited for numerical approximation, and it remains to replace the integration tn tn−1 by a numerical quadrature. In the following, we will define two variants of doing so, where the first approach is mentioned for comparison with existing literature [18,38] and the second is our suggested approach.
Definition 5.1 (Topology-insensitive time integration). In the topology-insensitive time integration, tn tn−1 is approximated by Gaussian quadrature. Let (t i , ω i ) i be the Gaussian quadrature rule on [t n−1 , t n ] of order k t,int . Then, (5.3) will be approximated by Note however, that the topology of Ω lin (t) ∩ T might change over t ∈ I n ; for instance the domain boundary can leave or enter an element. This will cause the integrand Ω lin (t)∩T f (x, t) dx to show a discontinuity in the first derivative even if f is smooth, which impedes the accuracy of the topology-insensitive time integration. To circumvent this, we define a time integration that takes the cut topology changes of Q lin,n into account. Define the scalar function φ v : I n → R, t → φ lin (v, t).

5:
Search for the roots R v of φ v .

6:
Set R ← R ∪ R v . 7: end for 8: Define R * as the set of intervals with endpoints according to R. 9: return R * .
Denote by I s n = (t l , t u ) ∈ R * a subinterval in R * . Then, (5.3) is approximated by the topology-preserving time integration as where (t s i , ω s i ) is the Gaussian integration rule transformed to the interval I s n . This time integration procedure exploits the fact that a (multi-)linear-in-space level set function is fully described by its function values on vertices of the finite element T . The functions φ v within Algorithm 5.1 will in general be of higher polynomial order.
The topology-preserving time integration is our method of choice and for the numerical studies in this work we choose integration rules in time of exactness degree 2(k t + 1). In Figure 10, we show the quadrature points generated by both approaches for simple examples in one and two space dimensions. Later in this paper, we investigate numerically differences to the topology-insensitive time integration in Subsection 6.2.3.

Numerical investigations.
In this section, we want to demonstrate the numerical properties of the proposed methods. We consider three problem geometries to investigate specific aspects, each in one of the following subsections. Recall that the notation k = r means k s = k t = r and analogously q = r means q s = q t = r.
All experiments are performed with ngsxfem [24], an unfitted finite element extension of ngsolve. Reproduction data are available at https://gitlab.gwdg.de/fabian. heimann/repro-ho-unf-space-time-fem.  The linear systems are solved with the direct solvers umfpack and pardiso of IntelMKL. As a basis for P kt ([t n−1 , t n ]), we implement a Lagrange basis with the Gauß-Lobatto points between t n−1 and t n of order k t as Lagrange points for DG and CG. These satisfy Assumption 4.2. For GCC of order 3, we implement cubic Hermite spline polynomials, which satisfy Assumption 4.4.
6.1.1. Geometry approximation quality. In a first experiment, we investigate the quality of the discrete domain approximation by measuring the distance to the exact interface, max t∈[0,T ] dist(∂Ω h (t), ∂Ω(t)), for the kite geometry with k = q and refinements in space and time, i = i t = i s . As can be seen in Figure 5b, the discrete regions satisfy the asymptotic error bound given in (3.8).
6.1.2. Convergence of methods of matching space and time order. First, we evaluate convergence of the numerical error for i s = i t = i, k = q. The range of parameters k ∈ {1, 2, . . . , 6} is applied to the DG and CG methods, whilst the GCC method is investigated for k = 3. In Figure 11 (left) and Figure 12, we display the numerical errors for simultaneous space-time refinements i = 1, 2, . . . . We measure the numerical error in terms of the following norms: u(x, t) 2 dx dt.
We find both errors to decrease as follows, where k = k s = k t = q s = q t , h ∼ ∆t,   This confirms the expected higher order convergence property of the suggested methods.
6.1.3. Choice of the stabilisation parameter. The stabilisation parameter of the ghost penalty has to be chosen large enough to fulfill sufficient conditions for stability of the method and to obtain reasonable bounds on the condition numbers of arising linear systems [10,35,34]. In Figure 11 (right) we illustrate the dependence of the numerical error on γ J for the problem considered in this subsection using i = i s = i t and k = q = 4 for the DG method. Overall, γ J influences the absolute values of the errors, whilst the asymptotic behaviour remains intact even if the parameter is chosen relatively large. A particular limit is that of stabilisation constant γ J → 0. There, we observe that the method is stable in general. However, condition numbers may blow up, which impedes the practical solution of the linear system at some refinement level/ polynomial degree. (In this case, i = 2 is the last level solved properly.) We explain this behavior as follows: For mere stability, i.e. assuming exact arithmetics, ghost penalty stabilization is not even necessary. Note that in [26] even robust, i.e. cutposition independent, a priori error bounds have been derived without ghost penalty stabilization in a comparable setting (involving the assumption of exact integration). However, taking computer arithmetics, i.e. round-off errors and their impact, into account the possible ill-conditioning of linear systems is practically hardly acceptable. A more detailed investigation of these observations (potentially also in combination with different time integration strategies) in terms of a full numerical analysis would be an interesting task for future research.
(Ω(t)),0,T ) Fig. 13. Superconvergence study of the DG and CG methods for kt = 1, 2 where ks = 2kt + 1. Refinement in time for a fine mesh. On the left-hand side, the error is measured in the L 2 (Ω(T ))norm, whilst on the right-hand side the same results are investigated in the L 2 (L 2 (Ω(t)), 0, T )-norm.
In conclusion, we may pick γ J small as long as a sufficiently accurate (and efficient) solution of linear systems is possible. We choose γ J = 0.05 for the remaining studies.
6.1.4. Superconvergence. The investigations in the first part of this subsection confirmed the numerical error to scale at least with O(∆t kt+1 ) in the mentioned norms. A well-known property of DG type discretisations is that of superconvergence, i.e. the convergence in time with a higher order (at specific time instances, e.g. nodes; cf. also [4] for a superconvergence result for a CG method of the heat equation.). We investigate numerically whether our method shows such a behaviour in a two-stage procedure: In a first experiment, we perform only time-refinements for a fine mesh and a high order discretisation in space. This experiment allows to estimate the order of (super)-convergence in time, although asymptotically this order will be hard to reach exactly because a saturation from the spatial error will occur. Hence, we validate the first estimate for the superconvergence in a second step by a space-time-refinement study involving the suiting spatial discretisation order. The results for the first study, involving the choices of order (k t , k s ) = (1, 3), (2,5), q = k s and a fixed mesh in space (i s = const) are shown in Figure 13. Whilst in the space-time norm L 2 (L 2 (Ω(t)), 0, T ) the estimate of (6.1) is confirmed, for L 2 (Ω(T )) we observe for the DG method and the CG method with k t = 2, u − u h L 2 (Ω(T )) ≈ O(∆t kt+1.8 ), which in light of the aforementioned saturation effect suggests (6.2) u − u h L 2 (Ω(T )) = O(∆t kt+2 ), for DG, CG with k t = 2.
For CG with k t = 1, the estimate (6.1) is also observed to be sharp, there is no superconvergence. In a second step, we confirm (6.2) to hold with a space-time refinement study involving (k t , k s ) = (1, 2), (2, 3), q = k s , see Figure 14. We mention that theoretical results of findings similar to (6.2) exist for static domains, c.f. e.g. [14,15], and leave a further theoretical study of these phenomena in this particular setting for future research.

6.2.
Moving n-sphere. The second example geometry is a moving circle in 2D/ sphere in 3D/ interval in 1D. These cases involve no change in deformation, but exemplify a moving domain setting. In dimension d, they are described by   The manufactured solution u remains u = cos(π r r0 ) · sin(πt), leading to slightly different expressions for f in each dimension. Moreover, h = 0.5 (is+1) , ∆t = 2 −it−2 , γ J = 0.05, f = 1.1.
6.2.1. On 1D and 3D. In this part, we want to show that the convergence results (6.1) also hold for spatially one-and three-dimensional problems. To this end, we apply the methods of consideration to the mentioned problems and study the decay of the numerical error for k = q, i s = i t . In Figure 15, the DG method is applied in three spatial dimensions and the CG method is considered in 1D. The respective choice of refinement levels and tested discretisation orders reflects the different computational demands. The results in the L 2 (L 2 (Ω(t)), 0, T )-norm are equivalent and shown in the supplement, as well as results for the DG method in 1D and the CG method in 3D. In all cases, the convergence orders of (6.1) are confirmed.

Non-zero entries of the matrices.
For an evaluation of the computational costs, we compare the number of non-zero entries (nzes) of the discrete problem matrix of the respective methods. We fix the setting of the moving circle, k = q, i s = 2. On each time-slice, after the system matrix is assembled, nzes are counted and the minimal and maximal value of this number over all slices is considered.
For the CG method, in comparison to the DG method, nzes will depend on the factor f , as we evaluate in a first study. To this end, we choose i t = 2 and evaluate nzes for the CG method with f ∈ {1, 1.1, 1.5, . . . , 100}, and for comparison for the DG method, for the polynomial degrees k = 1, 2, 3. The results are displayed in Figure 16. First, we observe that over all polynomial degrees the CG method leads   Table 1 Range for non-zero entries (nzes) of system matrix (depends on time step) and absolute L 2 (T )error (err) in comparison for DG, CG and GCC. In this whole table, is = 2.
to less non-zero entries as the DG method, whilst the relative difference is specifically strong for small orders k. This aligns with the elimination of one dof of the time finite element space. Second, this benefit of the CG method is reduced by a high choice of the factor f , as a larger domain of extension will induce more entries in the matrix. Third, the original CG method of Subsection 4.2 performs better than the variant mentioned in Remark 4.3. Hence, to exploit the computational benefits of the CG method fully, we suggest to use the non-tensor-product variant of the ghost penalty on a region as small as possible.
In a second study, we fix f = 1.1 and compare the non-zero entries for all three methods for different time step sizes ∆t and different polynomial orders. The numbers are given in Table 1. Summarizing the results within each cell, the number of non-zero entries decreases from DG to CG and in turn to GCC. In particular for discretisation parameter i t = 3, 5, this holds within a regime of similar absolute numerical errors. Within each row of the table, we note that smaller time steps decrease the numbers of non-zero entries for all methods, as the interface movement per time step decreases also. This effect is stronger for CG and GCC because of the extended regions of the ghost penalty stabilisation. For increasing polynomial order, we again observe a relatively smaller benefit of the CG and GCC method as compared to DG. We conclude that those methods are specifically interesting for orders such as k = 1, 2, 3 from a computational point of view.
6.2.3. Different time integration strategies. In this part we investigate numerical differences between the topology-insensitive and the topology-preserving method of numerical quadrature in time suggested in Section 5. For this purpose, we implement both approaches and compare the respective numerical errors. By and large, we observe competitive numerical results over most of the problem settings introduced so far. This is in agreement with the good results obtained in the literature, cf. [38,18]. However, we observed differences in the following setup with almost no ghost-penalty stabilisation: Consider the DG method applied to the moving interval in 1D with a radius R = 0.505, constant velocity ρ = 0.5t, manufactured right-hand side for u = (x − ρ + R) 2 · (x − ρ − R) 2 , and polynomial order k = q = 4. For this example, we fix i t = 0 and refine in space. From i s = 2 on, the discrete interface is captured exactly so that the method can approximate the solution exactly as a polynomial of order 4. The results of this study are shown in Figure 17. For the method with the topology-preserving time integration, we note that there appear only errors below a reasonable round-off error tolerance such as 10 −12 . Meanwhile, the method of topology-insensitive time integration leads to significant errors when very small stabilisation constants and mesh sizes are used. This is caused by cut configurations where an element is included in the active mesh because of a small cut, but the corresponding matrix entries vanish because of the numerical integration in time. In the presence of a moderate ghost penalty stabilisation, the dofs are still controlled by this stabilisation, otherwise stability issues appear. For an interpretation of this behaviour in terms of numerical analysis, we refer the reader to the comment made about γ J → 0 in Subsection 6.1.3. Moreover, from a numerical analysis point of view we expect that the topology-preserving time integration is important for the geometrical consistency. In the numerical experiments however, we did not observe corresponding issues with the topology-insensitive time integration, at least not in combination with sufficient ghost penalty stabilization. We are not sure if there is an underlying mechanism that can explain this behaviour or if the test case is still too mild and leave further investigations of this for future research. In addition we compared the computational effort for both time integration strategies for the two-dimensional configuration of subsection 6.1. As differences only occur on a small subset of elements in the computation domain, the overall computation time is only marginally affected by the choice of the integration strategy.
We conclude that the topology-preserving time integration represents the discrete geometry properly, and hence is our method of choice, although the topologyinsensitive time integration seems to perform equally well with sufficient stabilisation. Finally, we note that by doubling the order or introducing 10 evenly spaced substeps, the errors for the topology-insensitive time integration approach with degenerate stabilisation can be improved. 9 6.3. Colliding circles. Thirdly, we want to investigate the stability of the method in a topologically challenging example taken from [25]: Two circles collide and separate again afterwards. This setting can be described by the following level set function: w(x, y, t) = (0, −1) T if (y > 0 and t ≤ T /2) or (y < 0 and t > T /2) (0, 1) T if (y ≤ 0 and t ≤ T /2) or (y > 0 and t > T /2).
Here, u 0 is defined to be u 0 (x, y) := sign(y), i.e. 1 in the upper circle and -1 in the lower one. First, we investigate the numerical method with k = q = 4, i s = 2, ∆t = T /16 to generate a reference solution, which is shown on the left-hand side of Figure 18. As one would expect physically, first the circles are transported with the convection field. After they merged, the diffusion successively evens out the concentration difference. By and large, the simulation does not show any numerical instability or large error. As a small detail, we notice that some numerical diffusion will appear shortly before the contact of the circles, when the numerically extended regions merge. This was also observed in [25] for the time stepping method presented there. But for fine meshes or small stabilisation constants, this effect can be kept small. We conclude that the method performs well, despite a challenging topology change, without the need of any specific adaptation.
Second, we consider a more extreme version by setting ∆t = T , so that the whole time interval will be solved in one time step. We show the numerical solutions for the parameter choices k = q = 2, 4 in Figure 18. As the overall topology of the problem is still resolved properly, the method appears capable of handling topologically demanding configurations properly.

Conclusion.
In this final section, we want to summarise the findings of this work and highlight some questions that remain for future research.
We introduced unfitted space-time methods for a problem posed on a moving domain. These are based on discontinuous Galerkin, continuous Galerkin, and higher order continuous Galerkin-Collocation discretisations in time. We further introduced techniques to achieve a high order geometry handling of the space-time domains based on isoparametric mappings of the background mesh. In extensive numerical experiments we discussed and compared several important properties of the proposed methods. Especially, the expectations of robustness and higher order convergence in space and time have been confirmed. Even in topologically challenging configurations, our methods perform reliably. These results motivate the following directions for future investigations: First, the methods have been investigated from a computational point of view in this paper. This suggests to investigate whether the numerically established convergence estimates can be proven from the point of view of numerical analysis. As far as the DG method is concerned, the master's theses [34] and [19] present such error estimates, although in the former an exact handling of the geometries is assumed, which is weakened in the latter in regard to the spatial approximation. This leaves a numerical analysis relative to the spatially and temporally discrete geometry an open task. Moreover, for the CG and GCC method no such analysis has been performed to the author's best knowledge.
Next, the numerical experiments w.r.t. the different numerical integration approaches suggest that the topology-insensitive time integration is also a reasonable and accurate approach and the topology preserving time integration seems to be superfluous in the presence of ghost penalty stabilisation. This behavior is not fully understood and leaves room for future research.
Finally, let us mention that although we constructed a space-time mesh deformation that can be discontinuous between time slabs, it may be possible to construct a continuous mesh deformation that preserves the accuracy of the scheme. This would allow to remove the transfer operation Π u and hence simplify the method and potentially its analysis. The discrete problem is then given as: Find u = u 0 + u init with u 0 ∈ U n h,0 s.t.
Appendix B. Supplementary Plots. In the following, we present results of numerical experiments which we use to draw our final conclusions, which are however very similar to related graphs in the original paper.