<figure><img src="image_1.png" /></figure>

HAL
open science

# Discretization error analysis of a high-order unfitted space-time method for moving domain problems

Fabian Heimann, Christoph Lehrenfeld, Janosch Preuß

► To cite this version:

Fabian Heimann, Christoph Lehrenfeld, Janosch Preuß. Discretization error analysis of a high-order unfitted space-time method for moving domain problems. IMA Journal of Numerical Analysis, 2025, 10.1093/imanum/draf084. hal-05444496

HAL Id: hal-05444496

https://hal.science/hal-05444496v1

Submitted on 6 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Discretization error analysis of a high-order unfitted space–time method for moving domain problems

FABIAN HEIMANN*

Department of Mathematics, University College London, 25 Gordon Street, London WC1H 0AY, UK

*Corresponding author: f.heimann@ucl.ac.uk

CHRISTOPH LEHRENFELD

Institute for Numerical and Applied Mathematics, University of Göttingen, Lotzestraße 16-18, D-37083 Göttingen, Germany

AND

JANOSCH PREUSS

Project-Team Makutu, Inria, University of Pau and Pays de l'Adour, Bâtiment IPRA, Avenue de l'Université, F-64000 Pau, France

[Received on 11 April 2025; revised on 17 July 2025]

We present a numerical analysis of a higher-order unfitted space–time finite element method applied to a convection–diffusion model problem posed on a moving bulk domain. The method uses isoparametric space–time mappings for the geometry approximation of level set domains and has been presented and investigated computationally in Lehrenfeld & Preuß (2023, Geometrically higher order unfitted space–time methods for PDEs on moving domains. *SIAM J. Sci. Comp.* **45**, B139–B165). Recently, in Heimann & Lehrenfeld (2025, Geometry error analysis of a parametric mapping for higher order unfitted space–time methods. *IMA J. Numer. Anal.*) error bounds for the geometry approximation have been proven. In this paper we prove stability and accuracy including the influence of the geometry approximation.

**Keywords:** moving domains; unfitted FEM; isoparametric FEM; space–time; higher order.

## 1. Introduction

Partial differential equations (PDEs) posed on moving domains are ubiquitous in various fields such as physics, engineering, chemistry and biology. Traditional approaches to solve these problems often involve mesh adaptations or remeshings at every time step, which can be computationally expensive, especially for problems with strong deformations or topology changes. In recent years, finite element methods (FEMs) that are geometrically *unfitted* have emerged as an alternative paradigm, offering a promising solution to circumvent the challenges associated with mesh generation. Geometrically unfitted methods, such as CutFEM, finite cell method, XFEM and fictitious domain methods, decouple the computational mesh from the underlying geometry of the problem. Instead, the geometry is prescribed separately, typically implicitly by a level set function (Sethian, 1999; Burman *et al.*, 2015). By avoiding the need for mesh adaptation unfitted methods become particularly attractive for problems involving time-dependent domains.

This paper is the third part of a trilogy—following Heimann *et al.* (2023) and Heimann & Lehrenfeld (2025)—where we focus on the numerical analysis of an unfitted space–time method with higher-order accurate geometry handling for time-dependent smooth domains. Building upon the groundwork laid
in our previous works we extend our analysis to derive error bounds for our space-time unfitted FEM discretization based on a discontinuous Galerkin (DG) methods in time.

In the literature various approaches have been proposed to address the challenges associated with geometrically unfitted FEM and time integration for moving domains. These approaches encompass a range of techniques for numerical integration, stability enforcement and error analysis. We refrain from providing a comprehensive review of the literature, but instead refer the reader to the literature discussion in Heimann et al. (2023) and Heimann & Lehrenfeld (2025).

The foundation for the method discussed in this paper and its error analysis rests upon the principles established in our previous works (Preuß, 2018; Heimann, 2020; Wendler, 2022; Heimann et al., 2023; Heimann & Lehrenfeld, 2025). In Heimann et al. (2023) we introduced a class of geometrically unfitted space-time DG methods for the convection-diffusion problem on moving domains, including methods for high-order geometry approximation, but without a rigorous error analysis. In the Master's theses (Preuß, 2018; Heimann, 2020; Wendler, 2022) the DG-in-time framework has been analysed in detail, but without (full) consideration of the impact of geometric errors. In Heimann & Lehrenfeld (2025) we addressed the challenges of geometric accuracy in unfitted space-time methods, culminating in the derivation of rigorous bounds for the distance between realized and ideal mappings. Building upon these results we now turn our attention to the examination of error sources within the DG-in-time framework.

For comparison note that an error analysis for a similar method for a parabolic surface PDE on a moving surface, based as well on the space-time isoparametric unfitted geometry approximation, was presented recently in Reusken & Sass (2024).

## 1.1 *Structure of the paper*

The remainder of the paper is structured as follows: in Section 2 the model convection-diffusion problem is introduced in detail. In Section 3 the discrete space-time geometries and function spaces are introduced and we summarize the relevant results from Heimann & Lehrenfeld (2025) in the form of a geometry-related assumption. Further, we introduce some relevant discrete objects that are specific for the isoparametrically mapped discretization. Afterwards, in Section 4, the discrete method is introduced. Here, a constant-in-space and piecewise linear-in-time function is split up to be able to put emphasis on global mass conservation. The following Section 5 contains the stability analysis based on an inf-sup procedure and constitutes one of the main contributions of this work. This final result of that section is a Strang lemma bounding the discretization error by approximation and consistency errors. Both are then analysed in Section 6 yielding an *a priori* discretization error bound as our main result. In order to improve the readability of the paper some proofs, which are technical yet rather standard in the field, have been outsourced to Appendix A.

The results of this paper are in line with previous works (Preuß, 2018; Heimann, 2020, 2025; Heimann et al., 2023; Heimann & Lehrenfeld, 2025) on the same discretization. Specifically, in Preuß (2018) no consideration of the higher-order discrete space-time geometries is included, whereas Heimann (2020) contains a partial consideration insofar as a semidiscrete in space computational geometry is considered. We build upon these results and Heimann & Lehrenfeld (2025) in order to present an analysis including the full space-time discrete computational geometry in the following. Relatedly, note that Badia et al. (2023) contains a stability analysis for a related discretization, where robustness against small cuts is realized by element aggregation instead of Ghost-penalty (GP), and a weaker norm is chosen. Also, Lehrenfeld & Reusken (2013) contains a stability analysis in such a weaker norm for a two-domain problem involving a Nitsche stabilization with appropriate average operators, so that Ghost-penalties are not needed. We also mention reference Larson et al. (2024) in which stability is indeed
proven in a stronger norm similar to our approach. The context is however a bit different since this reference deals with a discretization on overlapping meshes. We comment on details in regard to these numerical analyses in Remarks 5.1 and 6.4.

## 2. Model problem

Our model problem of consideration is the convection-diffusion problem on moving domains:

$$ \partial_t u(x, t) + \text{div}(\vec{w}u(x, t)) - \alpha \Delta u(x, t) = f(x, t) \quad \text{in } \Omega(t), \quad (2.1) $$

where $u = u(x, t)$ represents the unknown field, $\alpha > 0$ a diffusion coefficient, $\vec{w} = \vec{w}(x, t)$ a convection velocity, $f(x, t)$ a source term and $\Omega(t)$ a time-dependent, open, bounded and smooth domain in $\mathbb{R}^d$, $d = 1, 2, 3$. The problem is complemented by initial conditions $u_0(x)$. We assume that the convection field $\vec{w}$ is divergence-free (div $\vec{w} = 0$), that $\|\vec{w}\|_{W^{1,\infty}} < \infty$ and that $\vec{w}$ also drives the evolution of the geometry. As boundary condition we use $\nabla u \cdot n_{\partial\Omega} = 0$, which describes that there is no transport across the (possibly moving) boundary. In the following text we rescale the time variable in the equation so that we can set w.l.o.g. $\alpha = 1$. Moreover, we fix the time interval to $t \in [0, T]$ with $T > 0$.

We arrive at the following full strong formulation of the problem for $u(x, t)$:

$$ \partial_t u + \vec{w} \cdot \nabla u - \Delta u = f \quad \text{in } \Omega(t) \text{ for } t \in (0, T], \quad (2.2a) $$

$$ \nabla u \cdot \vec{n}_{\partial \Omega} = 0 \quad \text{on } \partial \Omega(t) \text{ for } t \in (0, T], \quad (2.2b) $$

$$ u(\cdot, 0) = u_0 \quad \text{in } \Omega(0). \quad (2.2c) $$

Due to the no-flux boundary condition (2.2b) there holds a simple balance. With the notation $\bar{v}(t) = \int_{\Omega(t)} v(x, t) dx$ there holds $\bar{u}(t) = \bar{u}(0) + \int_{0}^{t} \bar{f}(\tau) d\tau$.

Let us note that the model problem may appear to be relatively simple as a scalar-valued problem on a moving domain, whose structure is assumed to be given a priori. However, it contains important parts of more involved problems, such as Navier-Stokes equations and related free boundary problems. Hence, the detailed study of the model problem serves as a preparatory step of these more involved problems.

## 3. Discrete space-time geometries and function spaces

In this section we want to introduce the discrete geometries used in our method, together with relevant results from their numerical analysis. These have been shown in detail in the paper Heimann & Lehrenfeld (2025) under regularity assumptions on the level set geometry.

A difficulty investigated in detail in Heimann & Lehrenfeld (2025) is the question of the choice of the blending, which describes the behaviour of the isoparametric mesh deformation between cut elements and the remainder of the mesh. Two blending options, denoted as finite Element blending (FE blending) and smooth blending, have been treated in Heimann & Lehrenfeld (2025). The FE blending that has also been considered in Heimann et al. (2023) is computationally attractive due to its simplicity. However, the analysis for the FE blending is technically more involved as discontinuous deformed meshes and corresponding transfer operations need to be analysed and stronger assumptions on the time steps are necessary. In this work we will work with assumptions on the geometry approximation that are fulfilled by the smooth blending, but not by the FE blending. For weaker assumptions on the geometry approximation including the FE blending we refer to the analysis in (Heimann, 2025, chapter 4).
We fix some notation. First, we assume that there exists a background domain $\tilde{\Omega}$, which is sufficiently large so that for all $t \in [0, T]$, $\Omega(t) \subseteq \tilde{\Omega}$. Moreover, we assume that a family of unfitted triangulations on this background domain, $\{\tilde{\mathcal{T}}_h\}$, is given, with mesh size $h$ and the set of facets $\mathcal{F}_h$. Throughout this paper we assume that the mesh $\tilde{\mathcal{T}}_h$ consists of shape regular simplices.<fnref n="1" /> The time interval $[0, T]$ is assumed to be subdivided into time steps of equal width $\Delta t$, $I_n = (t_{n-1}, t_n]$, $t_n - t_{n-1} = \Delta t$, $n = 1, \dots, N$, $t_0 = 0$, $t_N = T$ for simplicity. We fix an order of spatial geometry approximation order $q_s$ and a temporal geometry approximation order $q_t$ and define on each time interval $I_n$

$$W_h^{n,(q_s,q_t)} := V_h^{q_s} \otimes \mathcal{P}^{q_t}([t_{n-1}, t_n]), \quad V_h^{q_s} := \{v \in H^1(\tilde{\Omega}) | v|_{T} \in \mathcal{P}^{q_s}(T) \forall T \in \tilde{\mathcal{T}}_h\}. \quad (3.1)$$

In the remainder we will use $q = (q_s, q_t)$ for the geometrical accuracy in space and time, whereas $k = (k_s, k_t)$ refers to the orders of the discrete spaces used for approximating the solution.

We assume that the moving domain is given implicitly by a sufficiently smooth<fnref n="2" /> levelset function $\phi: \tilde{\Omega} \times [0, T] \to \mathbb{R}$, i.e., for all $t \in [0, T]$,

$$\Omega(t) = \{x \in \tilde{\Omega} | \phi(x, t) < 0\}. \quad (3.2)$$

We introduce $Q$ for the space-time version corresponding to $\Omega(t)$, i.e.,

$$Q := \{(x, t) \in \tilde{\Omega} \times [0, T] | \phi(x, t) < 0\} = \bigcup_{t \in (0, T)} \Omega(t) \times \{t\}. \quad (3.3)$$

To separate spatial and temporal continuity we will use the following notation for space-time Sobolev spaces (which are also called $t$-anisotropic Sobolev spaces):

$$H^{r,s}(Q) = H_{\square}^{r,s}(Q) := \{u | \partial_t^p D^\alpha u \in L^2(Q), p, q \in \mathbb{N}, q = |\alpha|, \frac{q}{r} + \frac{p}{s} \le 1\}. \quad (3.4)$$

$$H_{\square}^{r,s}(Q) := \{u | \partial_t^p D^\alpha u \in L^2(Q), p, q \in \mathbb{N}, q = |\alpha| \le r, p \le s\}. \quad (3.5)$$

Moreover, we abbreviate $H^{r,r} := H^r$. Within the assumption about the discrete isoparametric unfitted geometry two discrete approximations of the levelset function will be used in order to control the discrete geometry, a spatially element-wise linear approximation $\phi^{\text{lin}}$, and a higher-order approximation $\phi_h$. (See also Heimann & Lehrenfeld (2025) about the detailed assumptions on and role of the latter.) The element-wise linear approximation $\phi^{\text{lin}}$ corresponds to an approximation of $Q$, called $Q^{\text{lin}}$, which has the computationally highly convenient feature that restricted to one time instance (within one time slice); the domain can be decomposed into a finite number of simplices. In that way numerical integration on the linear reference geometry can be performed straightforwardly by an affine mapping of the standard Gaussian simplex rules, cf. (Heimann & Lehrenfeld, 2019; Heimann et al., 2023). In order to make use of these integration rules while ensuring spatially high-order accurate approximation at the same time an isoparametric mapping of the mesh is applied. On an abstract level it is related to strategies that have been developed for the boundary approximation with standard FEM (Ciarlet & Raviart, 1972; Bernardi, 1989). An isoparametric mapping for spatial unfitted FEM has been developed in Lehrenfeld (2016)

<footnote n="1">¹ For quadrilateral meshes, in general, most of the constructions should hold with according adaptations, which we however leave aside. In Heimann & Lehrenfeld (2019) the quadrilateral case was discussed for a spatial problem.</footnote>

<footnote n="2">² See Heimann & Lehrenfeld (2025) for technical details of the regularity assumed/needed in relation to $\phi$.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Mappings involved in the discrete space-time geometry construction and accuracy analysis.</figcaption></figure>

and Lehrenfeld & Reusken (2018) and generalized to the space-time setting in Heimann & Lehrenfeld (2025). Following the notation from these works we denote by $\Theta_h^{st}$ the isoparametric mapping, so that $Q^h := \Theta_h^{st}(Q^{\text{lin}})$ will be the higher-order discrete geometry approximation with spatial boundary $\partial_s Q^h$.

For the purposes of the numerical analysis it is of relevance to assess the quality of approximation of $\Theta_h^{st}$. To this end, an ideal counterpart $\Psi^{st}$ is introduced (Heimann & Lehrenfeld, 2025), which maps $Q^{\text{lin}}$ to the exact geometry $Q$. The approximation quality can then be expressed mathematically in terms of norms on the difference between $\Theta_h^{st}$ and $\Psi^{st}$. Finally, this difference can be interpreted equivalently by bounding the mapping $\Phi^{st} := \Psi^{st} \circ (\Theta_h^{st})^{-1}$, cf. Fig. 1 for a sketch.

In this work we will use $\lesssim (\gtrsim)$ to mean that the left-hand side is bounded above (below) by a constant times the right-hand side for constants that do not depend on $h$ or $\Delta t$, or the position of a cut within any element or geometry. We use $a \simeq b$ if $a \lesssim b$ and $b \lesssim a$ holds simultaneously.

In the following text we assume that a space-time isoparametric strategy is used for the geometry approximation that yields the following properties:

**ASSUMPTION 3.1** (Approximation of unfitted geometry). For each choice of approximation orders in space and time, $q_s, q_t \ge 1$, there exists a function $\phi^{\text{lin}}: \tilde{\Omega} \times [0, T] \to \mathbb{R}$, $\phi^{\text{lin}}|_{\tilde{\Omega} \times I_n} \in W_h^{n,(1,q_t)}$ for all $n = 1, \dots, N$, which is continuous along time slice boundaries, i.e., $\lim_{s \searrow t_n} \phi^{\text{lin}}(\cdot, s) = \lim_{s \nearrow t_n} \phi^{\text{lin}}(\cdot, s)$ ($n = 1, \dots, N-1$), as well as functions $\Psi^{\text{st}}: \tilde{\Omega} \times [0, T] \to \tilde{\Omega} \times [0, T]$, $\Theta_h^{\text{st}}: \tilde{\Omega} \times [0, T] \to \tilde{\Omega} \times [0, T]$ such that $Q^n = \{(x, t) \in \tilde{\Omega} \times I_n | \phi(x, t) < 0\} = \Psi^{\text{st}}(Q^{\text{lin},n})$, where $Q^{\text{lin},n} := \{(x, t) \in \tilde{\Omega} \times I_n | \phi^{\text{lin}}(x, t) < 0\}$. Further, these functions are sufficiently regular so that $\Phi^{\text{st}} := \Psi^{\text{st}} \circ (\Theta_h^{\text{st}})^{-1}$ is well defined. Then, these functions satisfy the following proximity, smallness and continuity properties:

$$ \| \partial_t^{l_t} D_x^{l_s} (\text{Id} - \Phi^{\text{st}}) \|_{\infty, \tilde{\Omega} \times I_n} \lesssim h^{q_s+1-l_s} + \Delta t^{q_t+1-l_s}, \quad l_t+l_s \in \{0, 1\}, \quad (3.6a) $$

$$ \|1 - \det D \Phi^{\text{st}}\|_{\infty, \tilde{\Omega} \times I_n} \lesssim h^{q_s} + \Delta t^{q_t+1}, \quad (3.6b) $$

$$ \|1 - J_{\Phi}\|_{\infty, \partial_s Q^h} \lesssim h^{q_s} + \Delta t^{q_t}, \quad (3.6c) $$

$$ \| \partial_t^{l_t} D_x^{l_s} (\text{Id} - \Theta_h^{\text{st}}) \|_{\infty, \tilde{\Omega} \times I_n} \lesssim h^{2-l_s} + \Delta t^{q_t+1-l_s}, \quad l_t+l_s \in \{0, 1\}, \quad (3.6d) $$

$$ \lim_{s \searrow t_n} \Theta_h^{\text{st}}(\cdot, s) = \lim_{s \nearrow t_n} \Theta_h^{\text{st}}(\cdot, s) := \Theta_h^{\text{st}}(t_n), \quad (3.6e) $$

$$ \sum_{l_s=0,1} \|D_x^{l_s} \Psi^{\text{st}}\|_{H^{0,q_t+1}(\tilde{\Omega} \times I_n)} + \sum_{l_t=0,1} \|\partial_t^{l_t} \Psi^{\text{st}}\|_{H^{q_s+1,0}(\tilde{\Omega} \times I_n)} \lesssim 1. \quad (3.6f) $$

Here, $J_{\Phi}$ denotes the ratio of surface measures for $\partial \tilde{\Omega}(t)$ and $\partial \tilde{\Omega}_h(t)$ for $t \in [0, T]$. In Heimann & Lehrenfeld (2025) the properties of Assumption 3.1 have been shown under certain assumptions, cf.
especially corollaries 5.1–5.3 in that paper. These assumptions include the choice of the blending step in the isoparametric mapping approach and the space–time level set geometry, specifically in the form of the regularity of the space–time level set function and—for the anisotropy of the error bounds—on the quality of semidiscrete approximations of the level set function in space and time, respectively. We do not explicate the details here, but refer to Heimann & Lehrenfeld (2025) instead. Moreover, let us summarize that for sufficiently smooth domains these assumptions are reasonable.

**ASSUMPTION 3.2.** Since the analysis in Heimann & Lehrenfeld (2025) is essentially asymptotic we will henceforth assume that the mesh size $h$ and time step $\Delta t$ are sufficiently small. That is, in all subsequent statements $h$ (or $\Delta t$) is assumed to be smaller than a certain constant, which may vary, depending on the context of each statement, but is independent of cut configurations and $\Delta t$ (or $h$, respectively).

An important implication of Assumption 3.1 is the following:

**LEMMA 3.1** Let $\hat{u}$ be a function defined on $\tilde{\Omega} \times [0, T]$ (or a tensor-product subregion), which is sufficiently smooth s.t. all expressions in this Lemma are well defined. Let $u := \hat{u} \circ (\Theta_h^{\text{st}})^{-1}$; then for all $S \subseteq \tilde{\Omega}$ with $\text{vol}_d(S) > 0$ and $t \in [0, T]$ there holds

$$ \|D^l \hat{u}(\cdot, t)\|_S \simeq \|D^l u(\cdot, t)\|_{\Theta_h^{\text{st}}(S,t)}, \quad l_s \in \{0, 1\} \quad \text{and} \quad \|\hat{u}(\cdot, t)\|_{\partial S} \simeq \|u(\cdot, t)\|_{\Theta_h^{\text{st}}(\partial S,t)}. \quad (3.7) $$

We introduce the following convention for denoting functions: while the space–time to space–time mappings (as above) are called $\Theta_h^{\text{st}}$, $\Psi^{\text{st}}$ and $\Phi^{\text{st}}$, the respective space–time to space mappings bear the same name without index st, i.e., $\Theta_h$, $\Psi$ and $\Phi$ are defined such that

$$ \Theta_h^{\text{st}}(x, t) = (\Theta_h(x, t), t)^T, \quad \Psi^{\text{st}}(x, t) = (\Psi(x, t), t)^T, \quad \Phi^{\text{st}}(x, t) = (\Phi(x, t), t)^T. $$

Based on these functions we introduce the following discrete domains:

$$ \Omega^{\text{lin}}(t) := \{x \in \tilde{\Omega} | \phi^{\text{lin}}(x, t) < 0\}, \quad \Omega^h(t) := \Theta_h(\Omega^{\text{lin}}(t), t), \quad (3.8a) $$

$$ Q^{\text{lin},n} := \bigcup_{t \in I_n} \Omega^{\text{lin}}(t) \times \{t\}, \quad Q^{h,n} := \Theta_h^{\text{st}}(Q^{\text{lin},n}), \quad Q^h = \bigcup Q^{h,n}. \quad (3.8b) $$

Next, we introduce notation for the smallest set of space–time-tensor-product elements with nonempty intersection with $Q^{\text{lin},n}$, those that are contained in $Q^{\text{lin},n}$ and corresponding mapped versions:

$$ Q_{\mathcal{E}}^{\text{lin},n} := \{T \times I_n \mid (T \times I_n) \cap Q^{\text{lin},n} \neq \emptyset, T \in \mathcal{T}_h\}, \quad Q_{\mathcal{E}}^{h,n} = \Theta_h^{\text{st}}(Q_{\mathcal{E}}^{\text{lin},n}), \quad (3.9a) $$

$$ Q_{\mathcal{F}}^{\text{lin},n} := \{T \times I_n \mid (T \times I_n) \subseteq Q^{\text{lin},n}, T \in \mathcal{T}_h\}, \quad Q_{\mathcal{F}}^{h,n} = \Theta_h^{\text{st}}(Q_{\mathcal{F}}^{\text{lin},n}). \quad (3.9b) $$

The related merely spatial element collections are introduced as

$$ \Omega_{\mathcal{E}}^{\text{lin},I_n} := \{T \in \mathcal{T}_h | T \times I_n \in Q_{\mathcal{E}}^{\text{lin},n}\}, \quad \Omega_{\mathcal{E}}^{h,I_n}(t) = \Theta_h^{\text{st}}(\Omega_{\mathcal{E}}^{\text{lin},I_n}, t), \quad (3.9c) $$

$$ \Omega_{\mathcal{F}}^{\text{lin},I_n} := \{T \in \mathcal{T}_h | T \times I_n \in Q_{\mathcal{F}}^{\text{lin},n}\}, \quad \Omega_{\mathcal{F}}^{h,I_n}(t) = \Theta_h^{\text{st}}(\Omega_{\mathcal{F}}^{\text{lin},I_n}, t), \quad (3.9d) $$

$t \in [0, T]$. With abuse of notation we will use the same notation for the sets of elements $\Omega_{\mathcal{E}}^{\cdot\cdot\cdot}, \Omega_{\mathcal{F}}^{\cdot\cdot\cdot}, Q_{\mathcal{E}}^{\cdot\cdot\cdot}$ and $Q_{\mathcal{F}}^{\cdot\cdot\cdot}$, as well as for the corresponding domains, depending on the context.
<figure><img src="image_4.png" /><figcaption>FIG. 2. Illustration of an example for $\Omega_{\mathcal{E}}^{\text{lin},I_n}$ and $\Omega_{\mathcal{J}}^{\text{lin},I_n}$ on the undeformed mesh for some time interval $I_n$ with $\partial\Omega^{\text{lin}}(t_{n-1})$ and $\partial\Omega^{\text{lin}}(t_n)$ illustrated to indicate the cut topologies.</figcaption></figure>

In Fig. 2 we give an example of the regions $\Omega_{\mathcal{E}}^{\text{lin},I_n}$ and $\Omega_{\mathcal{J}}^{\text{lin},I_n}$ for an example case. We refer the reader to Heimann et al. (2023) for a sketch of the corresponding space–time regions.

We introduce a variant of the discrete spaces on the whole mesh (cf. Equation (3.1)): for some space–time domain $\mathcal{Q} \subseteq \tilde{\Omega} \times [0, T]$ we want to introduce the subspace with activated degrees of freedom only on $\mathcal{Q}$, denoting $k = (k_s, k_t)$

$$W_{h,\text{cut}}^{n,k}(\mathcal{Q}) = \{v \in W_h^{n,k} \mid v \text{ vanishes outside of } \mathcal{Q} \text{ as element of } W_h^{n,k}\}. \quad (3.10)$$

We glue together these restricted spaces from all time-steps to obtain the following set of discrete functions on the linear reference configuration:

$$W_h^{\text{lin}} = \left\{ v \in L^2 \left( \bigcup_{n=1}^{N} Q_{\mathcal{E}}^{\text{lin},n} \right) \mid \forall n = 1, \dots, N: v|_{Q_{\mathcal{E}}^{\text{lin},n}} \in W_{h,\text{cut}}^{n,k}(Q_{\mathcal{E}}^{\text{lin},n}) \right\}. \quad (3.11)$$

Finally, we apply the isoparametric mapping:

$$W_h = W_h^{\text{lin}} \circ (\Theta_h^{\text{st}})^{-1}. \quad (3.12)$$

Note that discrete functions in $W_h$ are allowed to be discontinuous along time slice boundaries. We define the jump of discrete functions across time slice boundaries for $u \in W_h$ at some time-slice boundary $t_n$ as

$$[u]^n := u_+ - u_-, \text{ where } u_+ = \lim_{s \searrow t_n} u(s), u_- = \lim_{s \nearrow t_n} u(s). \quad (3.13)$$

For the purposes of the upcoming analysis we also want to introduce a variant of the time derivative $\partial_t$, namely a partial derivative in time in mesh coordinates, for argument functions $u$ defined on $Q^{hn}$.
**DEFINITION 3.1** (Derivative in mesh direction). Let $u: Q^{h,n} \to \mathbb{R}$. Then, its time derivative in mesh coordinates is defined as

$$ \partial_t^\Theta u = \partial_t u + (\partial_t \Theta_h^\text{st} \circ (\Theta_h^\text{st})^{-1}) \cdot \nabla u. \quad (3.14) $$

From the definition of $\partial_t^\Theta$ and the calculations above we obtain

**COROLLARY 3.1.** Denoting $\hat{u} := u \circ \Theta_h^\text{st}$ we have for $S \subseteq \tilde{\Omega}$ and $t \in [0, T]$

$$ \partial_t^\Theta u = (\partial_t \hat{u}) \circ (\Theta_h^\text{st})^{-1}, \text{ which implies } \|\partial_t^\Theta u\|_{\Theta_h(S,t)} \simeq \|\partial_t \hat{u}\|_S. \quad (3.15) $$

Further, for $u \in W_h$ we find $\partial_t^\Theta u \in W_h$ as $\partial_t \hat{u} \in W_h^\text{lin}$ for $\hat{u} \in W_h^\text{lin}$.

## 4. The discrete method

We now proceed by introducing the discrete discontinuous Galerkin in time (DG) method. In general, the presentation follows that of the presentation for computational purposes in Heimann et al. (2023). However, while we used a time-slice wise formulation there, which is more instructive for the implementation, we now directly give a formulation over all time slice contributions at once, which makes it more accessible for the numerical analysis. We first define the bilinear form corresponding to the PDE operator in (2.2):

$$ B_h(u, v) = (\partial_t u + \vec{w} \cdot \nabla u, v)_{Q^h} + (\nabla u, \nabla v)_{Q^h} + \sum_{n=1}^{N-1} ([u]^n, v_+^n)_{\Omega^h(t_n)} + (u_+^0, v_+^0)_{\Omega^h(0)}. $$

Here, we use $(\cdot, \cdot)_{Q^h}$ as a shorthand notation for $\sum_{n=1}^N (\cdot, \cdot)_{Q^{h,n}}$. The DG upwind penalty terms in time are included, which penalize jumps in time and ensure the transport of information between the time slabs. In addition, we define a variant that can be obtained from partial integration from $B_h(\cdot, \cdot)$ for exact geometry handling:

$$ B_{mc}(u, v) = -(u, \partial_t v + \vec{w} \cdot \nabla v)_{Q^h} + (\nabla u, \nabla v)_{Q^h} - \sum_{n=1}^{N-1} (u_-^n, [v]^n)_{\Omega^h(t_n)} + (u_+^N, v_-^N)_{\Omega^h(T)}. $$

The index mc shall indicate a mass conservation property. An important property of $B_h(\cdot, \cdot)$ and $B_{mc}(\cdot, \cdot)$ is obtained by evaluating the sum for $v = u$:

$$ (B_h + B_{mc})(u, u) = 2 \|\nabla u\|_{Q^h}^2 + \sum_{n=1}^{N-1} \|[u]^n\|_{\Omega^h(t_n)}^2 + \|u_+^0\|_{\Omega^h(0)}^2 + \|u_-^N\|_{\Omega^h(T^-)}^2. \quad (4.1) $$

In accordance to the previous definitions the right-hand side $f_h$ is given as

$$ f_h(v) := \sum_{n=1}^{N} (f^n, v)_{Q^{h,n}} + (u_0^e, v_+^0)_{\Omega^h(0)}. \quad (4.2) $$
<figure><img src="image_5.png" /><figcaption>FIG. 3. Details of the definition of the jump operator on the curved elements.</figcaption></figure>

Here, and in the following text, we use the construction of a degree-independent continuous Sobolev extension as in Stein (1971), which we indicate by a superscript e, e.g., $u_0^e$ is such an extension of $u_0$ from $\Omega(0)$ to $\mathbb{R}^d$ and $f^e$ denotes a smooth extension of $f$ from $Q$ to $Q^h$. After introducing a cut stabilization next we will discuss the mass conservation property together with the geometry approximation error, which renders $B_{mc}(\cdot, \cdot)$ different from $B_h(\cdot, \cdot)$.

## 4.1 *Ghost-penalties*

In order to ensure stability in the context of cut elements a GP stabilization is commonly used in unfitted methods (Burman, 2010; Burman et al., 2015). There are different variants of this technique, and we will follow Preuß (2018) and Heimann (2020) in using the direct variant. Other variants may certainly also work in practice. Let us note however that in the context of isoparametric mappings only the *derivative jump* GP and the *direct* GP variant have been analysed in Lehrenfeld (2017) and Lou & Lehrenfeld (2022) and only for the *stationary* domain case so far.

First, we define the transformed facet patch region for some facet $F$ from the uncurved mesh, $F = T_1 \cap T_2, \omega_F = T_1 \cup T_2, \omega_F^h(t) := \Theta_h(\omega_F, t)$. In addition, we need a jump operation on these curved elements. Fix a time $t \in I_n$ and let us assume a spatial discrete function $u = \hat{u} \circ (\Theta_h)^{-1}$ on $\omega_F^h$ is given. We are interested in the 'volumetric jump' for some point $x$ in the facet patch. To this end, we extend the mapped polynomials from both sides of the facet to $\mathbb{R}^d$ and subtract them. Let $\hat{u}_1 = \mathcal{E}^p(u \circ \Theta_h)|_{T_1}$ and $\hat{u}_2 = \mathcal{E}^p(u \circ \Theta_h)|_{T_2}$, where $\mathcal{E}^p : \mathcal{P}^{k_s}(S) \to \mathcal{P}^{k_s}(\mathbb{R}^d)$ is the canonical extension of polynomials of degree $k_s$ from a domain $S$ to $\mathbb{R}^d$.

Then, $u_1 = \hat{u}_1 \circ (\mathcal{E}^p \Theta_h|_{T_1})^{-1}$ and $u_2 = \hat{u}_2 \circ (\mathcal{E}^p \Theta_h|_{T_2})^{-1}$ are well defined in a neighbourhood of $\omega_F^h(t)$ and we define

$$[u]_{\omega_F^h(t)} := u_1 - u_2 = \mathcal{E}^p(u \circ \Theta_h)|_{T_1} \circ (\mathcal{E}^p \Theta_h|_{T_1})^{-1} - \mathcal{E}^p(u \circ \Theta_h)|_{T_2} \circ (\mathcal{E}^p \Theta_h|_{T_2})^{-1}. \quad (4.3)$$

This situation is sketched in Fig. 3.

With that notation we introduce the GP terms $J(u, v) = \sum_{n=1}^{N} j_{h}^{n}(u, v)$ with

$$j_h^n(u, v) = \sum_{F \in \mathcal{F}_R^n} j_F^n(u, v), \quad j_F^n(u, v) := \iint_{h=1}^{n} \iint_{\mathcal{B}_F^h(t)} \frac{\tilde{\gamma}_J}{h^2} [u]_{\omega_F^h(t)} [v]_{\omega_F^h(t)} dx dt. \quad (4.4)$$
Here, $\tilde{\gamma}_J = (1 + \Delta t/h) \gamma_J$ with $\gamma_J > 0$ and the set of facets $\mathcal{F}_R^n \subseteq \mathcal{F}_h$ is chosen to ensure Assumption 5.4. We state a preliminary discrete problem: find $u_h \in W_h$, s.t.

$$B_h(u_h, v_h) + J(u_h, v_h) = f_h(v_h) \quad \forall v_h \in W_h. \qquad (4.5)$$

For the mass conserving bilinear form we formulate the variant: find $u_{mc} \in W_h$, s.t.

$$B_{mc}(u_{mc}, v_h) + J(u_{mc}, v_h) = f_h(v_h) \quad \forall v_h \in W_h. \qquad (4.6)$$

## 4.2 Global conservation and a space splitting

Above, we introduced the bilinear forms $B_h(\cdot, \cdot)$ and $B_{mc}(\cdot, \cdot)$. Let $u_{mc}$ be the discrete solution to (4.5) if $B_h(\cdot, \cdot)$ is replaced by $B_{mc}(\cdot, \cdot)$, assuming—for now—that it exists. Further, let $\bar{u}_{mc}^\ell = \int_{\Omega_h^\ell(t)} u_{mc}^\ell(x, t_{\ell,-}) dx, \ell = 0, \dots, N$ be the notation for the total mass at time $t_\ell$. Then, for one time step $n \in \{1, \dots, N\}$, by testing with $1_{I_n} \in W_h$, where $1_{I_n}(x, t) = 1$ for $t \in I_n$ and $1_{I_n}(x, t) = 0$ else, we obtain the following balance law for the total mass:

$$\bar{u}_{mc}^n = \bar{u}_{mc}^{n-1} + \int_{t_{n-1}}^{t_n} \bar{f}_h(t) \, dt, \quad n = 1, \dots, N, \qquad (4.7)$$

where $\bar{f}_h(t) := \int_{\Omega_h^h(t)} f(x, t) \, dx$. This implies $\bar{u}_{mc}^N = \bar{u}_{mc}^0 + \int_0^T \bar{f}_h(t) \, dt$.

This conservation property is lost for the discretization (4.5) only due to geometry approximation errors, which we bound in the next lemma.

**LEMMA 4.1** For all $u, v \in H^1(\bigcup_{n=1}^N Q^{hn})$ there is

$$|B_h(u,v) - B_{mc}(u,v)| \le \iint_{\partial \Omega_h^h(t)} |\vec{w} \cdot n_{\partial \Omega_h^h} - \mathcal{Y}_h| uv \, ds \, dt, \qquad (4.8a)$$

where $\mathcal{Y}_h$ is the velocity in the direction $n_{\partial \Omega_h}$ that describes the motion of the discrete geometry. Furthermore, there is a constant $C_{4,1}$ (depending on $\|\vec{w}\|_{W^{1,\infty}} < \infty$) such that

$$\|\vec{w} \cdot n_{\partial \Omega_h^h} - \mathcal{Y}_h\|_{L^\infty(\partial_s Q^h)} \le C_{4,1} (h^{qs} + \Delta t^{qt}). \qquad (4.8b)$$

*Proof.* We only sketch the proof here and refer to (Heimann, 2025, lemma 4.3) for details. Applying partial integration w.r.t. time on each time slab we obtain as a difference only a space-time boundary integral on the lateral boundary as in (4.8a). As $\mathcal{Y}_h$ solely depends on the discrete space-time normal $n_h$ the result (4.8a) can be obtained from geometrical error bounds for $n_h$ and proximity bounds of $\Phi^{st}$ for the identity, cf. (Heimann & Lehrenfeld, 2025, theorem 5.8 and lemma 5.12). $\square$

We observe that the error term in (4.8a) stems from the fact that the space-time velocity $(\vec{w}, 1)^T$ that we use for the space-time convection operator $(\vec{w}, 1) \cdot (\nabla, \partial_t)^T = \partial_t \cdot + \vec{w} \nabla \cdot$ is not exactly tangential to the discrete space-time geometry, which perturbs the perfect global balance. Due to these geometry errors for the discretization (4.5) we only have

$$\bar{u}_h^n = \bar{u}_h^{n-1} + \int_{t_{n-1}}^{t_n} \bar{f}_h(t) \, dt + \underbrace{B_h(u_h, 1_{I_n}) - B_{mc}(u_h, 1_{I_n})}_{\le \Delta t(h^{qs} + \Delta t^{qt})} \qquad (4.9)$$
Note that due to the no-flux boundary condition the only mechanism that allows to control functions that are constant in space, e.g., $\int_S \bar{\phi} u_h(x, t) dx$ is the global mass balance. This is in contrast to the case of other boundary conditions, e.g., Dirichlet- or Robin-type boundary conditions, which would allow to control these constants more directly. In the analysis we will therefore apply a special treatment for the mean value at the end of each time slice. We introduce

$$ \overline{\overline{W}}_h = \{v_h \in W_h | v_h(x, t) = \phi(t), \phi(t) \in C^0([0, T]); \phi|_{I_n} \in \mathcal{P}^1(I_n)\} \quad (4.10) $$

as the space of functions that are constant in space and piecewise linear in time. Note that $\overline{\overline{W}}_h$ has $N+1$ degrees of freedoms that correspond to the mean value at the $N+1$ times $t_n$. Let $\overline{\overline{u}}_h \in \overline{\overline{W}}_h$ be a piecewise linear approximation of the mean value on $\Omega^h(t_n)$ for $n = 1, \dots, N$ based on the initial data and the source term on the discretely approximated domain. We define

$$ \tilde{W}_h = \{v \in W_h | (v^n, 1)_{\Omega^h(t_n)} = \int_S \phi(t_n) v^n = 0\} \quad (4.11) $$

as the space of functions with mean value zero on $\Omega^h(t_n)$ for $n = 1, \dots, N$, s.t. $\bar{v}_h^n = 0$, $n = 1, \dots, N$ for all $v_h \in \tilde{W}_h$. Note that $W_h = \tilde{W}_h + \overline{\overline{W}}_h$. We then formulate a discrete problem for given $\overline{\overline{u}}_h \in \overline{\overline{W}}_h$ as follows: find $u_h \in \tilde{W}_h + \overline{\overline{u}}_h$ s.t.

$$ (B_h + J)(u_h, v_h) = f_h(v_h) \quad \forall v_h \in \tilde{W}_h. \quad (4.12) $$

Let $f_S v dx = \int_S v dx/|S|$ denote the mean value of $v$ over a set $S$ with measure $|S|$; then a straightforward choice for $\overline{\overline{u}}_h$ is $\overline{\overline{u}}_h = \int_{\Omega^h(t_n^-)} u_h(t_n^-) dx$ for $u_h$, the solution to (4.5) (assuming it exists) rendering (4.12) equivalent to (4.5). In practice, numerical experiments indicate that (4.5) has a unique solution and hence, the splitting $u_h = \tilde{u}_h + \overline{\overline{u}}_h$ can be regarded merely as an analysis artefact, cf. Section 4.3 where we discuss the realization of the different methods.

In the remainder we discuss the rigorous approach assuming the splitting. One reasonable choice to prescribe the total mass is based on the balance law considerations from above by setting

$$ \underbrace{|\Omega^h(t_n)| \cdot \overline{\overline{u}}_h(t_n)}_{=\overline{\overline{u}}_h^n} = \underbrace{|\Omega^h(t_{n-1})| \cdot \overline{\overline{u}}_h(t_{n-1})}_{=\overline{\overline{u}}_h^{n-1}} + \int_{t_{n-1}}^{t_n} \int_{t_n}^{t} \overline{f}_h(t) dt, \quad n = 1, \dots, N. \quad (4.13) $$

For the subsequent analysis we choose $\overline{\overline{u}}_h$ as in (4.13) and consider solving for the mean value free part through the following problem: find $\tilde{u}_h \in \tilde{W}_h$ such that

$$ (B_h + J)(\tilde{u}_h, v_h) = f_h(v_h) - B_h(\overline{\overline{u}}_h, v_h) \quad \forall v_h \in \tilde{W}_h, \quad (4.14) $$

where we note that in $B_h(\overline{\overline{u}}_h, v_h)$ only the time derivative integral remains (spatial gradients and time jumps vanish). Especially, $J(\overline{\overline{u}}_h, v_h) = 0$ for $v_h \in \tilde{W}_h$.
## 4.3 Realization of methods

As we have shown at the beginning of the previous subsection solving the discrete problem on $W_h$ with the mass conserving bilinear form $B_{mc}(\cdot, \cdot)$ actually yields (4.13) for the mean value. (4.6) can hence be analysed using the splitting without the need to explicitly realize the splitting in an implementation. Note however that the solution will be different from (4.14) as $B_h \neq B_{mc}$. In the analysis below we will treat the mass conserving discretization as a variant for which we conclude a priori error bounds, cf. Corollary 6.1.

The main method for the analysis however shall be (4.14) with (4.13). To realize it in practice we do not explicitly want to work with $\tilde{W}_h$ in an implementation and, hence, one will usually work with $W_h$ in practice. In the remainder of this section let us briefly discuss how to implement the method and relate this to a consistent penalty approach as it has been considered in Heimann (2025).

Let $\overline{\tilde{w}}_h^n = \overline{\tilde{u}}_h(t_n)$ be given as in (4.13) and $\overline{\tilde{w}}_h \in \overline{\tilde{W}}_h$ be the corresponding piecewise linear function. Denote by $(w_h, \lambda) \in W_h \times \mathbb{R}^N$ the solution to

$$ (B_h + J)(w_h, v_h) + \sum_{n=1}^{N} (v_{h,-}^n, \lambda_n)_{\Omega^h(t_n)} = f_h(v_h) \quad \forall v_h \in W_h, \quad (4.15a) $$

$$ \sum_{n=1}^{N} (w_{h,-}^{n}, \mu_{n})_{\Omega^{h}(t_{n})} = \sum_{n=1}^{N} (\overline{\tilde{w}}_{h}^{n}, \mu_{n})_{\Omega^{h}(t_{n})} \quad \forall \mu \in \mathbb{R}^{N}. \quad (4.15b) $$

Then, $f_h \int_{\Omega^h(t_n)} w_h dx = \overline{\tilde{u}}_h(t_n)$ for $n = 1, \dots, N$ and $\tilde{w}_h = w_h - \overline{\tilde{u}}_h \in \overline{\tilde{W}}_h$ solves (4.14). Considering standard/saddlepoint theory we find that problem (4.15) is well-posed iff $(B_h + J)(\cdot, \cdot)$ is stable (and continuous) on $\tilde{W}_h$, i.e., if (4.14) is well-posed (the operator to the scalar constraint is obviously surjective).

In Heimann (2025) a penalty approach w.r.t. the mean value has been considered, which takes the following form: with $I(v, w) := \sum_{n=1}^{N} (v_{-}^n, 1)_{\Omega^h(t_n)} (w_{-}^n, 1)_{\Omega^h(t_n)}$ and for $K > 0$ find $u_h^K \in W_h$ s.t.

$$ (B_h + J + K \cdot I)(u_h^K, v_h) = f_h(v_h) + K \cdot I(\overline{\tilde{w}}_h, v_h). \quad (4.16a) $$

We note that an implementation of $I(\cdot, \cdot)$ should exploit the rank-one nature of this term to avoid setting up dense system matrices. Alternatively, the problem can be reformulated via a scalar auxiliary variable as follows: find $(u_h^K, \lambda) \in W_h \times \mathbb{R}^N$, s.t.

$$ (B_h + J)(u_h^K, v_h) + \sum_{n=1}^{N} (v_{h,-}^n, \lambda_n)_{\Omega^h(t_n)} = f_h(v_h) \quad \forall v_h \in W_h, \quad (4.16b) $$

$$ \sum_{n=1}^{N} (u_{h,-}^{K,n}, \mu_n)_{\Omega^h(t_n)} - K^{-1} \sum_{n=1}^{N} \lambda_n \mu_n = \sum_{n=1}^{N} (\overline{\tilde{w}}_h^n, \mu_n)_{\Omega^h(t_n)} \quad \forall \mu \in \mathbb{R}^N, \quad (4.16c) $$

where we note that (4.16c) implies $\lambda_n = K \cdot (u_{h,-}^{K,n} - \overline{\tilde{w}}_h^n, 1)_{\Omega^h(t_n)}$. Further, we note that (4.16a), when restricted to $\tilde{W}_h$, ends up in (4.14) (for $\overline{\tilde{u}}_h = \overline{\tilde{u}}_h^K$). We conclude that in all these methods the understanding of (4.14), especially in view of stability, plays a crucial role.

**REMARK 4.1** (Literature on related mass conserving variant methods). Counterparts of the method variant (4.6) in slightly different set-ups have been discussed in the literature. In Preuß (2018), (4.5) and (4.6) were analysed under the assumption of an exact handling of geometries, which implies exact equality between both formulations. (C.f. (Preuß, 2018, lemma 2.1, remark 3)). Furthermore, the paper Myrbäck & Zahedi (2024) contains a computational study on a formulation similar to (4.6) under the assumption of exact geometries and including a different variant for the GP stabilization.
## 5. Stability and continuity analysis

One main ingredient for the numerical analysis of the suggested method is the relevant norms. We define the following discrete norms:

$$ \|u\|_J^2 := \|u\|_c^2 + \|u\|_J^2, \quad \|u\|_c^2 := \|\nabla u\|_{Q^h}^2 + \|u\|_J^2, \quad \|u\|_J^2 := J(u, u) \quad (5.1) $$

$$ \|u\|_J^2 := \|u\|^2 + \|u\|_J^2, \quad \|u\|^2 := \Delta t \|\partial_t^\Theta u\|_{Q^h}^2 + \|u\|_c^2, \quad (5.2) $$

$$ \text{with } \|u\|^2 := \sum_{n=1}^{N-1} \|[u]^n\|_{\Omega^h(t_n)}^2 + (u_+^0, u_+^0)_{\Omega^h(0)} + (u_-^N, u_-^N)_{\Omega^h(T)}, $$

where $\|\cdot\|_J$ and $\|\cdot\|$ correspond to the GP and upwind treatment, respectively.

For continuity another norm is used to bound the pairing of $\partial_t u$ and $v$, as well as the pairing of $[u]^n$ and $v$ on time slice boundaries. It contains an $L^2$-norm that is scaled reciprocal to the time derivative part in $\|\cdot\|$ and full control on the time step traces. We define $\|u\|_{J,*}^2 := \|u\|_*^2 + \|u\|_J^2$ with

$$ \|u\|_*^2 := \Delta t^{-1} \|u\|_{Q^h}^2 + \|\nabla u\|_{Q^h}^2 + \|u\|_*^2, \quad \text{with } \|u\|_*^2 := \sum_{n=1}^{N} \|u^n\|_{\Omega^h(t_n)}^2. \quad (5.3) $$

For the subsequent analysis, the norms $\|\cdot\|$, respectively, $\|\cdot\|_J$ and $\|\cdot\|_*$, respectively are the most important $\|\cdot\|_{J,*}$, $\|\cdot\|_c$ and $\|\cdot\|_J$ represent weaker norms that occur naturally in some of our subresults, as well as in the literature. Let us give a brief outline of the stability analysis. Fundamentally, we obtain inf-sup stability by considering for some $u \in \tilde{W}_h$ a linear combination $\tilde{v}(u) = \alpha u + \Pi^{\perp}(\Delta t \partial_t^{\Theta} u)$ for $\alpha \in \mathbb{R}$, where $\Pi^{\perp}$ is an operator subtracting the mean value at the end of each time slice to ensure that $\Pi^{\perp}v \in \tilde{W}_h$ for $v \in W_h$. Note that $B_h(u, u)$ is close to $B_{mc}(u, u)$, cf. Lemma 4.1, so that with (4.1) $B_h(u, u)$ yields important contributions to control the norm $\|\cdot\|_c$ up to geometry approximation errors. Similarly, the GP terms yield control on the GP norm. The more important and technically more involved part is to control the time derivative norm parts. These parts are then obtained from $B_h(u, \Pi^{\perp}(\Delta t \partial_t^{\Theta} u))$. When controlling the time derivative part with $B_h(u, \Pi^{\perp}(\Delta t \partial_t^{\Theta} u))$ possibly negative terms may occur, which are then absorbed by the contributions that are controlled by the symmetric contributions and, hence, choosing $\alpha$ sufficiently large.

This section is organized as follows: firstly, Section 5.1 outlines the fundamental assumptions regarding velocity, the space-time geometry and space-time anisotropy. In Section 5.2 we present a specific trace inequality that allows to bound terms arising from geometry approximation errors. This estimate facilitates control over the error term from Lemma 4.1 and provides the initial component of the stability estimate in lemma 5.2, indicating the necessity of controlling the time derivative component. To prepare for this we explore several results on GP stabilization in Section 5.3. These results are somewhat unconventional and more technical than usual due to the space-time setting and the parametric space-time mapping involved in the discretization method. After analysing the splitting of integral constants and remainder as in Section 4.2 we are ready to control the time derivative in Section 5.4, at the cost of some terms that are however controlled by symmetric contributions. Finally, Section 5.5 synthesizes the previously established results to derive the main stability and continuity estimates.

**REMARK 5.1** (Related stability analyses in the literature). In Lehrenfeld & Reusken (2013) a two-domain convection-diffusion problem with an unfitted moving interface is considered. The interface condition
is imposed by a Nitsche method that in combination with a special weighting allows for a stable discretization without GP. In the error analysis a weaker norm, similar to $||\cdot||_c$, is used that can be controlled solely by coercivity arguments as in (4.1). Hence, it does not yield control of the time derivate parts and further requires a stronger continuity norm. Compared with the analysis presented here the analysis in Lehrenfeld & Reusken (2013) leads to a bound for the error in a weaker norm in terms of an approximation error in a stronger norm. The analysis in Badia et al. (2023) proceeds similarly, but is adapted to their setting where stability is achieved through element aggregation. Neither analysis considers geometry approximation errors. In Remark 6.4 we will compare our final error bound with those obtained in Badia et al. (2023) and Lehrenfeld & Reusken (2013) again. With regard to stability our analysis is closer to the one in Larson et al. (2024). Both analyses exploit a linear combination of the discrete function and its time derivative to obtain stability in a norm that contains a weighted time derivative, cp. (Larson et al., 2024, lemma 4.5) and Theorem 5.1. This idea already appears in Schieweck (2010) in a different context.

## 5.1 General assumptions

We introduce some mild assumptions concerning the relation between time step and mesh width.

**ASSUMPTION 5.1.** Time step $\Delta t$ and mesh size $h$ satisfy $h^2 \lesssim \Delta t$.

**ASSUMPTION 5.2.** Time step $\Delta t$ and mesh size $h$ satisfy $\Delta t^{q_t+1} \lesssim h^{1+\frac{1}{2}}$.

Moreover, we pose the following assumption on a ratio between boundary and (delayed) volume measure within each time slab.

**ASSUMPTION 5.3.** There are $C_\Omega$, $C_Q > 0$, s.t. for sufficiently small $h$ and $\Delta t$, cf. Assumption 3.2,

$$ \frac{1}{\Delta t} \max_{n=1,...,N} \iint_{\Omega} \frac{|\partial \Omega^h(t)|}{|\Omega^h(t_n)|} dt \le C_\Omega \quad \text{and} \quad \frac{1}{\Delta t} \max_{n=1,...,N} \frac{|Q^{h,n}|}{|\Omega^h(t_n)|} \le C_Q. \qquad (5.4) $$

Furthermore, we assume that $|\Omega(t_n)|$ and, hence, $|\Omega^h(t_n)|$ is bounded from above and below uniformly in $n = 1, ..., N$, i.e., $|\Omega(t_n)| \simeq |\Omega^h(t_n)| \simeq 1$.

## 5.2 Controlling boundary integrals and control by symmetric testing

We require specific special trace inequalities to bound the terms arising from geometry approximation errors as in Lemma 4.1.

**LEMMA 5.1 (Special trace inequality).** For all $u \in \tilde{W}_h$ it holds

$$ \iint_{\Omega} \|u(\cdot, t)\|_{\partial \Omega^h(t)}^2 dt \lesssim \|u\|_{\partial_s Q^h}^2 \lesssim \|\nabla u\|_{Q^h}^2 + \Delta t^2 \|\partial_t^\Theta u\|_{Q^h}^2 \\ \le C_{5,1} \min\{\|u\|^2, \|u\|_*^2\}, \qquad (5.5a) $$

$$ \|u\|_{Q^h} \lesssim \min\{\|u\|, \|u\|_*\}. \qquad (5.5b) $$

$$ \text{Moreover, for } u \in H^1(Q^h), \quad \|u\|_{\partial_s Q^h} \lesssim \|u\|_* \qquad (5.5c) $$
**Proof.** We start with (5.5a) and the bound for the $||\cdot||$-norm. We consider an arbitrary fixed time interval $I_n = (t_{n-1}, t_n]$ and aim at estimating $\int_{I_n} ||u(\cdot, t)||^2_{\partial\Omega^{h}(t)} dt$ for $u \in \tilde{W}_h$. Let $\theta_t$ be the affine mapping that maps the unit interval $\hat{I} := (0, 1]$ onto $I_n$ and $\Theta_t(x, t) = (x, \theta_t(t))$ the corresponding space-time mapping. Further, let $Q_i^{h,n} = \Theta_t^{-1}(Q^{h,n})$ be the space-time domain $Q^{h,n}$ rescaled in time. We obtain the desired result after the application of 1. the transformation formula (in time), 2. boundedness of the geometry motion, 3. a trace inequality in space-time, 4. a version of the Poincaré inequality (with boundary value constraints) for which we exploit $u \in \tilde{W}_h$, 5. the chain rule, 6. Fubini's theorem and transformation formula in time and 7. a triangle inequality based on (3.14) and the boundedness of $\partial_t \Theta_h^{st}$, cf. Equation (3.6d) for $(l_t, l_s) = (1, 0)$:

$$ 
\begin{aligned}
\Delta t^{-1} \int_{I_n} ||u(\cdot, t)||^2_{\partial\Omega^{h}(t)} dt &\stackrel{1}{=} \int_{I} ||u(\cdot, \theta_t(t))||^2_{\partial\Omega^{h}(\theta_t(t))} dt \\
&\stackrel{2}{\lesssim} ||u \circ \Theta_t||^2_{\partial_s Q_i^{h,n}} \\
&\stackrel{3}{\lesssim} ||u \circ \Theta_t||^2_{H^1(Q_i^{h,n})} \\
&\stackrel{4}{\lesssim} ||\nabla(u \circ \Theta_t)||^2_{Q_i^{h,n}} + ||\partial_t(u \circ \Theta_t)||^2_{Q_i^{h,n}} \\
&\stackrel{5}{\lesssim} ||(\nabla u) \circ \Theta_t||^2_{Q_i^{h,n}} + \Delta t^2 ||(\partial_t u) \circ \Theta_t||^2_{Q_i^{h,n}} \\
&\stackrel{6}{\lesssim} \Delta t^{-1} ||\nabla u||^2_{Q_i^{h,n}} + \Delta t ||\partial_t u||^2_{Q_i^{h,n}} \\
&\stackrel{7}{\lesssim} \Delta t^{-1} ||\nabla u||^2_{Q_i^{h,n}} + \Delta t ||\partial_t^{\Theta} u||^2_{Q_i^{h,n}}.
\end{aligned}
$$

The same arguments have been discussed with more detail in (Wendler, 2022, lemma 2.13) and Heimann (2025). Now, we can bound $\Delta t^2 ||\partial_t^{\Theta} u||^2_{Q_i^{h,n}}$ obviously with $\Delta t ||\partial_t^{\Theta} u||^2_{Q_i^{h,n}}$ to see the $||\cdot||$-norm bound. For the $||\cdot||_*$-norm we bound the time derivative norm with an $L^2$-norm in space-time using corollary 3.1 and $\hat{u} = u \circ \Phi^{st}$

$$ \Delta t^2 ||\partial_t^{\Theta} u||^2_{Q_i^{h,n}} \simeq \Delta t^2 ||\partial_t \hat{u}||^2_{Q_i^{\text{lin},n}} \simeq ||\hat{u}||^2_{Q_i^{\text{lin},n}} \simeq ||u||^2_{Q_i^{h,n}} \lesssim \Delta t^{-1} ||u||^2_{Q_i^{h,n}}. $$

(5.5b) follows from similar arguments (put $||u \circ \Theta_t||^2_{Q_i^{h,n}}$ in front of 3. above).

In relation to (5.5c) we argue similarly: note that a merely spatial trace theorem on Sobolev spaces allows to bound $||u||_{\partial\Omega^{h}(t)}$ by $||u||_{H^1(\Omega^{h}(t))}$. Integrating over time allows to establish a bound by $||u||_*$, as it contains the relevant summands. $\square$

We can now state the first result for a contribution of the stability analysis.

**LEMMA 5.2** For all $u \in W_h$ there holds

$$ B_h(u, u) \geq \frac{1}{2} \|u\|_c^2 - \frac{C_{4,1}}{2} (h^{q_s} + \Delta t^{q_t}) \cdot \|u\|^2. \quad (5.6) $$

**Proof.** With Lemma 4.1 and (5.5a) we find

$$ B_h(u, u) \geq \frac{1}{2} (B_h(u, u) + B_{mc}(u, u)) - \frac{C_{4,1}}{2} (h^{q_s} + \Delta t^{q_t}) \cdot \|u\|^2. \quad (5.7) $$

Applying (4.1) yields the claim. $\square$
We observe that, although the perturbation is of higher order, the geometry approximation error cannot be controlled by the $||\cdot||_c$ norm. To obtain control on the stronger norm $||\cdot||$ we need to control the time derivative part for which we require the GP stabilization.

**REMARK 5.2.** To circumvent the difficulties stemming from the nonpositive term in (5.7) one could consider another variant of the method based on the bilinear form $B'_h = \frac{1}{2}B_h + \frac{1}{2}B_{mc}$ for which we—by construction—obtain $B'_h(u, u) \ge \frac{1}{2} \|u\|_c^2$. This would allow to follow a coercivity-based analysis w.r.t. the norm $||\cdot||_c$ even with inexact geometry handling and without GP stabilization.

## 5.3 GP stabilization on deformed space-time meshes

The analysis of the GP stabilization relies on a few technical assumptions, cf. also Preuß (2018). The general idea (see also Burman (2010) and Burman et al. (2015)) is that in order to yield control on cut elements with potentially small cut parts of the discrete domain neighbouring elements with sufficient support in the discrete domain are used. The coordination between interior and cut elements in the space-time setting is the purpose of the following exterior-interior mapping between elements.

**ASSUMPTION 5.4** (Existence of an exterior-interior mapping for the GP). There exists a mapping $\mathcal{B}: \Omega_{\mathcal{E}}^{\text{lin},I_n} \to \Omega_{\mathcal{J}}^{\text{lin},I_n}$ such that

1. The number of exterior elements mapped onto $T \in \Omega_{\mathcal{J}}^{\text{lin},I_n}$ is bounded

$$ \#\{\mathcal{B}^{-1}(T)\} \le C \quad \text{for } T \in \Omega_{\mathcal{J}}^{\text{lin},I_n}. \qquad (5.8a) $$

2. For $T \in \Omega_{\mathcal{E}}^{\text{lin},I_n} \setminus \Omega_{\mathcal{J}}^{\text{lin},I_n}$ there exists a path $\{T_i\}_{i=0}^M$ between $T_0 = T$ and $T_M = \mathcal{B}(T) \in \Omega_{\mathcal{J}}^{\text{lin},I_n}$ across the set of facets $\mathcal{F}_T^n$ such that $\mathcal{F}_T^n \subseteq \mathcal{F}_R^n$.

3. Moreover, the maximum occurrences of a facet in the paths of $\mathcal{B}$, which we denote by $N_F := \#\{T \in \Omega_{\mathcal{E}}^{\text{lin},I_n} \setminus \Omega_{\mathcal{J}}^{\text{lin},I_n} | F \in \mathcal{F}_T^n\}$, is bounded:

$$ \max_{F \in \mathcal{F}_R^n} N_F \le C_F (1 + \Delta t/h), \qquad (5.8b) $$

where $C_F$ is independent of $\Delta t$ and $h$.

The paths referred to in point 2, which give rise to the mapping $\mathcal{B}$, are illustrated in Fig. 4. While it is reasonable to assume that the number of elements mapped back to a specific element $T \in \Omega_{\mathcal{J}}^{\text{lin},I_n}$ is bounded uniformly the same is not true for the maximum number of occurrences of a facet in the paths. This is essentially due to the fact that the thickness of the number of cut elements depends on the ratio of $\Delta t$ and $h$. To account for this behaviour we introduced the scaling in (5.8b), which further motivates the scaling in $\tilde{\gamma}_j$. For additional discussion of Assumption 5.4 we refer to (Preuß, 2018, remark 4).

The paths between full inner elements and cut elements alluded to in Assumption 5.4 are important for the stability proof. The mechanism of moving along facet elements is explicated in the following lemma:

**LEMMA 5.3** Let $T_1$ and $T_2$ be the elements of the facet patch $\omega_F = T_1 \cup T_2$. Let $v = \hat{v} \circ (\Theta_h(\cdot, t))^{-1}$, where $\hat{v}$ is a piecewise polynomial of degree at most $k_s$ defined on $\omega_F$, and $t$ is a fixed time, $t \in I_n$, in a time interval $I_n$. Then, there holds

$$ \|D^{l_s} v\|_{\Theta_h(T_1,t)}^2 \lesssim \|D^{l_s} v\|_{\Theta_h(T_2,t)}^2 + h^{-2l_s} \|[v]_{\omega_F^h}\|_{\Theta_h(T_1,t)}^2, \quad l_s \in \{0,1\}. \qquad (5.9a) $$
<figure><img src="image_6.png" /><figcaption>FIG. 4. Illustration of examples of paths from cut elements to the interior as required by Assumption 5.4. C.f. Fig. 2 for a sketch of the same geometric setting with a highlight on $\Omega_{\varepsilon}^{\text{lin},I_n}$.</figcaption></figure>

*Proof.* Can be found in Appendix A.1. □

This result can be applied in combination with the interior-exterior mapping from Assumption 5.4 in order to derive a result that allows to bound norms on the whole active domain by corresponding norms on the interior domain plus GP contributions. This will be relevant later for the stability analysis:

LEMMA 5.4 It holds for every $u \in W_h, \tau \in I_n$

$$ \|D^l u\|_{Q_{\varepsilon}^{h,n}}^2 \lesssim \|D^l u\|_{Q_{\varepsilon}^{h,n}}^2 + \frac{h^2-2l_s}{\gamma_J} j_h^n(u, u), \quad l_s \in \{0, 1\}. \quad (5.10) $$

*Proof.* Can be found in Appendix A.2. □

In the next lemma we will derive a bound for the norm of a discrete function on the temporal parts of the space-time boundary against the norm on the space-time volume:

LEMMA 5.5 For any $u \in W_h$ it holds

$$ \|u_{+}^{n-1}\|_{\Omega^h(t_{n-1})}^2 + \|u_{-}^{n}\|_{\Omega^h(t_n)}^2 \lesssim \frac{1}{\Delta t} \|u\|_{Q_{\varepsilon}^{h,n}}^2 \le \frac{C_{5,5}}{\Delta t} \left( \|u\|_{Q_{\varepsilon}^{h,n}}^2 + \frac{h^2}{\gamma_J} j_h^n(u, u) \right). \quad (5.11) $$

*Proof.* Let $\hat{u} = u \circ (\Theta_h^{\text{st}})^{-1}$ be the discrete function on the undeformed mesh. Then, at each point $\hat{x} \in T \in \Omega_{\varepsilon}^{\text{lin},I_n}$ the function $\hat{u}(\cdot, t)$ is polynomial and we can apply an inverse trace inequality again, cf. (Warburton & Hesthaven, 2003, theorem 2), yielding

$$ \|u_{+}^{n-1}\|_{\Omega^h(t_{n-1})}^2 \lesssim \|\hat{u}_{+}^{n-1}\|_{\Omega_{\varepsilon}^{\text{lin},I_n}}^2 \lesssim \frac{1}{\Delta t} \|\hat{u}\|_{Q_{\varepsilon}^{\text{lin},n}}^2 \lesssim \frac{1}{\Delta t} \|u\|_{Q_{\varepsilon}^{h,n}}^2, $$
where we further exploited the equivalence between the norms on $Q^{\text{lin}}$ and $Q^h$. In combination with Lemma 5.4 we obtain the result for $u_+^{n-1}$. The proof for $u^n$ follows the same lines. $\square$

This finishes the estimates for discrete functions $u$ itself. Next, we will also consider a similar estimate for the time derivative of the discrete function. Both will be of importance for later application in the stability proof.

LEMMA 5.6 For any $u \in W_h$ it holds

$$ \Delta t \|D^l \partial_t^\Theta u\|_{Q_h^{h,n}}^2 \lesssim \frac{C_{5,6}}{\Delta t} \left( \|D^l u\|_{Q_h^{h,n}}^2 + \frac{h^{2-2l_s}}{\gamma_J} j_h^n(u, u) \right), \quad l_s \in \{0, 1\}. \quad (5.12) $$

*Proof.* First, let us consider $l_s = 0$. Fix $u \in W_h$. Then, by Corollary 3.1 for $\hat{u} = u \circ (\Theta_h^{\text{st}})$, it holds $\Delta t \|\partial_t^\Theta u\|_{Q_h^{h,n}}^2 \lesssim \Delta t \|\partial_t \hat{u}\|_{Q_h^{\text{lin},n}}^2$. For $\hat{u}$, we again apply inverse estimates in time, see e.g. (Brenner & Scott, 2008, theorem 4.5.11), yielding

$$ \Delta t \|\partial_t^\Theta u\|_{Q_h^{h,n}}^2 \lesssim \frac{1}{\Delta t} \|\hat{u}\|_{Q_h^{\text{lin},n}}^2 \lesssim \frac{1}{\Delta t} \|u\|_{Q_h^{h,n}}^2. $$

Applying Lemma 5.4 yields the result for $u \in W_h$. Next, for $l_s = 1$ we observe that because $\partial_t$ and $\nabla$ commute on the reference configuration

$$ \|\nabla \partial_t^\Theta u\|_{Q_h^{h,n}}^2 \stackrel{(3.16)}{=} \|\nabla ((\partial_t \hat{u}) \circ (\Theta_h^{\text{st}})^{-1})\|_{Q_h^{h,n}}^2 \stackrel{(3.7)}{\simeq} \|\nabla \partial_t \hat{u}\|_{Q_h^{\text{lin},n}}^2 = \|\partial_t \nabla \hat{u}\|_{Q_h^{\text{lin},n}}^2. $$

Now, exploiting that $\nabla \hat{u}$ is an elementwise polynomial again to enable inverse inequalities again yields

$$ \Delta t \|\nabla \partial_t^\Theta u\|_{Q_h^{h,n}}^2 \lesssim \Delta t \|\partial_t \nabla \hat{u}\|_{Q_h^{\text{lin},n}}^2 \lesssim \frac{1}{\Delta t} \|\nabla \hat{u}\|_{Q_h^{\text{lin},n}}^2 \lesssim \frac{1}{\Delta t} \|\nabla u\|_{Q_h^{h,n}}^2. $$

Applying Lemma 5.4 again concludes the proof. $\square$

Finally, we show that a relation between the GP of the time derivative and the function itself can be established:

LEMMA 5.7 There holds for all $u \in W_h$

$$ \Delta t^2 \|\partial_t^\Theta u\|_J^2 \le C_{5,7} \left( \|u\|_J^2 + \gamma_J (h + \Delta t) \|\nabla u\|_{Q_h^h}^2 \right). \quad (5.13) $$

*Proof.* We fix a time interval $I_n$ and a facet $F$ from $\mathcal{F}_R^n$ and observe

$$ \begin{aligned} j_F^n(\partial_t^\Theta u, \partial_t^\Theta u) &= \iint_{h-1}^{h} \iint_{\partial_F(t)} \frac{\tilde{\gamma}_J}{h^2} [\partial_t^\Theta u]_{\omega_F(t)}^{2h} \, dx \, dt \\ &\stackrel{(4.3)}{=} \iint_{h-1}^{h} \sum_{i,j=1,2, i \neq j} \iint_{\partial_h(T_i, t)} \frac{\tilde{\gamma}_J}{h^2} \left( \partial_t^\Theta u(x) - \mathcal{E}^p(\partial_t \hat{u}|_{T_j}) \left( \left( \mathcal{E}^p \Theta_h(\cdot, t) \right)_{|T_j}^{-1}(x) \right) \right)^2 \, dx \, dt \lesssim \dots \end{aligned} $$
Recalling $\partial_t^\Theta u(x) = (\partial_t \hat{u}|_{T_i}) \circ (\Theta_h(\cdot, t)|_{T_i})^{-1}$, and abbreviating $\mathcal{E}^p \Theta_h(t)|_{T_j}$ as $\Theta_h^{\mathcal{E}}(t)$, $j = 1,2$, we continue

$$ \begin{aligned} & \dots \lesssim \iint_{T_{i-1}}^{T_i} \sum_{i,j=1,2, \atop i \neq j} \iint_i \frac{\tilde{\gamma}_J}{h^2} \left( \partial_t \hat{u}|_{T_i} - \mathcal{E}^p (\partial_t \hat{u}|_{T_j}) \circ \underbrace{(\Theta_{h,T_j}^{\mathcal{E}}(t))^{-1} \circ (\Theta_h(t)|_{T_i})}_{:=\psi} (x) \right)^2 dx \\ & \lesssim \iint_{T_{i-1}}^{T_i} \sum_{i,j=1,2, \atop i \neq j} \iint_i \frac{\tilde{\gamma}_J}{h^2} (\partial_t \hat{u}|_{T_i} - \mathcal{E}^p (\partial_t \hat{u}|_{T_j}))^2 + \frac{\tilde{\gamma}_J}{h^2} (\mathcal{E}^p (\partial_t \hat{u}|_{T_j}) - \mathcal{E}^p (\partial_t \hat{u}|_{T_i}) \circ \psi)^2 dx. \end{aligned} $$

The first expression $\partial_t \hat{u}|_{T_i} - \mathcal{E}^p(\partial_t \hat{u}|_{T_j}) = \partial_t(\hat{u}|_{T_i} - \mathcal{E}^p(\hat{u}|_{T_j}))$ is a polynomial in time and we can again apply an inverse inequality in time. For the second term we can find a bound as through proximity of $\psi$ to the identity:

$$ \mathcal{E}^p(\partial_t \hat{u}|_{T_j}) - \mathcal{E}^p(\partial_t \hat{u}|_{T_i}) \circ \psi \lesssim \|\nabla \partial_t \mathcal{E}^p(\hat{u}|_{T_j})\|_{\infty, T_i} \cdot \underbrace{\|\mathrm{id} - \psi\|_{\infty, T_i}}_{\lesssim h^2 + \Delta t^{q_i+1}} \text{ by (3.6d)} $$

The term $\nabla \partial_t \mathcal{E}^p(\hat{u}|_{T_j})$ is also a polynomial so that applying inverse inequalities gives

$$ \begin{aligned} j_F^n (\partial_t^\Theta u, \partial_t^\Theta u) \lesssim & \iint_{T_{i-1}}^{T_i} \sum_{i,j=1,2, \atop i \neq j} \iint_i \frac{\tilde{\gamma}_J}{h^2 \cdot \Delta t^2} (\hat{u}|_{T_i} - \mathcal{E}^p (\hat{u}|_{T_j}))^2 dx \\ & + \sum_{i,j=1,2, \atop i \neq j} \frac{\tilde{\gamma}_J \cdot (h^2 + \Delta t^{q_i+1})^2}{h^2 \Delta t^2} |T_i| \left( \|\nabla \mathcal{E}^p (\hat{u}|_{T_j})\|_{\infty, T_i \times I_n} \right)^2 dx. \end{aligned} $$

To bring the expression $(\hat{u}|_{T_i} - \mathcal{E}^p(\hat{u}|_{T_j}))^2$ back into the form of the GP jump we apply a triangle inequality to get

$$ (\hat{u}|_{T_i} - \mathcal{E}^p(\hat{u}|_{T_j}))^2 \lesssim (\hat{u}|_{T_i} - \mathcal{E}^p(\hat{u}|_{T_j}) \circ \psi)^2 + \underbrace{(\mathcal{E}^p(\hat{u}|_{T_j}) \circ \psi - \mathcal{E}^p(\hat{u}|_{T_j}))^2}_{\lesssim (h^2 + \Delta t^{q_i+1})^2 \|\nabla \mathcal{E}^p(\hat{u}|_{T_j})\|_{\infty, T_i}^2} $$

With Assumption 5.2 we can bound $(h^2 + \Delta t^{q_i+1})^2/h^2 \lesssim h$, which together with the last two results yields

$$ \begin{aligned} & \Delta t^2 j_F^n (\partial_t^\Theta u, \partial_t^\Theta u) \lesssim j_F^n (u, u) + \iint_{T_{i-1}}^{T_i} \sum_{i,j=1,2, \atop i \neq j} \tilde{\gamma}_J h |T_i| (\|\nabla \mathcal{E}^p (\hat{u}|_{T_j})\|_{\infty, T_i})^2 \\ & \lesssim j_F^n (u, u) + \iint_{T_{i-1}}^{T_i} \tilde{\gamma}_J h \|\nabla u\|_{\omega_F^h(t)}^2 \lesssim j_F^n (u, u) + \iint_{T_{i-1}}^{T_i} \gamma_J (h + \Delta t) \|\nabla u\|_{\omega_F^h(t)}^2. \end{aligned} $$
Summing over all facet contributions in the total GP term and applying applying Lemma 5.4 another time then yields

$$ \Delta t^2 \|\partial_t^\Theta u\|_J^2 \lesssim \|u\|_J^2 + \gamma_J(h + \Delta t) \|\nabla u\|_{Q^h}^2. \quad \square $$

## 5.4 Controlling the time derivative

In the stability proof it will be relevant when comparing discrete norms to apply an operator that subtracts the final integral value on the respective time slice, which we will denote as $\tilde{\Pi}(\cdot)$ from a discrete function. This construction follows Wendler (2022). We call this operator $\Pi^\perp$ and define it as follows.

**DEFINITION 5.1.** We define $\tilde{\Pi}, \Pi^\perp: L^2(Q^h) \to L^2(Q^h)$ as

$$ \tilde{\Pi}(u) := \sum_{n=1}^{N} \int_{\Omega^h(t_n)} u_n^- dx \cdot 1_{I_n} \quad \text{and} \quad \Pi^\perp(u) := u - \tilde{\Pi}(u), \quad u \in L^2(Q^h), \qquad (5.14) $$

where $1_{I_n}$ denotes the time interval indicator function and $f_S$ denotes the mean value operation, cf. Section 4.2.

By construction, there holds $\Pi^\perp(u) \in \tilde{W}_h$ for all $u \in W_h$. In the stability proof we will exploit the compatibility of the $\|\cdot\|_J$ norms between a function $u$ and a specially tailored function $\tilde{v}(u)$, which contains a summand of $\Pi^\perp(\Delta t \partial_t^\Theta u)$. To prepare this we show the following:

**LEMMA 5.8** For all $u \in W_h$ there holds for sufficiently small mesh sizes $h$ and time steps $\Delta t$ that

$$ \|\Pi^\perp(\Delta t \partial_t^\Theta u)\|_J \lesssim \|u\|_J. \qquad (5.15) $$

*Proof.* We start by writing out the squared left-hand side and exploiting that slice-wise constants vanish for all norm contributions except for the time jump terms:

$$ \|\Pi^\perp(\Delta t \partial_t^\Theta u)\|_J^2 = \Delta t^3 \|(\partial_t^\Theta)^2 u\|_{Q^h}^2 + \Delta t^2 \|\nabla \partial_t^\Theta u\|_{Q^h}^2 + \Delta t^2 \|\partial_t^\Theta u\|_J^2 + \|\Pi^\perp(\Delta t \partial_t^\Theta u)\|^2. $$

We bound the first three terms:

$$ \begin{align*} \Delta t^3 \|(\partial_t^\Theta)^2 u\|_{Q^h}^2 &\lesssim \overbrace{h^2 \Delta t}^{(5.12)} \gamma_J^{-1} \|\partial_t^\Theta u\|_J^2 + \Delta t \|\partial_t^\Theta u\|_{Q^h}^2, \\ &\lesssim \gamma_J^{-1} \|u\|_J^2 + (h + \Delta t) \|\nabla u\|_{Q^h}^2 + \Delta t \|\partial_t^\Theta u\|_{Q^h}^2 \lesssim \|u\|_J^2, \\ \Delta t^2 \|\nabla \partial_t^\Theta u\|_{Q^h}^2 &\lesssim \overbrace{\gamma_J^{-1} \|u\|_J^2 + \|\nabla u\|_{Q^h}^2}^{(5.12)} \lesssim \|u\|_J^2, \\ \Delta t^2 \|\partial_t^\Theta u\|_J^2 &\lesssim \overbrace{\|u\|_J^2 + \gamma_J(h + \Delta t) \|\nabla u\|_{Q^h}^2}^{(5.13)} \lesssim \|u\|_J^2. \end{align*} $$
It only remains to bound the jump norm summand, which we split into $\Delta t \partial_t^\Theta u$ and the piecewise constant term:

$$ \|\Pi^{\perp}(\Delta t \partial_t^{\Theta} u)\|^2 \lesssim \|\Delta t \partial_t^{\Theta} u\|^2 + \left\|\sum_{n=1}^{N} \int_{\Omega^h(t_n)} (\Delta t \partial_t^{\Theta} u)_-^n \, dx \cdot 1_{I_n}\right\|^2. $$

With $\|\Delta t \partial_t^{\Theta} u\|^2 = \Delta t^2 \|\partial_t^{\Theta} u\|^2$ and recalling $\partial_t^{\Theta} u \in W_h$ (cf. (3.1)) we get

$$ \begin{aligned} \|\Delta t^2 \|\partial_t^{\Theta} u\|^2 &\lesssim \Delta t^2 \sum_{n=1}^{N} \left( \| (\partial_t^{\Theta} u)_+^{n-1} \|_{\Omega^h(t_{n-1})}^2 + \| (\partial_t^{\Theta} u)_-^n \|_{\Omega^h(t_n)}^2 \right) \\ &\lesssim \Delta t \|\partial_t^{\Theta} u\|_{Q^h}^2 + \gamma_j^{-1} \underbrace{h^2 \Delta t \|\partial_t^{\Theta} u\|_{J}^2}_{\text{(5.11)}} \lesssim \|u\|_{J}^2. \\ &\lesssim \Delta t^2 \text{ by Ass. 5.1} \end{aligned} $$

Similarly,

$$ \begin{aligned} \left\|\left[\sum_{n=1}^{N} \int_{\Omega^h(t_n)} (\Delta t \partial_t^{\Theta} u)_-^n \, dx \cdot 1_{I_n}\right]\right\|^2 &\lesssim \frac{1}{\Delta t} \left\|\sum_{n=1}^{N} \int_{\Omega^h(t_n)} (\Delta t \partial_t^{\Theta} u)_-^n \, dx \cdot 1_{I_n}\right\|^2_{Q^h} \\ &\lesssim \frac{1}{\Delta t} \sum_{n=1}^{N} \| (\Delta t \partial_t^{\Theta} u)_-^n \|_{\Omega^h(t_n)}^2 \cdot |Q^{hn}|/|\Omega^h(t_n)| \stackrel{\text{Ass. 5.3}}{\lesssim} \sum_{n=1}^{N} \| (\Delta t \partial_t^{\Theta} u)_-^n \|_{\Omega^h(t_n)}^2 \\ &\lesssim \gamma_j^{-1} h^2 \Delta t \|\partial_t^{\Theta} u\|_{J}^2 + \Delta t \|\partial_t^{\Theta} u\|_{Q^h}^2 \lesssim \|u\|_{J}^2, \end{aligned} $$

where in the first estimate we exploited $\|1_{I_n}\|_J = 0$. Taking all the results together yields the claim. □

In the stability proof we will consider a linear combination between $u$ and $\Pi^{\perp}(\Delta t \partial_t^{\Theta} u)$ as the test function $v$. The symmetric contribution $B_h(u, u)$ has been estimated in Lemma 5.2 already. Now, we also want to gain control of the time derivative with $v = \Pi^{\perp}(\Delta t \partial_t^{\Theta} u)$. This is the aim of the following lemma:

LEMMA 5.9 There exists a constant $C_{5,9}$, independent of $h$ and $\Delta t$, such that for any $u \in \tilde{W}_h$ there holds

$$ (B_h + J)(u, \Pi^{\perp}(\Delta t \partial_t^{\Theta} u)) \geq \frac{\Delta t}{2} \|\partial_t^{\Theta} u\|_{Q^h}^2 - C_{5,9} \|u\|_J^2. $$

*Proof.* We define $U := \bar{\Pi}(\Delta t \partial_t^{\Theta} u)$ and with $J(u, U) = 0$ we see

$$ (B_h + J)(u, \Pi^{\perp}(\Delta t \partial_t^{\Theta} u)) = \underbrace{B_h(u, \Delta t \partial_t^{\Theta} u)}_{①} + \underbrace{J(u, \Delta t \partial_t^{\Theta} u)}_{②} - \underbrace{B_h(u, U)}_{③}. $$

We will bound the three r.h.s. terms from below one after the other.
Term ①: $B_h(u, \Delta t \partial_t^\Theta u)$ We start by writing out $B_h(u, \Delta t \partial_t^\Theta u)$:

$$B_h(u, \Delta t \partial_t^\Theta u) = (\partial_t u, \Delta t \partial_t^\Theta u)_{Q^h} + (\tilde{w} \cdot \nabla u, \Delta t \partial_t^\Theta u)_{Q^h} + (\nabla u, \Delta t \nabla \partial_t^\Theta u)_{Q^h} \\ + \sum_{n=1}^{N-1} \Delta t ([u]^n, (\partial_t^\Theta u)_+^n)_{\Omega^h(t_n)} + \Delta t (u_+^0, (\partial_t^\Theta u)_+^0)_{\Omega^h(0)}.$$

From Definition 3.1 we obtain $\partial_t u = \partial_t^\Theta u - \underbrace{(\partial_t \Theta_h^{st} \circ (\Theta_h^{st})^{-1})}_{=:\tilde{w}} \cdot \nabla u$, leading to

$$B_h(u, \Delta t \partial_t^\Theta u) - \Delta t \|\partial_t^\Theta u\|_{Q^h}^2 = \underbrace{(\tilde{w} - \tilde{w}) \cdot \nabla u, \Delta t \partial_t^\Theta u)_{Q^h}}_{\textcircled{a}} \\ + \underbrace{\sum_{n=1}^{N-1} \Delta t([u]^n, (\partial_t^\Theta u)_+^n)_{\Omega^h(t_n)} + \Delta t(u_+^0, (\partial_t^\Theta u)_+^0)_{\Omega^h(0)}}_{\textcircled{b}} + \underbrace{(\nabla u, \Delta t \nabla \partial_t^\Theta u)_{Q^h}}_{\textcircled{c}}.$$

All those terms will now be estimated separately. In general, we apply a combination of Cauchy–Schwarz and Youngs inequality, combined with several previous results. We start with ①:

$$\begin{aligned} \textcircled{a} &= \sum_{n=1}^{N} ((\tilde{w} - \tilde{w}) \cdot \nabla u, \Delta t \partial_t^\Theta u)_{Q^h} \le 2\Delta t \|(\tilde{w} - \tilde{w}) \cdot \nabla u\|_{Q^h}^2 + \frac{\Delta t}{8} \|\partial_t^\Theta u\|_{Q^h}^2 \\ &\le 2\Delta t \|(\tilde{w} - \tilde{w})\|_{\infty}^2 \|\nabla u\|_{Q^h}^2 + \frac{\Delta t}{8} \|\partial_t^\Theta u\|_{Q^h}^2 \le C_{\textcircled{b}} \Delta t \|u\|_{c}^2 + \frac{\Delta t}{8} \|\partial_t^\Theta u\|_{Q^h}^2, \end{aligned}$$

where we set $C_{\textcircled{b}} = 2 \|(\tilde{w} - \tilde{w})\|_{\infty}^2$, which is bounded uniformly in $h$ and $\Delta t$.

Next, we continue with ⑥ applying Cauchy–Schwarz and Young’s inequality:

$$\textcircled{b} \le 2C_{5,5} \|u\|^2 + \frac{\Delta t^2}{8C_{5,5}} \sum_{n=1}^{N} \| (\partial_t^\Theta u)_+^{n-1} \|_{\Omega^h(t_{n-1})}^2.$$

For the sum we exploit $\partial_t^\Theta u \in W_h$ and apply lemma 5.5 and lemma 5.7 where we make use of Assumption 5.1 to bound $h^2 \lesssim \Delta t$. To keep the reference to this implicit constant, which has the role of an inverse CFL condition, let us denote by $C_{iCFL}$ the constant such that $h^2 \le C_{iCFL} \Delta t$. With that,

$$\begin{aligned} \Delta t^2 \sum_{n=1}^{N} \| (\partial_t^\Theta u)_+^{n-1} \|_{\Omega^h(t_{n-1})}^2 &\stackrel{(5.11)}{\le} C_{5,5} \cdot \Delta t \left( \| \partial_t^\Theta u \|_{Q^h}^2 + \frac{h^2}{\gamma_J} \| \partial_t^\Theta u \|_{J}^2 \right) \\ &\stackrel{(5.13)}{\le} C_{5,5} \left( C_{5,7} C_{iCFL} \left\{ \gamma_J^{-1} \| u \|_{J}^2 + (h + \Delta t) \| \nabla u \|_{Q^h}^2 \right\} + \Delta t \| \partial_t^\Theta u \|_{Q^h}^2 \right), \end{aligned} \quad (5.16)$$
which yields with $C_⊕ := \max\{2C_{5,5}, \frac{1}{8}C_{5,7}C_{iCFL} \cdot \max\{\gamma_J^{-1}, h + \Delta t\}$, which is bounded uniformly in $h$ and $\Delta t$, the estimate for the term

$$ \Rightarrow \oplus \le C_⊕ \|u\|_j^2 + \frac{\Delta t}{8} \|\partial_t^\Theta u\|_{Q^h}^2. \quad (5.17) $$

We estimate $\oplus$

applying Cauchy–Schwarz and Young’s inequality with $\varepsilon = \sqrt{C_{5,6}}$:

$$ \odot = (\nabla u, \Delta t \nabla \partial_t^\Theta u)_{Q^h} \stackrel{(5.12)}{\le} \frac{\varepsilon}{2} \|\nabla u\|_{Q^h}^2 + \frac{C_{5,6}}{2\varepsilon} \left( \|\nabla u\|_{Q^h}^2 + \frac{1}{\gamma_J} \|u\|_J^2 \right) \le C_⊕ \|u\|_j^2 $$

with $C_⊕ = \sqrt{C_{5,6}} \cdot \max\{1, \frac{1}{2\gamma_J}\}$ independent of $h$ and $\Delta t$. Putting the estimates for $\oplus, \ominus$ and $\odot$ together we obtain

$$ B_h(u, \Delta t \partial_t^\Theta u) \ge \frac{3}{4} \Delta t \|\partial_t^\Theta u\|_{Q^h}^2 - C_\oplus \|u\|_j^2 $$

with $C_\oplus := C_⊕\Delta t + C_⊕ + C_\odot$, which is uniformly bounded in $h$ and $\Delta t$.

**Term 2:** $J(u, \Delta t \partial_t^\Theta u)$ Next, we consider the GP term applying Lemma 5.7, Cauchy–Schwarz and a Young’s inequality with $\beta := \sqrt{C_{5,7}}$ (s.t. $C_{5,7}/\beta = \beta$):

$$ \begin{aligned} J(u, \Delta t \partial_t^\Theta u) &\ge -\frac{\beta}{2} \|u\|_J^2 - \frac{\Delta t^2}{2\beta} \|\partial_t^\Theta u\|_J^2 \\ &\ge -\frac{\beta}{2} \|u\|_J^2 - \frac{C_{5,7}}{2\beta} \|u\|_J^2 - \frac{C_{5,7}}{\beta} \frac{\gamma_J(h + \Delta t)}{2} \|\nabla u\|_{Q^h}^2 \ge -C_\oplus \|u\|_j^2, \end{aligned} $$

with $C_\odot := \beta \cdot \max\{1, \gamma_J(h + \Delta t)/2\}$.

**Term 3:** $B_h(u, U)$ Applying Lemma 4.1 and a Cauchy–Schwarz inequality we obtain

$$ B_h(u, U) \le B_{mc}(u, U) + C_{4,1} \cdot (h^{q_s} + \Delta t^{q_t}) \cdot \left( \int_0^T \|u\|_{\partial \Omega^{h(t)}}^2 dt \right)^{\frac{1}{2}} \left( \int_0^T \|U\|_{\partial \Omega^{h(t)}}^2 dt \right)^{\frac{1}{2}}. $$

For $B_{mc}(u, U)$ we note that $U$ is a slice-wise constant function in both space and time, hence the summands involving $\nabla U$ and $\partial_t U$ vanish and we obtain

$$ B_{mc}(u, U) = \sum_{n=1}^{N-1} (u_{-n}^{n}, [\bar{\Pi}(\Delta t \partial_t^\Theta u)]^n)_{\Omega^{h}(t_n)} + (u_{-n}^{N}, \bar{\Pi}(\Delta t \partial_t^\Theta u)_{-n})_{\Omega^{h}(T^{-})} = 0 $$

for $u \in \tilde{W}_h$ through the defining property of $\tilde{W}_h$, cf. (4.11). With a Young’s inequality (with weight $\Delta t^{\frac{1}{2}}$) and Equation (5.5a) we hence have

$$ \frac{-2}{C_{4,1}} B_h(u, U) \ge -(h^{q_s} + \Delta t^{q_t}) \Delta t^{\frac{1}{2}} \cdot C_{5,1} \|u\|_J^2 - (h^{q_s} + \Delta t^{q_t}) \Delta t^{-\frac{1}{2}} \cdot \iint f \|U\|_{\partial \Omega^{h(t)}}^2 dt. $$
We hence turn our attention to $\int_0^T \|U\|_{\partial\Omega^h(t)}^2 dt$, which we first split into the time step contributions and plug in the definition of $\tilde{\Pi}$:

$$ \begin{aligned} \iint_T \|U\|_{\partial\Omega^h(t)}^2 dt &= \sum_{n=1}^N \iint_h \|U\|_{\partial\Omega^h(t)}^2 dt = \sum_{n=1}^N \left( \int_{\Omega^h(t_n)} (\Delta t \partial_t^\Theta u)_-^n dx \right)^2 \iint_n \|1_{I_n}\|_{\partial\Omega^h(t)}^2 dt \\ &\le \sum_{n=1}^N \|(\Delta t \partial_t^\Theta u)_-^n\|_{\partial\Omega^h(t_n)}^2 \int_h^T \frac{\partial \Omega^h(t)}{|\Omega^h(t_n)|} dt \stackrel{\text{Ass. 5.3}}{\le} C_\Omega \sum_{n=1}^N \Delta t^3 \|(\partial_t^\Theta u)_-^n\|_{\partial\Omega^h(t_n)}^2 \le C_\Omega C_{(5.16)} \Delta t \|u\|^2, \end{aligned} $$

where in the last step we applied an estimate as in (5.16) again, but bounded directly (slightly cruder) by $\|u\|^2$. Finally, we have with $C_\circled{3} = \frac{1}{2} C_{4.1} (C_{5.1} + C_\Omega C_{(5.16)})$

$$ -B_h(u, U) \geq -C_\circled{3} (h^{q_s} + \Delta t^{q_t}) \Delta t^{\frac{1}{2}} (\|u\|_c^2 + \Delta t \|\partial_t^\Theta u\|_{Q^h}^2). $$

Putting it all together:

Collecting all the subresults we arrive at

$$ (B_h + J)(u, \tilde{\Pi}(\Delta t \partial_t^\Theta u)) \geq \left(3/4 - C_\circled{3} (h^{q_s} + \Delta t^{q_t}) \Delta t^{\frac{1}{2}}\right) \Delta t \|\partial_t^\Theta u\|_{Q^h}^2 \\ - \left(C_\circled{1} + C_\circled{2} + C_\circled{3} (h^{q_s} + \Delta t^{q_t}) \Delta t^{\frac{1}{2}}\right) \|u\|_j^2. $$

Now, choosing $h$ and $\Delta t$ sufficiently small, cf. Assumption 3.2, such that $C_\circled{3} (h^{q_s} + \Delta t^{q_t}) \Delta t^{\frac{1}{2}} \le 1/4$ and defining $C_{5.9} := C_\circled{1} + C_\circled{2} + C_\circled{3} (h^{q_s} + \Delta t^{q_t}) \Delta t^{\frac{1}{2}}$ we obtain the desired result. $\square$

## 5.5 Stability and continuity

The presented method is stable in the sense of inf-sup-stability. This ensures well-posedness of the discrete problem and is shown in the following:

**THEOREM 5.1.** For all $u \in \tilde{W}_h$ there exists a $\tilde{v}(u) \in \tilde{W}_h$ such that

$$ (B_h + J)(u, \tilde{v}(u)) \geq C_{\text{stab}} \|u\|_j \cdot \|\tilde{v}(u)\|_j. \quad (5.18) $$

*Proof.* Fix $u \in \tilde{W}_h$ and consider $\tilde{v}(u) := (2C_{5.9} + 1)u + \Pi^\perp(\Delta t \partial_t^\Theta u) \in \tilde{W}_h$. Then, we observe

$$ \begin{aligned} (B_h + J)(u, \tilde{v}(u)) &= (2C_{5.9} + 1) \cdot (B_h + J)(u, u) + (B_h + J)(u, \Pi^\perp(\Delta t \partial_t^\Theta u)) \\ &\geq \left(C_{5.9} + \frac{1}{2}\right) \left(\|u\|_j^2 - \frac{C_{4.1}}{2}(h^{q_s} + \Delta t^{q_t}) \cdot \|u\|^2\right) + \frac{\Delta t}{2} \|\partial_t^\Theta u\|_{Q^h}^2 - C_{5.9} \|u\|_j^2 \\ &= \frac{1}{2} \|u\|_j^2 - \left(C_{5.9} + \frac{1}{2}\right) \frac{C_{4.1}}{2} (h^{q_s} + \Delta t^{q_t}) \|u\|^2 \geq \frac{1}{4} \|u\|_j^2, \end{aligned} $$
where in the last step we assumed $\Delta t$ and $h$ sufficiently small, cf. Assumption 3.2, so that the latter summand can be bounded by $\frac{1}{4}\|u\|_j^2$. With (5.15), we observe

$$ \|\tilde{v}(u)\|_j \le (2C_{5,9} + 1)\|u\|_j + \|\Pi^{\perp}(\Delta t \partial_t^{\Theta} u)\|_j \le (2C_{5,9} + 1 + C_{5,8}) \cdot \|u\|_j, \quad (5.19) $$

which implies the claim with $C_{\text{stab}} = \frac{1}{4}(2C_{5,9} + 1 + C_{5,8})^{-1}$. $\square$

We can extend the stability statement to the mass conserving formulation as well:

COROLLARY 5.1. For all $u \in \tilde{W}_h$ there exists a $\tilde{v}(u) \in \tilde{W}_h$ such that

$$ (B_{mc} + J)(u, \tilde{v}(u)) \ge C_{\text{stab}} \|u\|_j \cdot \|\tilde{v}(u)\|_j. \quad (5.20) $$

*Proof.* Reconsidering the previous proof we can rewrite $(B_{mc} + J)(u, \tilde{v}(u))$ as $(B_h + J)(u, \tilde{v}(u)) + (B_{mc} - B_h)(u, \tilde{v}(u))$, where the latter part is bounded by $(h^{qs} + \Delta t^{q_i})\|u\| \cdot \|\tilde{v}(u)\|$ due to Lemma 4.1, Lemma 5.1 and (5.19). Hence, with Assumption 3.2 the claim follows accordingly. $\square$

Moreover, the bilinear forms are also continuous, which we show as follows:

LEMMA 5.10. For all $u \in W_h + H^1(Q^h)$ and all $v \in W_h$

$$ B_h(u, v) \lesssim \|u\|_* \|v\|, \quad (5.21a) $$

$$ B_{mc}(u, v) \lesssim \|u\|_* \|v\|, \quad (5.21b) $$

$$ J(u, v) \lesssim \|u\|_J \|v\|_J. \quad (5.21c) $$

*Proof.* The third claim follows from Cauchy–Schwarz. For the first claim we first apply Lemma 4.1 together with Cauchy–Schwarz and Lemma 5.1

$$ B_h(u, v) \le B_{mc}(u, v) + C_{4,1} (h^{qs} + \Delta t^{q_i}) \|u\|_* \|v\|. $$

For $B_{mc}(u, v)$ we can directly apply Cauchy–Schwarz and then obtain the claim. $\square$

## 6. Strang-type analysis

The results of the previous sections allow us now to prove a Strang-like result. This is of high importance as it gives us a first upper bound on the numerical error.

As the solution of the continuous problem is defined on the exact geometry our mapping $\Psi^{\text{st}}$ translating between $Q^{\text{lin}}$ and $Q$ becomes relevant here. Remember that it was defined such that $\Psi^{\text{st}}(Q^{\text{lin}}) = Q$, paralleling $\Theta_h^{\text{st}}(Q^{\text{lin}}) = Q^h$. Last, also a mapping $\Phi^{\text{st}}: Q^h \to Q, \Phi^{\text{st}} = \Psi^{\text{st}} \circ (\Theta_h^{\text{st}})^{-1}$ was defined.

We define the lifting of the exact solution $u \in H^2(Q)$ to the discrete geometry as

$$ u^l := u \circ \Phi^{\text{st}} = u \circ \Psi^{\text{st}} \circ (\Theta_h^{\text{st}})^{-1}. \quad (6.1) $$
In addition, we lift some functions $q_h$ in the other direction with the definition

$$q_h^{-1} := q_h \circ (\Phi^{\text{st}})^{-1} = q_h \circ \Theta_h^{\text{st}} \circ (\Psi^{\text{st}})^{-1}. \qquad (6.2)$$

As counterparts to the discrete bi- and linear forms we define

$$B(u, v) := (\partial_t u + \mathbf{w} \cdot \nabla u, v)_Q + (\nabla u, \nabla v)_Q + (u_+^0, v_+^0)_{\Omega(0)} \qquad (6.3a)$$

$$f(v) := (f, v)_Q + (u_0, v_+^0)_{\Omega(0)}. \qquad (6.3b)$$

For comparison with the discrete formulation we also split $u = \bar{u} + \tilde{u}$ with $\bar{u} \in \bar{W}_h$ and $\tilde{u}(t_n) = \int_{\Omega^h(t_n)} u \, dx$ and $\int_{\Omega^h(t_n)} \tilde{u} \, dx = \int_{\Omega^h(t_n)} u \, dx = 0$ for $n = 0, \dots, N$, but we do not have $\int_{\Omega^h(t_n)} \tilde{u} \, dx = 0$.

One easily checks that for every test function $v \in W_h + H^2(Q)$ we have the consistency property $B(u, v) = f(v)$ for $u$ the solution to (2.2) and, hence, also $B(\tilde{u}, v) = f(v) - B(\bar{u}, v)$, where $B(\bar{u}, v) = (\partial_t \bar{u}, v)_Q$. This allows us to measure the difference between $\tilde{u}$ and the solution to the discrete problem $\tilde{u}_h = u_h - \bar{u}_h$ as follows:

**THEOREM 6.1** (Strang-type result). Denoting by $u$ the solution to (2.2) with $u = \tilde{u} + \bar{u}$ the unique decomposition introduced above and by $u_h = \tilde{u}_h + \bar{u}_h$ the solution to (4.14) we obtain

$$\|\tilde{u}^l - \tilde{u}_h\| \lesssim \inf_{\tilde{v}_h \in \tilde{W}_h} \|\tilde{u}^l - \tilde{v}_h\| + \|\tilde{u}^l - \tilde{v}_h\|_* + \|\tilde{v}_h\|_J + \sup_{q_h \in W_h} \frac{|f_h(q_h) - f(q_h^{-l})|}{\|q_h\|_J} \\ + \sup_{q_h \in W_h} \frac{|B_h(\bar{u}_h, q_h) - B(\bar{u}, q_h^{-l})|}{\|q_h\|_J} + \sup_{q_h \in W_h} \frac{|B_h(\tilde{u}^l, q_h) - B(\tilde{u}, q_h^{-l})|}{\|q_h\|_J}. \qquad (6.4)$$

**REMARK 6.1.** Note that we have opted for a formulation of Theorem 6.1 involving suprema over $q_h \in W_h$, whereas in a first instance these terms would involve all $q_h \in \tilde{W}_h$. This is because we are not going to exploit the particular structure of $\tilde{W}_h$ later on, and as $\tilde{W}_h \subset W_h$, a straightforward upper bound holds.

*Proof.* Fix an arbitrary $\tilde{v}_h \in \tilde{W}_h$. Then, it holds by the triangle inequality

$$\|\tilde{u}^l - \tilde{u}_h\| \le \|\tilde{u}^l - \tilde{v}_h\| + \|\tilde{v}_h - \tilde{u}_h\|.$$

Concerning the second summand we observe that for some $q_h \in \tilde{W}_h$ chosen as $q_h = \tilde{v}(\tilde{u}_h - \tilde{v}_h)$ as in Theorem 5.1 there holds

$$\|\tilde{u}_h - \tilde{v}_h\| \cdot \|q_h\|_J \le \|\tilde{u}_h - \tilde{v}_h\|_J \cdot \|q_h\|_J \lesssim (B_h + J)(\tilde{u}_h - \tilde{v}_h, q_h) \\ = (B_h + J)(\tilde{u}_h, q_h) - B_h(\tilde{v}_h, q_h) - J(\tilde{v}_h, q_h) = \dots \\ = f_h(q_h) - B_h(\bar{u}_h, q_h)$$
As $u$ was the solution to the continuous problem we can add $0 = B(\tilde{u}, q_h^{-l}) - f(q_h^{-l}) + B(\bar{u}, q_h^{-l})$ and
$0 = -B_h(u^l, q_h) + B_h(u^l, q_h)$ to obtain

$$
\begin{align*}
\dots = f_h(q_h) - B_h(\bar{u}_h, q_h) - B_h(\tilde{v}_h, q_h) - J(\tilde{v}_h, q_h) & \\
& + B(\tilde{u}, q_h^{-l}) - f(q_h^{-l}) + B(\bar{u}, q_h^{-l}) + B_h(\tilde{u}^l, q_h) - B_h(\tilde{u}^l, q_h) & \\
& \le |f_h(q_h) - f(q_h^{-l})| + |B_h(\tilde{u}^l - \tilde{v}_h, q_h)| + |B_h(\bar{u}_h, q_h) - B(\bar{u}, q_h^{-l})| & \\
& \lesssim \|\tilde{u}^l - \tilde{v}_h\|_{\star} \|q_h\|_{\|} \text{ by (5.21a)} & \\
& + |B(\tilde{u}, q_h^{-l}) - B_h(\tilde{u}^l, q_h)| + |J(\tilde{v}_h, q_h)|. &
\end{align*}
$$

Note that the last summand, $J(\tilde{v}_h, q_h)$, can be bounded as $\lesssim \|\tilde{v}_h\|_J \cdot \|q_h\|_j$ by Equation (5.21c). Dividing the overall inequality by $\|q_h\|_j$ and combining it with the first triangle inequality step of this proof we arrive at the result. $\square$

The last three terms in Theorem 6.1 are consistency errors, which are bounded one by one in the next section.

## 6.1 Geometrical consistency errors

In the next Lemmata we will derive bounds for the summands in the previous result, which are related to geometrical consistency error. For notational convenience in the following text we use lifted functions $v^l = v \circ \Phi^{\text{st}}$ for any function $v$ on $Q$ and $v^l = v \circ \Phi^{\text{st}}(\cdot, t)$ for any function $v$ on $\Omega(t)$, $t \in (0, T)$ and with $v^{-l}$ we denote the corresponding inverse lifting. Further, we recall the chain rule:

$$
\begin{align*}
\partial_t v^l &= \partial_t (v \circ \Phi^{\text{st}}) = \partial_t v \circ \Phi^{\text{st}} + (\nabla v \circ \Phi^{\text{st}}) \partial_t \Phi = (\partial_t v)^l + (\nabla v)^l \partial_t \Phi \\
\nabla (v^l) &= D_x (\Phi^{\text{st}})^T (\nabla v)^l, \quad \nabla (w^{-l}) = D_x (\Phi^{\text{st}})^{-T} (\nabla w)^{-l}.
\end{align*}
$$

We now treat the consistency error terms one after another in the following lemmas.

**LEMMA 6.1** (Estimate on right-hand-side consistency error). It holds

$$
\sup_{q_h \in W_h} \frac{|f_h(q_h) - f(q_h^{-l})|}{\|q_h\|_j} \lesssim (h^{q_s} + \Delta t^{q_t+1}) \cdot \mathcal{R}_f \quad \text{with} \quad (6.5a)
$$

$$
\mathcal{R}_f := \|f\|_{W^{1,\infty}(Q)} + \|u_0\|_{W^{1,\infty}(\Omega(t_0))}. \quad (6.5b)
$$

*Proof.* Can be found in Appendix A.3. $\square$

**LEMMA 6.2** (Estimate on consistency of total mass error). With $\mathcal{R}_f$ as in (6.5b) there holds

$$
\|\bar{u}_h - \bar{u}\| = \Delta t^{\frac{1}{2}} \|\partial_t^\Theta (\bar{u}_h - \bar{u})\|_{Q^h} \lesssim \Delta t^{\frac{1}{2}} (h^{q_s} + \Delta t^{q_t}) \cdot \mathcal{R}_f. \quad (6.6a)
$$

$$
\sup_{q_h \in \bar{W}_h} \frac{|B_h(\bar{u}_h, q_h) - B(\bar{u}, q_h^{-l})|}{\|q_h\|_j} \lesssim (h^{q_s} + \Delta t^{q_t}) \cdot \mathcal{R}_f. \quad (6.6b)
$$
*Proof.* As most of the terms in $B_h(\cdot, \cdot)$ and $B(\cdot, \cdot)$ vanish for $\bar{\bar{u}}_h$ and $\bar{u}$ we have

$$B_h(\bar{\bar{u}}_h, q_h) - B(\bar{u}, q_h^{-l}) = (\partial_t \bar{\bar{u}}_h, q_h)_{Q^h} - (\partial_t \bar{u}, q_h^{-l})_{Q^h}.$$

Further, as $\Phi^{\text{st}}$ keeps the time component unaffected we have $\bar{\bar{u}} = \bar{u}^l = \bar{u} \circ \Phi^{\text{st}}$.

$$\begin{aligned}
(\partial_t \bar{\bar{u}}_h, q_h)_{Q^h} - (\partial_t \bar{u}, q_h^{-l})_Q &= (\partial_t \bar{\bar{u}}_h - |\det D\Phi^{\text{st}}| \partial_t \bar{u}, q_h)_{Q^h} \\
&\le |(\partial_t \bar{\bar{u}}_h - \bar{u}), q_h)_{Q^h}| + |((1 - |\det D\Phi^{\text{st}}|) \partial_t \bar{\bar{u}}, q_h)_{Q^h}|
\end{aligned}
\quad (6.7)$$

We start with the second term and introduce the notation $\Omega_n = \Omega(t_n)$ and $\Omega_n^h = \Omega^h(t_n)$. On each time interval $I_n$ we find that $\partial_t \bar{\bar{u}}$ is constant and recalling $(\cdot)^n$ as the spatial integral at time $t^n$ over $\Omega_n$ as in Section 4.2 we find

$$\begin{aligned}
\Delta t \partial_t \bar{\bar{u}}|_{I_n} &= \int_{\Omega_n} u(\cdot, t_n) \, dx - \int_{\Omega_{n-1}} u(\cdot, t_{n-1}) \, dx \\
&= |\Omega_n|^{-1} \bar{u}^n - |\Omega_{n-1}|^{-1} \bar{u}^{n-1} \\
&= \underbrace{|\Omega_{n-1}|^{-1} |\Omega_n|^{-1}}_{\lesssim 1} \left( \underbrace{|\bar{u}^n - \bar{u}^{n-1}| |\Omega_{n-1}| + |\bar{u}^{n-1}| (|\Omega_{n-1}| - |\Omega_n|)}_{\underset{\approx}{=} \mathcal{W}} \right).
\end{aligned}$$

From Reynolds transport theorem we find

$$\mathcal{W} = |\Omega_n| - |\Omega_{n-1}| = \int_{t_{n-1}}^{t_n} \frac{d}{dt} \int_{\Omega_n} 1 \, dx = \int_{t_{n-1}}^{t_n} \int_{\partial \Omega_n} \vec{w} \cdot n_{\partial \Omega} \, ds \, dt,$$

which implies $|\mathcal{W}| \lesssim \Delta t$ and define $\mathcal{W}^h := |\Omega_n^h| - |\Omega_{n-1}^h| = \int_{t_{n-1}}^{t_n} \int_{\partial \Omega_n} \mathcal{W}_h \, ds \, dt$ correspondingly with $\mathcal{W}_h$ as in Section 4.2. Finally, noting that the balance of mass yields $\bar{u}^n - \bar{u}^{n-1} = \int_{\Omega_n} \vec{f}(t) \, dt$, as well as $\bar{u}^n = \bar{u}^0 + \int_0^T \bar{f}(t) \, dt$ we also find $|\bar{u}^n - \bar{u}^{n-1}|$ and $|\bar{u}^n|$ to be bounded by $\mathcal{R}_f$, which yields the bound for the second term in (6.7). We now turn our attention to the first term in (6.7). There holds

$$\begin{aligned}
\Delta t \partial_t (\bar{u} - \bar{u}_h)|_{I_n} &= (\bar{u}^n - \bar{u}^{n-1}) - (\bar{u}_h^n - \bar{u}_h^{n-1}) \\
&= (|\Omega_n|^{-1} \bar{u}^n - |\Omega_{n-1}|^{-1} \bar{u}^{n-1}) - (|\Omega_n^h|^{-1} \bar{u}_h^n - |\Omega_{n-1}^h|^{-1} \bar{u}_h^{n-1}) \\
&= |\Omega_n|^{-1} \underbrace{(\bar{u}^n - \bar{u}^{n-1})}_{= \int_{n-1}^{n} \bar{f}(t) \, dt} + |\Omega_n^h|^{-1} \underbrace{(\bar{u}_h^n - \bar{u}_h^{n-1})}_{= \int_{n-1}^{n} \bar{f}_h(t) \, dt} + \bar{u}^{n-1} \underbrace{\left(\frac{1}{|\Omega_n|} - \frac{1}{|\Omega_{n-1}|}\right)}_{\le h^{qs} + \Delta t^{q,t+1}} + \bar{u}_h^{n-1} \underbrace{\left(\frac{1}{|\Omega_n^h|} - \frac{1}{|\Omega_{n-1}^h|}\right)}_{= \Delta t \le \|f\|_{W^{1,\infty}(Q)}} \\
&\le \underbrace{|\Omega_n^h|^{-1}}_{\le 1} \underbrace{\int_{t_{n-1}}^{t} \max_{t \in I_n} |\bar{f}_h(t) - \bar{f}(t)| \, dt}_{\circled{1}} + \underbrace{|\Omega_n^h|^{-1} - |\Omega_n|^{-1}}_{\le h^{qs} + \Delta t^{q,t+1}} \underbrace{\int_{t_{n-1}}^{t} \max_{t \in I_n} |\bar{f}(t)| \, dt}_{= \Delta t \le \|f\|_{W^{1,\infty}(Q)}} \\
&\quad + \underbrace{|\bar{u}^{n-1} - \bar{u}_h^{n-1}|}_{\circled{2}} \underbrace{|\Omega_n^h|^{-1} - |\Omega_{n-1}^h|^{-1}|}_{\le \Delta t} + \underbrace{|u^{n-1}|}_{\le \mathcal{R}_f} \underbrace{\left(|\Omega_n|^{-1} - |\Omega_{n-1}|^{-1}| - (|\Omega_n^h|^{-1} - |\Omega_{n-1}^h|^{-1})\right)}_{= \circled{3}} \\
&\le \Delta t \cdot (\circled{1} + \circled{2} + (h^{qs} + \Delta t^{q,t+1}) \|f\|_{W^{1,\infty}(Q)}) + \circled{3} \cdot \mathcal{R}_f.
\end{aligned}$$
① we can easily bound using the same techniques as in (A.2) yielding

$$ \textcircled{1} = \max_{t \in I_n} |\bar{f}_h(t) - \bar{f}(t)| \lesssim (h^{qs} + \Delta t^{q_t+1}) \|f\|_{W^{1,\infty}(Q^{hn})}. \quad (6.8) $$

For ② we consider the balance equation and find with the help of (A.1) and (6.8)

$$ \textcircled{2} = |\bar{u}^{n-1} - \bar{u}_h^{n-1}| \le \underbrace{|\bar{u}^0 - \bar{u}_h^0|}_{\lesssim h^{qs} + \Delta t^{q_t+1}} + \iint_{\mathbb{R}}^{t-1} \underbrace{|\bar{f}(t) - \bar{f}_h(t)|}_{\lesssim h^{qs} + \Delta t^{q_t+1}} dt \lesssim (h^{qs} + \Delta t^{q_t+1}) \mathcal{R}_f. $$

For ③ we find

$$ \underbrace{|\Omega_n^h| \cdot |\Omega_{n-1}^h| \cdot |\Omega_n| \cdot |\Omega_{n-1}|}_{\gtrsim 1} \textcircled{3} = |\Omega_n^h| |\Omega_{n-1}^h| \underbrace{(|\Omega_n| - |\Omega_{n-1}|)}_{=\mathcal{W}} - |\Omega_n| |\Omega_{n-1}| \underbrace{(|\Omega_n^h| - |\Omega_{n-1}^h|)}_{=\mathcal{W}^h} \\ \le \underbrace{|\Omega_n^h| |\Omega_{n-1}^h| - |\Omega_n| |\Omega_{n-1}|}_{\textcircled{4}} \underbrace{|\mathcal{W}|}_{\lesssim \Delta t} + \underbrace{|\Omega_n| |\Omega_{n-1}|}_{\lesssim 1} \underbrace{|\mathcal{W} - \mathcal{W}^h|}_{\textcircled{5}} \lesssim \Delta t \cdot \textcircled{4} + \textcircled{6} \text{ with} $$

$$ \textcircled{4} = | |\Omega_n^h| |\Omega_{n-1}^h| - |\Omega_n| |\Omega_{n-1}| | \le \underbrace{|\Omega_n^h| |\Omega_{n-1}^h| - |\Omega_{n-1}| |\Omega_n|}_{\lesssim 1} + \underbrace{| |\Omega_n^h| - |\Omega_n| |\Omega_{n-1}|}_{\lesssim h^{qs} + \Delta t^{q_t+1}} \lesssim 1 $$

$$ \textcircled{5} = |\mathcal{W} - \mathcal{W}^h| \le \iint_{I_h} \iint_{J_{\Omega_h(t)}} |J_{\Phi}(\vec{w} \cdot n_{\partial \Omega})|^l - \gamma_h |ds dt \lesssim \Delta t \|J_{\Phi}(\vec{w} \cdot n_{\partial \Omega})|^l - \gamma_h \|_{\infty, \partial_s Q^{hn}} $$

where $J_\Phi$ is the ratio of the surface measures between $\partial \Omega^h$ and $\partial \Omega$ with $\|1 - J_\Phi\|_\infty \lesssim h^{qs} + \Delta t^{q_t}$, cf. (3.6c) and for any point $s \in \partial_s Q^{hn}$, $t \in I_n$ we have

$$ \begin{aligned} |\gamma_h - (\vec{w} \cdot n_{\partial \Omega})|^l & \le \underbrace{|\gamma_h - (\vec{w} \cdot n_{\partial \Omega^h})| + |\vec{w} \cdot (n_{\partial \Omega^h} - n_{\partial \Omega})|}_{\lesssim h^{qs} + \Delta t^{q_t} \text{ by (4.1)}} + \underbrace{|(\vec{w} - \vec{w}^l)|}_{\lesssim (h^{qs} + \Delta t^{q_t+1}) \|\vec{w}\|_{W^{1,\infty}}} \underbrace{|n_{\partial \Omega}^l|}_{\lesssim 1} \\ & \lesssim h^{qs} + \Delta t^{q_t}. \end{aligned} $$

Hence, $\textcircled{3} \lesssim \Delta t(h^{qs} + \Delta t^{q_t})$, which concludes the proof. □

LEMMA 6.3 (Estimate on left-hand-side consistency error in Strang). With $\mathcal{R}_f$ as in (6.5b) it holds

$$ \sup_{q_h \in W_h} \frac{|B_h(\tilde{u}^l, q_h) - B(\tilde{u}, q_h^{-l})|}{\|q_h\|_j} \lesssim (h^{qs} + \Delta t^{q_t}) (\|u\|_{H^1(Q)} + \mathcal{R}_f). \quad (6.9) $$

Proof. We first split

$$ B_h(\tilde{u}^l, q_h) - B(\tilde{u}, q_h^{-l}) = (B_h(u^l, q_h) - B(u, q_h^{-l})) - (B_h(\bar{u}, q_h) - B(\bar{u}, q_h^{-l})). $$

The latter part is handled as the second r.h.s. term in (6.7) and the former part follows more standard, but more tedious, arguments. For completeness this part is given in Appendix A.4. □
## 6.2 Approximation results

We continue with the overall task of finding upper bounds on all summands of Theorem 6.1. We will analyse an interpolator $I_h$ that yields an upper bound for the approximation term in the Strang lemma

$$ \inf_{\tilde{v}_h \in \bar{W}_h} \|\tilde{u}^l - \tilde{v}_h\| + \|\tilde{u}^l - \tilde{v}_h\|_* + \|\tilde{v}_h\|_J \le \|\tilde{u}^l - I_h \tilde{u}^l\| + \|\tilde{u}^l - I_h \tilde{u}^l\|_* + \|I_h \tilde{u}^l\|_J. \quad (6.10) $$

Assuming $I_h \overline{w}_h = \overline{w}_h$ for all $w_h \in \overline{W}_h$ we observe $\tilde{u}^l - I_h \tilde{u}^l = u^l - \bar{u} - I_h u^l + I_h \bar{u} = u^l - I_h u^l$, i.e., the interpolation error applied on $u^l$ coincides with the one applied on $\tilde{u}^l$. Further, we have $\|I_h \overline{w}_h\|_J = \|\overline{w}_h\|_J = 0$ for all $\overline{w}_h \in \overline{W}_h$ and find $\|I_h \tilde{u}^l\|_J = \|I_h u^l\|_J$ and, hence,

$$ \inf_{\tilde{v}_h \in \bar{W}_h} \|\tilde{u}^l - \tilde{v}_h\| + \|\tilde{u}^l - \tilde{v}_h\|_* + \|\tilde{v}_h\|_J \le \|u^l - I_h u^l\| + \|u^l - I_h u^l\|_* + \|I_h u^l\|_J. \quad (6.11) $$

It hence remains to establish approximation bounds on $W_h$ with an interpolation operator $I_h$ that has $I_h \overline{w}_h = \overline{w}_h$ for all $\overline{w}_h \in \overline{W}_h$.

To exploit the tensor-product structure of the underlying FE space on each time slab we will transfer the problem to the undeformed reference geometry where the geometry description follows a linear-in-space level set function defined on the active tensor-product mesh $Q_{\mathcal{E}}^{\text{lin},n}$.

We will make the following assumption concerning regularity and extendability<fnref n="3" /> of the exact solution.

**ASSUMPTION 6.1** (Regularity and extension of cont. solution $u$). We assume that with $k_s^* := \min\{k_s, q_s\}$, $k_t^* := \min\{k_t, q_t\}$ and $k_{\max} := \max\{k_s^*, k_t^*\}$ the continuous solution $u$ of (2.2) satisfies $u \in H^{k_{\max}+2}(Q)$. Moreover, its extension $u^e$ is bounded in the $H^{k_{\max}+2}$-norm on $\Psi^{\text{st}}(\bigcup_n Q_{\mathcal{E}}^{\text{lin},n})$.

We recapitulate some results from the geometry analysis (Heimann & Lehrenfeld, 2025, section 6) using the notation of t-anisotropic Sobolev spaces introduced in Equation (3.5).

**LEMMA 6.4** There exists an interpolation operator

$$ \Pi_W^n : L^2(Q_{\mathcal{E}}^{\text{lin},n}) \to W_{h, \text{cut}}^{n, k}(Q_{\mathcal{E}}^{\text{lin},n}), \quad (6.12) $$

s.t. for all functions $\hat{u} \in L^2(Q_{\mathcal{E}}^{\text{lin},n})$ with sufficient regularity, s.t. the r.h.s. expressions in the following estimates are well-defined and bounded, there holds for the interpolation error $\hat{e}_u = \hat{u} - \Pi_W^n \hat{u}$ and $\ell_s \in \{1, \dots, k_s+1\}$, $\ell_t \in \{1, \dots, k_t+1\}$:

$$ \begin{align*} \|\hat{e}_u\|_{Q_{\mathcal{E}}^{\text{lin},n}} &\lesssim \Delta t^{\ell_t} \|\hat{u}\|_{H_{\square}^{0, \ell_t}(Q_{\mathcal{E}}^{\text{lin},n})} + h^{\ell_s} \|\hat{u}\|_{H_{\square}^{\ell_s, 0}(Q_{\mathcal{E}}^{\text{lin},n})}, \\ \|\partial_t \hat{e}_u\|_{Q_{\mathcal{E}}^{\text{lin},n}} &\lesssim \Delta t^{\ell_t-1} \|\hat{u}\|_{H_{\square}^{0, \ell_t}(Q_{\mathcal{E}}^{\text{lin},n})} + h^{\ell_s} \|\hat{u}\|_{H_{\square}^{\ell_s, 1}(Q_{\mathcal{E}}^{\text{lin},n})} \quad \text{if } \ell_t \ge 2, \\ \|\nabla \hat{e}_u\|_{Q_{\mathcal{E}}^{\text{lin},n}} &\lesssim \Delta t^{\ell_t} \|\hat{u}\|_{H_{\square}^{1, \ell_t}(Q_{\mathcal{E}}^{\text{lin},n})} + h^{\ell_s-1} \|\hat{u}\|_{H_{\square}^{\ell_s, 0}(Q_{\mathcal{E}}^{\text{lin},n})} \quad \text{if } \ell_s \ge 2, \\ \|\hat{e}_u(\cdot, t_n)\|_{\Omega_{\mathcal{E}}^{\text{lin},n}} &\lesssim \Delta t^{\ell_t-\frac{1}{2}} \|\hat{u}\|_{H_{\square}^{0, \ell_t}(Q_{\mathcal{E}}^{\text{lin},n})} + \Delta t^{-\frac{1}{2}} h^{\ell_s} \|\hat{u}\|_{H_{\square}^{\ell_s, 0}(Q_{\mathcal{E}}^{\text{lin},n})} \quad \text{if } \ell_t \ge 2. \end{align*} $$

<footnote n="3">³ See (Stein, 1971, theorem 5, chapter VI) for a possible particular construction.</footnote>
In the construction of this operator in Preuß (2018) an $L^2$ interpolation in time is concatenated with an $L^2$ interpolation in space, s.t. $\Pi_W^n \overline{w}_h = \overline{w}_h$ holds readily for $\overline{w}_h \in \overline{W}_h$.

In order to obtain a suitable candidate for interpolation on the undeformed domain we define a linear reference geometry counterpart of $u^e$ with $\hat{u} := u^e \circ \Psi^{\text{st}}$, which is defined on (at least) the relevant region of $Q_{\mathbb{G}}^{\text{lin},n}$.

As the previous lemma illustrates it is relevant to obtain boundedness of higher-order derivatives of $\hat{u}$. Fundamentally, this is achieved by making use of the relevant bounds on $\Psi^{\text{st}}$, cf. (3.6f). To show this relation in detail a general Faa di Bruno formula has been applied in order to obtain the following result, cf. (Heimann & Lehrenfeld, 2025, lemma 8.2).

**LEMMA 6.5** Bounds on higher-order spatial or temporal derivatives of $\Psi^{\text{st}}$ imply bounds on higher-order derivatives of $u^e \circ \Psi^{\text{st}}$ by norms of $u^e$ as follows for $\ell_t, \ell_s \in \mathbb{N}$:

$$ \begin{align*} \|\Psi^{\text{st}}\|_{H_{\square}^{0, \ell_t}(Q_{\mathbb{G}}^{\text{lin},n})} &\lesssim 1 \Rightarrow \| u^e \circ \Psi^{\text{st}}\|_{H_{\square}^{0, \ell_t}(Q_{\mathbb{G}}^{\text{lin},n})} \lesssim \|u^e\|_{H^{\ell_t}(\Psi^{\text{st}}(Q_{\mathbb{G}}^{\text{lin},n}))}, \\ \|\Psi^{\text{st}}\|_{H_{\square}^{\ell_s, 0}(Q_{\mathbb{G}}^{\text{lin},n})} &\lesssim 1 \Rightarrow \| u^e \circ \Psi^{\text{st}}\|_{H_{\square}^{\ell_s, 0}(Q_{\mathbb{G}}^{\text{lin},n})} \lesssim \|u^e\|_{H^{\ell_s}(\Psi^{\text{st}}(Q_{\mathbb{G}}^{\text{lin},n}))}, \\ \|\Psi^{\text{st}}\|_{H_{\square}^{\ell_s, 1}(Q_{\mathbb{G}}^{\text{lin},n})} &\lesssim 1 \Rightarrow \| \partial_t(u^e \circ \Psi^{\text{st}})\|_{H_{\square}^{\ell_s, 0}(Q_{\mathbb{G}}^{\text{lin},n})} \lesssim \|u^e\|_{H^{\ell_s+1}(\Psi^{\text{st}}(Q_{\mathbb{G}}^{\text{lin},n}))}, \\ \|\Psi^{\text{st}}\|_{H_{\square}^{1, \ell_t}(Q_{\mathbb{G}}^{\text{lin},n})} &\lesssim 1 \Rightarrow \| \nabla(u^e \circ \Psi^{\text{st}})\|_{H_{\square}^{0, \ell_t}(Q_{\mathbb{G}}^{\text{lin},n})} \lesssim \|u^e\|_{H^{\ell_t+1}(\Psi^{\text{st}}(Q_{\mathbb{G}}^{\text{lin},n}))}. \end{align*} $$

Note that by the boundedness of the (degree-independent) Sobolev extension operator norms such as the Sobolev norm $H^{\ell_t}(\Psi^{\text{st}}(Q_{\mathbb{G}}^{\text{lin},n}))$ can be bounded by their counterparts on the physical domain. Taking together the results from the different time-slices let us denote by $I_h^{\text{lin}}u$ the function defined by $\Pi_W^n \hat{u} = \Pi_W^n(u^e \circ \Psi^{\text{st}})$ on each time slice and $I_h$ the corresponding interpolator on $W_h$:

$$ (I_h^{\text{lin}}\hat{u})|_{Q_{\mathbb{G}}^{\text{lin},n}} = \Pi_W^n(\hat{u}), \quad I_h u^l := I_h^{\text{lin}}\hat{u} \circ (\Theta_h^{\text{st}})^{-1}. \qquad (6.13) $$

We will now derive an upper bound on summands in (6.11):

**LEMMA 6.6** Let $k_t^* := \min\{q_t, k_t\}$, $k_s^* := \min\{q_s, k_s\}$, $k_{\max} := \max\{k_s^*, k_t^*\}$ and $u \in H^{k_{\max}+2}(Q)$. Then,

$$ \|u^l - I_h u^l\| + \|u^l - I_h u^l\|_* \lesssim (\Delta t^{k_t^*+1/2} + \Delta t^{-1/2} h^{k_s^*+1} + h^{k_s^*}) \|u\|_{H^{k_{\max}+2}(Q)}. \qquad (6.14) $$

*Proof.* With $\hat{u} = u^e \circ \Psi^{\text{st}} = u^e \circ \Phi^{\text{st}} \circ \Theta_h^{\text{st}} = (u^e \circ \Phi^{\text{st}}) \circ \Theta_h^{\text{st}} = u^l \circ \Theta_h^{\text{st}}$ and $I_h u^l \circ \Theta_h^{\text{st}} = I_h^{\text{lin}} \hat{u}$ and $e_u^l = u^l - I_h u^l$, $\hat{e}_u = \hat{u} - I_h^{\text{lin}} \hat{u}$ we have $e_u^l \circ \Theta_h^{\text{st}} = \hat{e}_u$ and we aim to bound $\|e_u^l\| + \|e_u^l\|_*$, which we
break down into its contributions:

$$
\begin{align*}
\|e_u^l\|_{Q^h}^2 &\stackrel{(3.7)}{\approx} \sum_{n=1}^N \|\hat{e}_u\|_{Q^{\text{lin},n}}^2 \stackrel{\text{L.6.4\&6.5}}{\lesssim} (\Delta t^{k_t^*+1} \|u^e\|_{H^{k_t^*+1}} + h^{k_s^*+1} \|u^e\|_{H^{k_s^*+1}})^2, \\
\|\partial_t^\Theta e_u^l\|_{Q^h}^2 &\stackrel{(3.16)}{\approx} \sum_{n=1}^N \|\partial_t \hat{e}_u\|_{Q^{\text{lin},n}}^2 \stackrel{\text{L.6.4\&6.5}}{\lesssim} (\Delta t^{k_t^*} \|u^e\|_{H^{k_t^*+1}} + h^{k_s^*+1} \|u^e\|_{H^{k_s^*+2}})^2, \\
\|\nabla e_u^l\|_{Q^h}^2 &\stackrel{(3.7)}{\approx} \sum_{n=1}^N \|\nabla \hat{e}_u\|_{Q^{\text{lin},n}}^2 \stackrel{\text{L.6.4\&6.5}}{\lesssim} (\Delta t^{k_t^*+1} \|u^e\|_{H^{k_t^*+2}} + h^{k_s^*} \|u^e\|_{H^{k_s^*+1}})^2, \\
\|\|e_u^l\|\|^2 &\lesssim \|\|e_u^l\|\|\|_*^2 \lesssim \sum_{i=1}^{N-1} (\| (e_u^l)_+^n \|_{\Omega^h(t_n)}^2 + \| (e_u^l)_-^n \|_{\Omega^h(t_n)}^2) + \sum_{s \in \{0,T\}} \|e_u^l(s)\|_{\Omega^h(s)}^2 \\
&\stackrel{(3.7)}{\approx} \sum_{i=1}^{N-1} (\| (\hat{e}_u)_+^n \|_{\Omega^{\text{lin}}(t_n)}^2 + \| (\hat{e}_u)_-^n \|_{\Omega^{\text{lin}}(t_n)}^2) + \sum_{s \in \{0,T\}} \|\hat{e}_u(s)\|_{\Omega^{\text{lin}}(s)}^2 \\
&\stackrel{\text{L.6.4\&6.5}}{\lesssim} (\Delta t^{k_t^*+1/2} \|u^e\|_{H^{k_t^*+1}} + \Delta t^{-1/2} h^{k_s^*+1} \|u^e\|_{H^{k_s^*+1}})^2.
\end{align*}
$$

Combining these bounds with the continuity of the extension operator we obtain the claim. □

REMARK 6.2. It is also possible to obtain an interpolation bound in terms of $\|u\|_{H^{k_{\text{max}}+1}(Q)}$ at the price of losing half an order in $\Delta t$, i.e.,

$$ \|u^l - I_h u^l\| + \|u^l - I_h u^l\|_* \lesssim (\Delta t^{k_t^*} + \Delta t^{-1/2} h^{k_s^*+1} + h^{k_s^*}) \|u\|_{H^{k_{\text{max}}+1}(Q)} \quad (6.15) $$

This finishes the approximation results involved in Equation (6.11) in relation to the norms $\|\dots\|$ and $\|\dots\|_*$. Hence, we proceed with the last step of bounding $\|I_h u^l\|_J$:

LEMMA 6.7 Let $k_t^* := \min\{q_t, k_t\}$, $k_s^* := \min\{q_s, k_s\}$, $k_{\text{max}} := \max\{k_s^*, k_t^*\}$ and $u \in H^{k_{\text{max}}+1}(Q) \cap W^{2,\infty}(Q)$. Then,

$$
\|I_h u^l\|_J \lesssim (\Delta t^{k_t^*+1} + h^{k_s^*}) \|u\|_{H^{k_{\text{max}}+1}(Q)} \\
+ (\Delta t^{q_t+1} + h^{q_s+1}) \|u\|_{W^{2,\infty}(Q)} + h^{q_s} \|u\|_{W^{1,\infty}(Q)}. \quad (6.16)
$$

*Proof.* First, we write out the l.h.s. expression

$$ \|I_h u^l\|_J^2 = \sum_{n=1}^N \sum_{F \in \mathcal{F}_R^n} \iint_{h_{t-1}}^h \iint_{\omega_F^h(t)} \frac{\tilde{\gamma}_J}{h^2} [I_h u^l]_{\omega_F^h(t)}^2 dx dt. $$

Several technical difficulties in the upcoming proof stem from the fact that we consider isoparametric elements and, hence, $[I_h u^l]_{\omega_F^h(t)}^2$ is in general not a polynomial function and as the mapping changes
between the facet patches several standard arguments such as norm equivalences on finite-dimensional spaces do not apply readily. For the case of straight elements we refer to Preuß (2018) for a simpler proof.

**Strategy.** Before turning to the details of the proof let us briefly outline the strategy. First, we will consider only the spatial integral, fixing one time instance and furthermore restriction to one facet patch. The overall result will then be obtained by summing over all facet patches and integration in time. Secondly, we will bound the $L^2$ norm of the jump that is scaled with $h^{-2}$ by the $L^2$ norm of the *gradient* jump without an additional h-scaling. To achieve that we pay the price of going to a larger neighbourhood of each facet patch. Thirdly, the jump of the interpolants on both sides of the facet will be split into approximation errors between both side's interpolants and the solution for which proper bounds are then derived. As $u^l$ is only piecewise smooth, but not globally on the neighbourhood of the facet patch for which we need to derive bounds, we will furthermore exchange $u^l$ and $u^e$ exploiting geometrical proximity.

**Step 1.** To ease the presentation of the proof we will start with a bit of notation. We fix time step $n$, facet $F \in \mathcal{F}_R^n$, time $t \in I_n$ and in the remainder of the proof will often skip the temporal argument (which is fixed to $t$). We denote the (undeformed) elements aligned to $F$ as $T_1$ and $T_2$ with $\omega := \omega_F^h(t) = \Theta_h^{st}(T_1 \cup T_2, t)$. We extend the mappings from each element to $\mathbb{R}^d$ and define $\Theta_\ell := \mathcal{E}^p \Theta_h^{st}(\cdot, t)|_{T_\ell} \in [\mathcal{P}^{qs}(\mathbb{R}^d)]^d$, $\ell \in \{1, 2\}$. Correspondingly, we define $\Psi_\ell$ as the smooth Sobolev extension of $\Psi^{st}(\cdot, t)|_{T_\ell}$ to $\mathbb{R}^d$. We find $\omega = \Theta_1(T_1) \cup \Theta_2(T_2)$ and note that we will make use of both $\Theta_\ell$ and $\Psi_\ell$, in a neighbourhood of $\omega$.

Recalling the hat-notation for expressions on the undeformed elements, e.g., $\hat{u} = u^e \circ \Psi^{st}$ we extend $I_h^{\text{lin}} \hat{u}|_{T_\ell} = (I_h u^l \circ \Theta_\ell)|_{T_\ell}$ and introduce the short hand notation $\hat{u}_\ell = \mathcal{E}^p I_h^{\text{lin}} \hat{u}|_{T_\ell}$ for the polynomial representation of the interpolant on one undeformed element. With $u_\ell := \hat{u}_\ell \circ \Theta_\ell^{-1}$, $\ell = 1, 2$ we find

$$[I_h u^l]_\omega = u_1 - u_2 = \hat{u}_1 \circ \Theta_1^{-1} - \hat{u}_2 \circ \Theta_2^{-1}$$

and stress that, while $\hat{u}_1$ and $\hat{u}_2$ are polynomials, $u_1$ and $u_2$ and, hence, $[I_h u^l]$ are in general not polynomials, but are nevertheless smooth in a neighbourhood of $\omega$ (for sufficiently small $h$ and $\Delta t$, s.t. $\|id - \Theta_\ell\|$ is small, cf. Assumption 3.2). Further, $[I_h u^l]$ vanishes on $\Theta_1(F) = \Theta_2(F)$.

**Step 2.** Applying the transformation $\Theta_1^{-1}$ on $\omega$ we find $\omega' := \Theta_1^{-1}(\omega) = T_1 \cup \Theta_1^{-1}(\Theta_2(T_2))$, which contains the straight facet $F$. On this domain we define $\xi := [I_h u^l]_\omega \circ \Theta_1 = \hat{u}_1 - \hat{u}_2 \circ \Theta_2^{-1} \circ \Theta_1$. Next, we transform further to a domain of size $\simeq 1$. Let $\Lambda$ be the affine map that maps the reference simplex $\hat{T}$ with diameter $\simeq 1$ onto $T_1$ while simultaneously mapping a unique reference facet $\hat{F}$ (of $\hat{T}$) with diameter $\simeq 1$ onto $F$. Combining $\Lambda$, $\Theta_1$ and shape regularity we find that $\hat{\omega} := \Lambda^{-1}(\omega')$ is a domain that is contained in a ball around the barycentre of $\hat{F}$ of radius $R \simeq 1$, which we denote as $B_{\hat{F}}$ where we stress that $R$ can be chosen uniformly in $F \in \mathcal{F}_R^n$. Then, by shape regularity and Lemma 3.1

$$h^{-2} \| [I_h u^l]_\omega \|_{\omega}^2 \simeq h^{d-2} \| \xi \|_{\omega'}^2 \simeq h^{d-2} \| \xi \circ \Lambda \|_{\hat{\omega}}^2 \le h^{d-2} \| \xi \circ \Lambda \|_{B_{\hat{F}}}^2 \lesssim \dots$$

Now, $\xi = [I_h u^l]_\omega \circ \Theta_1$ is smooth, $\Lambda$ is affine linear, i.e., $\xi \circ \Lambda = [I_h u^l]_\omega \circ \Theta_1 \circ \Lambda$ is smooth and further vanishes on $\hat{F}$. Hence, we can apply a Poincaré-type inequality and scaling arguments and obtain with
$$B'_{F} := \Lambda(B_{\hat{F}}) \text{ and } B_{F} := \Theta_{1}(B'_{F})$$

$$\dots \lesssim h^{d-2} \|\nabla(\xi \circ \Lambda)\|_{B_{\hat{F}}}^{2} \lesssim \|\nabla \xi\|_{B'_{F}}^{2} \approx \|\nabla(\xi \circ \Theta_{1}^{-1})\|_{B_{F}}^{2} = \|\nabla[I_{h}u^{l}]_{\omega}\|_{B_{F}}^{2} \lesssim \|\nabla[I_{h}u^{l}]_{\omega}\|_{B_{\omega}}^{2}.$$

Here, to symmetrize in $\ell \in \{1, 2\}$ in the following text we used $B_{\omega}$ for a small neighbourhood of the patch $\omega$ that contains $B_{F}$ for $\ell \in \{1, 2\}$.

We made the transition from the $L^2$ norm of the difference on $\omega$ to the $H^1$-seminorm of the difference on a neighborhood of $B_{\omega}$ and gained the factor $h^2$ that compensates for the scaling (in $h$) of the GP stabilization.

**Step 3.** Now, we will apply a triangle inequality shifting in three terms: a special (spatial) interpolant $u_1^l$ of $u^e$ based on $\Theta_1(T_1)$, a spatial interpolant $u_2^l$ of $u^e$ based on $\Theta_2(T_2)$ (for both it is essential that they are discrete functions) and $u^e$ itself:

$$\|\nabla[I_h u^l]_\omega\|_{B_\omega}^2 = \|\nabla(u_1 - u_2)\|_{B_\omega}^2 \le \sum_{\ell=1}^2 \left( \underbrace{\|\nabla(u_\ell - u_\ell^l)\|_{B_\omega}^2}_{\textcircled{1}} + \underbrace{\|\nabla(u_\ell^l - u^e)\|_{B_\omega}^2}_{\textcircled{2}} \right).$$

Note that we choose $u^e$ here instead of $u^l$ as $u^l$ is only piecewise smooth, while $u^e$ is assumed to be smooth also across element boundaries. We define $u_\ell^l$ as $u_\ell^l = \hat{u}_\ell^l \circ \Theta_\ell^{-1}$ with $\hat{u}_\ell^l = \mathfrak{F}_h(u^e \circ \Psi_\ell)$ the averaged Taylor polynomial of degree $k_s$ of $u^e \circ \Psi_\ell$ w.r.t. a ball within $L^2(T_\ell)$. This yields a bound for **②**, cf. (Brenner & Scott, 2008, chapter 4):

$$\begin{align*}
\textcircled{2} &= \|\nabla(u_\ell^l - u^e)\|_{B_\omega}^2 \approx \|\nabla(\hat{u}_\ell^l - u^e \circ \Theta_\ell)\|_{\Theta_\ell^{-1}(B_\omega)}^2 \\
&\lesssim \|\nabla(\hat{u}_\ell^l - u^e \circ \Psi_\ell)\|_{\Theta_\ell^{-1}(B_\omega)}^2 + \|\nabla(u^e \circ \Psi_\ell - u^e \circ \Theta_\ell)\|_{\Theta_\ell^{-1}(B_\omega)}^2 \\
&\lesssim h^{2k_s} \|u^e \circ \Psi_\ell\|_{H^{k_s+1}(\Theta_\ell^{-1}(B_\omega))}^2 + \|\Theta_\ell - \Psi_\ell\|_{\infty, \Theta_\ell^{-1}(B_\omega)}^2 \|u^e\|_{W^{2,\infty}(\tilde{B}_\omega)}^2 \\
&\quad + \|D(\Theta_\ell - \Psi_\ell)\|_{\infty, \Theta_\ell^{-1}(B_\omega)}^2 \|\nabla u^e \circ \Psi_\ell\|_{\Theta_\ell^{-1}(B_\omega)}^2 \\
&\stackrel{(3.6a),(3.6f)}{\lesssim} h^{2k_s} \|u^e\|_{H^{k_s+1}(B_\omega)}^2 + (h^{2q_s+2} + \Delta t^{2q_s+2}) \|u^e\|_{W^{2,\infty}(\tilde{B}_\omega)}^2 \\
&\quad + (h^{2q_s} + \Delta t^{2q_s+2}) \|u^e\|_{W^{1,\infty}(\tilde{B}_\omega)}^2,
\end{align*}$$

where we made use of the proximity result for $\|\Theta_\ell - \Psi_\ell\|_{\infty}^2$ and $\|D(\Theta_\ell - \Psi_\ell)\|_{\infty}^2$, which follows from (3.6a) and denoted $\tilde{B}_\omega$ as a small neighbourhood to $B_\omega$. We can now bound **①**, exploiting that $\hat{u}_\ell - \hat{u}_\ell^l = (u_\ell - u_\ell^l) \circ \Theta_\ell$ is a polynomial:

$$\begin{align*}
\textcircled{1} &= \|\nabla(u_\ell - u_\ell^l)\|_{B_\omega}^2 \approx \|\underbrace{\nabla(\hat{u}_\ell - \hat{u}_\ell^l)}_{\in[\mathcal{P}^{k_s-1}]^d} \|_{\Theta_\ell^{-1}(B_\omega)}^2 \approx \|\nabla(\hat{u}_\ell - \hat{u}_\ell^l)\|_{\tilde{T}_\ell}^2 \approx \|\nabla(u_\ell - u_\ell^l)\|_{\Theta_\ell(T_\ell)}^2 \\
&\lesssim \|\nabla(I_h u^l - u^l)\|_{\Theta_\ell(T_\ell)}^2 + \underbrace{\|\nabla(u_\ell^l - u^e)\|_{\Theta_\ell(T_\ell)}^2}_{\textcircled{2}} + \underbrace{\|\nabla(u^e - u^l)\|_{\Theta_\ell(T_\ell)}^2}_{\textcircled{3}}.
\end{align*}$$
With lemma 6.6 yielding a bound for the first term (after integration over $I_n$ and summation over $F \in \mathcal{F}_R^n$ and $n = 1, \dots, N$) we exploit the boundedness of the extension $u^e$, and apply finite overlap arguments to collect all terms together and obtain the claim. □

## 6.3 A priori error bounds

Taking all results from the previous subsections together we obtained suitable bounds for all terms in Theorem 6.1. Hence, we conclude with the following estimate.

**THEOREM 6.2** (A priori error bound in problem specific norm). Let $k_s^* := \min\{q_s, k_s\}$, $k_t^* := \min\{q_t, k_t\}$, $k_{\max} := \max\{k_s^*, k_t^*\}$, and $u \in H^{k_{\max}+2}(Q) \cap W^{2,\infty}(Q)$ be the solution to (2.2) and $u_h = \bar{u}_h + \tilde{u}_h$ the solution to (4.14). There holds

$$ \|u^l - u_h\| \lesssim (h^{q_s} + \Delta t^{q_t}) \mathcal{R}_{f,u} + (h^{k_s^*} + \Delta t^{k_t^*+1/2}) \|u\|_{H^{k_{\max}+2}(Q)} \quad (6.17) $$

with $\mathcal{R}_{f,u} := \mathcal{R}_f + (h + \Delta t)\|u\|_{W^{2,\infty}(Q)}$.

*Proof.* Splitting $\|u^l - u_h\| \le \|\bar{u} - \bar{u}_h\| + \|\tilde{u}^l - \tilde{u}_h\|$ we obtain the sufficient bound (6.6a) for the first term. For the second term, $\|\tilde{u}^l - \tilde{u}_h\|$, we apply Theorem 6.1, together with Lemmas 6.1 to 6.3 for the Strang lemma consistency error summands, and Lemmas 6.6 and 6.7 in relation to the approximation. Further, we note that the regularity in all the terms related to the geometrical consistency can be bounded by the strongest contribution $\|u\|_{W^{2,\infty}(Q)}$ and that the term $\Delta t^{-1/2}h^{k_s+1}$ is bounded by $h^{k_s}$ under Assumption 5.1. □

This result provides reasonable bounds for the discretization (4.14), which are optimal in terms of convergence rates in the specific norm involving $k_t$ and $k_s$, assuming sufficient geometric consistency. However, the result is slightly suboptimal with respect to regularity, at least in part due to the previously chosen space-time interpolation. Furthermore, although higher-order $L^2$-bounds can be deduced directly from the result, these are obviously not optimal.

(4.14) has the computational disadvantage that one needs to solve on $\tilde{W}_h$ for which one would typically introduce a scalar Lagrange multiplier. Although possible, in practice, one would prefer to rather compute on $W_h$ as in (4.5) directly. For this case we have the following remark.

**REMARK 6.3.** If in the splitting $u_h = \bar{u}_h + \tilde{u}_h$ the total mass $\bar{u}_h$ is not chosen as in (4.13), but with a deviation of $\lesssim \Delta t(h^{q_s} + \Delta t^{q_t})$ in the integral approximation of $f_h$, as in (4.9), the results of Lemma 6.3 and then also Theorem 6.2 remain unaffected. This implies that the standard method (4.5) yields the same error estimates once a discrete solution exists.

A sweet spot between a rigorously analysable and computable method is the mass-conserving variant (4.6) as we will treat next.

**COROLLARY 6.1.** Let $k_s^* := \min\{q_s, k_s\}$, $k_t^* := \min\{q_t, k_t\}$, $k_{\max} := \max\{k_s^*, k_t^*\}$, and $u \in H^{k_{\max}+2}(Q) \cap W^{2,\infty}(Q)$ be the solution to (2.2) and $u_{mc}$ the solution to (4.6). There holds the error bound

$$ \|u^l - u_{mc}\| \lesssim (h^{q_s} + \Delta t^{q_t}) (\mathcal{R}_f + \|u\|_{W^{2,\infty}(Q)}) + (h^{k_s^*} + \Delta t^{k_t^*+1/2}) \|u\|_{H^{k_{\max}+2}(Q)}. \quad (6.18) $$
*Proof.* In the Strang lemma using Lemmas 4.1 and 5.1 we can replace $B_h$ by $B_{mc}$ again at the price of a term of the form $(h^{q_s} + \Delta t^{q_t})(\|\tilde{u}\| + \|\overline{\tilde{u}}_h\|)$ where both parts are bounded by (up to a constant) $(h^{q_s} + \Delta t^{q_t}) \cdot \|u\|_{W^{2,\infty}(Q)}$. $\square$

With Corollary 6.1 we find that the mass conserving formulation yields a method that is completely analysable with reasonable error bounds *without the need* to implement the splitting $W_h = \overline{W}_h + \tilde{W}_h$ explicitly.

**COROLLARY 6.2.** The bounds in Theorem 6.2 and Corollary 6.1 also hold for $\|(u^l - u_h)(T)\|_{\Omega^h(T)}$ and $\|(u^l - u_{mc})(T)\|_{\Omega^h(T)}$, as, by definition, the summand is contained in $\|\cdot\|$.

Before turning to a comparison with related results in the literature in the next remark let us first state a very basic observation: in the simplest case $k_t = q_t = k_s = q_s$ the estimates is optimal for the case of simultaneous space-time refinements. Note that in this case the interpolation bound from Remark 6.2 may be employed, which lowers the regularity assumption on the exact solution from $k_{\max}+2$ to $k_{\max}+1$.

**REMARK 6.4** (Comparison to the literature). Let us finally compare these theoretical *a priori* estimates to the results of the counterpart analysis with assumption of an exact handling of geometries, Preuß (2018), as well as to related results from the literature, Lehrenfeld & Reusken (2013) and Badia *et al.* (2023), as already mentioned in Remark 5.1.

Starting with Preuß (2018) we observe that the fundamental inf-sup-stability proof approach conceptually applies both with exact or discrete geometries assumed. However, the fact that $B_h(\cdot, \cdot) \neq B_{mc}(\cdot, \cdot)$ leads to a perturbation, cf. (4.8a), that is hard to control with the $\|\cdot\|$ norm as it does not control slabwise constants. This led to the introduction of the splitting $u_h = \overline{u}_h + \tilde{u}_h$ in Section 4.2 and the use of Lemma 5.8. A different approach to deal with the constants, discussed in Heimann (2025) and Wendler (2022), are weak penalty or Lagrange multiplier approaches for the control of slabwise constants. For the approximation results the fact that we have to deal with *mapped* polynomials in the GP stabilization increased the complexity of the proof. Nevertheless, driving the geometry errors sufficiently small we retrieve the same error bound as in Preuß (2018) $\|u - u_h\| \lesssim (h^{k_s} + \Delta t^{k_t+\frac{1}{2}}) \|u\|_{H^{k_{\max}+2}(Q)}$.

In the works Lehrenfeld & Reusken (2013) and Badia *et al.* (2023) similar problems have been considered, which however circumvent some of the above-mentioned problems by assuming exact geometry handling and different boundary conditions, which would allow to control slabwise constants directly. Moreover, Badia *et al.* (2023) and Lehrenfeld & Reusken (2013) do not apply GPs and consider weaker norms, similar to $\|\cdot\|_c$, which are tailored for a much simpler coercivity based stability analysis. The bounds obtained are also optimal for $\Delta t \sim h$, but less flexible in space-time anisotropic cases. The regularity assumptions in Badia *et al.* (2023) on the solution $u$ are however explicated in more detail insofar as spatial and temporal derivatives are distinguished. Finally, in Reusken & Sass (2024), for the case of a parabolic *surface* PDE on a moving surface, a coercivity-based analysis with $\Delta t \sim h$ was considered to derive *a priori* error bounds.

**REMARK 6.5** (Numerical examples). As it relates to the numerical investigation of the presented *a priori* error bounds we want to refer the reader to the literature: first, in our paper Heimann *et al.* (2023), we presented numerical examples for different variants of the method. Concerning the imposition of the mass conservation Heimann (2025) contains numerical studies of the realization by a penalty approach (cf. (4.16)). An implementation of a Lagrange multiplier can be found furthermore in Heimann *et al.* (2025). We emphasize that these variant methods only differed in minor details in those works.
Acknowledgements

We thank two reviewers for their helpful suggestions to improve the presentation.

Funding

Engineering and Physical Sciences Research Council (EP/X042650/1 to F.H.).

REFERENCES

BADIA, S., DILIP, H. & VERDUGO, F. (2023) Space-time unfitted finite element methods for time-dependent problems on moving domains. *Comput. Math. Appl.*, **135**, 60–76.

BERNARDI, C. (1989) Optimal finite-element interpolation on curved domains. *SIAM J. Numer. Anal.*, **26**, 1212–1240.

BRENNER, S. C. & SCOTT, L. R. (2008) *The Mathematical Theory of Finite Element Methods*. New York, NY: Springer.

BURMAN, E. (2010) Ghost penalty. *C.R. Math. Acad. Sci. Paris*, **348**, 1217–1220.

BURMAN, E., CLAUS, S., HANSBO, P., LARSON, M. G. & MASSING, A. (2015) CutFEM: discretizing geometry and partial differential equations. *Int. J. Numer. Methods Eng.*, **104**, 472–501.

CIARLET, P. & RAVIART, P.-A. (1972) Interpolation theory over curved elements, with applications to finite element methods. *Comput. Methods Appl. Mech. Eng.*, **1**, 217–249.

HEIMANN, F. (2020) *On Discontinuous- and Continuous-In-Time Unfitted Space-Time Methods for PDEs on Moving Domains* Master’s thesis. GRO.data, Göttingen: University of Göttingen. https://doi.org/10.25625/CDCMYT.

HEIMANN, F. (2025) *Higher Order Unfitted Space-Time Finite Element Methods for Moving Domain Problems* PhD thesis. SUB Göttingen, Göttingen: University of Göttingen. https://doi.org/10.53846/goediss-11003.

HEIMANN, F. & LEHRENFELD, C. (2019) Numerical integration on hyperrectangles in isoparametric unfitted finite elements. *European Conference on Numerical Mathematics and Advanced Applications* (F. A. RADU, K. KUMAR, I. BERRE, J. M. NORDBOTTEN, I. S. POP eds). Cham: Springer, pp. 193–202.

HEIMANN, F. & LEHRENFELD, C. (2025) Geometry error analysis of a parametric mapping for higher order unfitted space-time methods. *IMA J. Numer. Anal.*, drae098.

HEIMANN, F., LEHRENFELD, C. & PREUSS, J. (2023) Geometrically higher order unfitted space-time methods for PDEs on moving domains. *SIAM J. Sci. Comput.*, **45**, B139–B165.

HEIMANN, F., LEHRENFELD, C. & PREUSS, J. (2025) *Discretization Error Analysis of a High Order Unfitted Space-Time Method for Moving Domain Problems*. Archive of software accompanying this paper. GRO.data, Göttingen. https://doi.org/10.25625/K9HFBR.

LARSON, M. G., LOGG, A. & LUNDHOLM, C. (2024) Space-time CutFEM on overlapping meshes I: simple continuous mesh motion. *Numer. Math.*, **156**, 1015–1054.

LEHRENFELD, C. (2016) High order unfitted finite element methods on level set domains using isoparametric mappings. *Comput. Methods Appl. Mech. Eng.*, **300**, 716–733.

LEHRENFELD, C. (2017) A higher order isoparametric fictitious domain method for level set domains. *Geometrically Unfitted Finite Element Methods and Applications* (S. P. A. BORDAS, E. BURMAN, M. G. LARSON & M. A. OLSHANSKII eds). Cham: Springer, pp. 65–92.

LEHRENFELD, C. & REUSKEN, A. (2013) Analysis of a Nitsche XFEM-DG discretization for a class of two-phase mass transport problems. *SIAM J. Numer. Anal.*, **51**, 958–983.

LEHRENFELD, C. & REUSKEN, A. (2018) Analysis of a high order unfitted finite element method for an elliptic interface problem. *IMA J. Numer. Anal.*, **38**, 1351–1387.

LOU, Y. & LEHRENFELD, C. (2022) Isoparametric unfitted BDF—finite element method for PDEs on evolving domains. *SIAM J. Numer. Anal.*, **60**, 2069–2098.

MYRBÄCK, S. & ZAHEDI, S. (2024) A high-order conservative cut finite element method for problems in time-dependent domains. *Comput. Methods Appl. Mech. Eng.*, **431**, 117245.

PREUSS, J. (2018) *Higher Order Unfitted Isoparametric Space-Time FEM on Moving Domains* Master’s thesis. GRO.data, Göttingen: University of Göttingen. https://doi.org/10.25625/UACWXS.
REUSKEN, A. & SASS, H. (2024) *Analysis of a Space-Time Unfitted Finite Element Method for PDEs on Evolving Surfaces*. Preprint: https://arxiv.org/abs/2401.01215.

SCHIEWECK, F. (2010) A-stable discontinuous Galerkin–Petrov time discretization of higher order. *J. Numer. Math.*, **18**, 25–57.

SETHIAN, J. A. (1999) *Level Set Methods and Fast Marching Methods*, vol. **98**. Cambridge: Cambridge UP.

STEIN, E. M. (1971) *Singular Integrals and Differentiability Properties of Functions (PMS-30)*, vol. **30**. Princeton: Princeton University Press.

WARBURTON, T. & HESTHAVEN, J. (2003) On the constants in hp-finite element trace inverse inequalities. *Comput. Methods Appl. Mech. Eng.*, **192**, 2765–2773.

WENDLER, A. C. (2022) *Monolithic Unfitted Space-Time Fem for an Osmotic Cell Swelling Problem* Master’s thesis. GRO.data, Göttingen: University of Göttingen. https://doi.org/10.25625/0KPEON.

## Appendix A. Detailed proofs of some lemmata

### A.1 *Proof of lemma 5.3*

To $t \in I_n$ we define $T_i^h := \Theta_h(T_i, t)$, $i = 1, 2$. Let us start with $l_s = 0$. With $v = \hat{v} \circ (\Theta_h(\cdot, t))^{-1}$ and (3.7) we find $\|v\|_{T_i^h} \simeq \|\hat{v}\|_{T_i}$. In particular, $\|v\|_{T_i^h} \lesssim \|\hat{v}\|_{T_i}$. That motivates to start our proof by deriving a representation for $\hat{v}|_{T_1}$ with the help of the facet patch jump: On $T_1$,

$$[v]_{\omega_F^h} \circ (\Theta_h(\cdot, t)) = \hat{v} - \mathcal{E}^p(\hat{v}|_{T_2}) \circ (\mathcal{E}^p \Theta_h(\cdot, t)|_{T_2})^{-1} \circ (\Theta_h(\cdot, t)).$$

In the following text we abbreviate $\Theta_h(\cdot, t)$ as $\Theta_h(t)$. We will also write $w_{T_2}^{\mathcal{E}}$ to denote the extension $\mathcal{E}^p w|_{T_2}$ of a polynomial $w$ on $T_2$ to $\mathbb{R}^d$. In particular, we denote $\mathcal{E}^p \Theta_h(t)|_{T_2}$ as $\Theta_{h,T_2}^{\mathcal{E}}(t)$. Then, rearranging and integrating over $T_1$ combined with Young’s inequality yields

$$\|v\|_{T_1^h}^2 \lesssim \|\hat{v}\|_{T_1}^2 \lesssim \underbrace{\|[v]_{\omega_F^h} \circ (\Theta_h(t))\|_{T_1}^2}_{\simeq \|[v]_{\omega_F^h} \|_{T_1^h}^2 \text{ by (3.7)}} + \underbrace{\|\hat{v}_{T_2}^{\mathcal{E}} \circ (\Theta_{h,T_2}^{\mathcal{E}}(t))^{-1} \circ (\Theta_h(t))\|_{T_1}^2}_{=:I}.$$

For $I$ we exploit shape regularity and standard scaling arguments to get

$$I \lesssim h^{d/2} \|\hat{v}_{T_2}^{\mathcal{E}} \circ (\Theta_{h,T_2}^{\mathcal{E}}(t))^{-1} \circ (\Theta_h(t))\|_{\infty, T_1} = h^{d/2} \|\hat{v}_{T_2}^{\mathcal{E}}\|_{\infty, T_1^*} =: II,$$

where $T_1^* = (\Theta_{h,T_2}^{\mathcal{E}}(t))^{-1} \circ (\Theta_h(t))(T_1)$. Next, we define two balls centred around the barycentre $x_2^*$ of $T_2$, $B_A = B(x_2^*, R_A)$ and $B_B = B(x_2^*, R_B)$ with radii $R_A, R_B$ s.t. $B_A \subset T_2 \subset (T_1^* \cup T_2) \subseteq B_B$. Then, we can conclude

$$h^{-d/2} \cdot (II) \lesssim \|\hat{v}_{T_2}^{\mathcal{E}}\|_{\infty, B_B} \stackrel{(*)}{\lesssim} \|\hat{v}_{T_2}^{\mathcal{E}}\|_{\infty, B_A} \lesssim \|\hat{v}_{T_2}^{\mathcal{E}}\|_{\infty, T_2} = \|\hat{v}\|_{\infty, T_2} \stackrel{(**)}{\lesssim} h^{-d/2} \|v\|_{T_2^h},$$

where we exploited norm equivalences on finite dimensional spaces in (*) and scaling arguments and (3.7) in (**). This concludes the proof for the $L^2$ norm.
For the case $l_s = 1$ we introduce $v^0 = \frac{1}{T_2} \int_{T_2} \hat{v} \, ds \in \mathbb{R}$ (with $[v^0]_{\omega_F^h} = 0$ and $\nabla v^0 = 0$) and use standard inverse inequalities and interpolation results to obtain

$$ 
\begin{aligned}
h^2 \|\nabla v\|_{\Theta_h(T_1,t)}^2 &= h^2 \|\nabla(v-v^0)\|_{T_1^h}^2 \lesssim \|(v-v^0)\|_{T_1^h}^2 \lesssim \|[v-v^0]_{\omega_F^h}\|_{T_1^h}^2 + \|v-v^0\|_{T_2^h}^2 \\
&\lesssim \|[v-v^0]_{\omega_F^h}\|_{\Theta_h(T_1,t)}^2 + \|\hat{v}-v^0\|_{T_2}^2 \lesssim \|[v]_{\omega_F^h}\|_{\Theta_h(T_1,t)}^2 + h^2 \|\nabla \hat{v}\|_{T_2}^2 \\
&\lesssim \|[v]_{\omega_F^h}\|_{\Theta_h(T_1,t)}^2 + h^2 \|\nabla v\|_{T_2^h}^2,
\end{aligned}
\quad \square $$

where we used the scaled Poincaré estimate $\|\hat{v} - v^0\|_{T_2} \lesssim h \|\nabla \hat{v}\|_{T_2}^2$.

## A.2 Proof of lemma 5.4

We start with (5.10) and $l_s = 0$: as $Q_{\mathcal{E}}^{h,n}$ was defined as the image of $Q_{\mathcal{E}}^{\text{lin},n}$ under the mapping $\Theta_h^{\text{st}}$ we can write

$$ \|u\|_{Q_{\mathcal{E}}^{h,n}}^2 = \|u\|_{Q_{\mathcal{F}}^{h,n}}^2 + \sum_{T \in \Omega_{\mathcal{E}}^{\text{lin},I_n} \setminus \Omega_{\mathcal{F}}^{\text{lin},I_n}} \iint_{I_n} \|u\|_{\Theta_h(T, \tau)}^2 d\tau. $$

Comparing this with the final result it remains to estimate the second term. To this end fix an element $T \subseteq \Omega_{\mathcal{E}}^{\text{lin},I_n} \setminus \Omega_{\mathcal{F}}^{\text{lin},I_n}$. By Assumption 5.4 we find an element $\mathcal{B}(T) \subseteq \Omega_{\mathcal{F}}^{\text{lin},I_n}$ and a path between $T$ and $\mathcal{B}(T)$ of finite length, $\{T_i\}_{i=0}^M$, $T_0 = T$, $T_M = \mathcal{B}(T)$. If we traverse along this path we can apply Lemma 5.3 in each step in order to obtain for each $T \subseteq \Omega_{\mathcal{E}}^{\text{lin},I_n} \setminus \Omega_{\mathcal{F}}^{\text{lin},I_n}$ and time $\tau \in I_n$

$$ \|u\|_{\Theta_h(T, \tau)}^2 \lesssim \|u\|_{\Theta_h(\mathcal{B}(T), \tau)}^2 + \sum_{F \in \mathcal{F}_T^n} \|[u]_{\omega_F^h}\|_{\omega_F^h(\tau)}^2, $$

where $\mathcal{F}_T^n$ is the set of facets that are crossed on the path from $T$ to $\mathcal{B}(T)$. If we integrate this equation in time $\tau$ over $I_n$ and sum up over all elements from $\Omega_{\mathcal{E}}^{\text{lin},I_n} \setminus \Omega_{\mathcal{F}}^{\text{lin},I_n}$ we end up with

$$ 
\begin{aligned}
\|u\|_{Q_{\mathcal{E}}^{h,n}}^2 &\lesssim \sum_{T \in \Omega_{\mathcal{F}}^{\text{lin},I_n}} (1 + \#\mathcal{B}^{-1}(T)) \iint_{I_n} \|u\|_{\Theta_h(T, \tau)}^2 d\tau + \sum_{F \in \mathcal{F}_R^n} N_F \iint_{I_n} \|[u]_{\omega_F^h}\|_{\omega_F^h(\tau)}^2 d\tau \\
&\lesssim \|u\|_{Q_{\mathcal{F}}^{h,n}}^2 + \left(1 + \frac{\Delta t}{h}\right) \sum_{F \in \mathcal{F}_R^n} \iint_{I_n} \|[u]_{\omega_F^h}\|_{\omega_F^h(\tau)}^2 d\tau \lesssim \|u\|_{Q_{\mathcal{F}}^{h,n}}^2 + \frac{h^2}{\gamma_J} j_h^n(u, u).
\end{aligned}
$$

This finishes the proof for the first estimate of $\|u\|_{Q_{\mathcal{E}}^{h,n}}^2$. The case $l_s = 1$ follows similar lines.
### A.3 Proof of Lemma 6.1

Recalling $q_h = q_h^{-l} \circ \Phi^{\text{st}}$ and $f^l = f \circ \Phi^{\text{st}}$, $u_0^l = u_0 \circ \Phi^{\text{st}}(\cdot, t_0)$ we have

$$ 
\begin{aligned}
f(q_h^{-l}) &= (f, q_h^{-l})_Q + (u_0, (q_h^{-l})_+^0)_{\Omega(t_0)} \\
&= (|\det D\Phi^{\text{st}}[f, q_h]_{Q_h}| + (|\det D_x\Phi(\cdot, t_0)|u_0^l, q_{h,+}^0)_{\Omega^{h}(t_0)}) \\
f_h(q_h) &= (f^e, q_h)_{Q^h} + (u_0, q_{h,+}^0)_{\Omega^{h}(t_0)}.
\end{aligned}
$$

We start by estimating the initial data difference

$$ 
\begin{aligned}
|(\det D_x \Phi(\cdot, t_0) |u_0^l - u_0, q_h^0)_{\Omega^{h}(t_0)}| &= \underbrace{(1 - \det D_x \Phi(\cdot, t_0))}_{|\cdot| \lesssim h^{qs} + \Delta t^{q_t+1} \text{ with (3.6b)}} |u_0^l + (\underbrace{(u_0 - u_0^l)}_{|\cdot| \lesssim \|u_0\|_{W^{1,\infty}} \|\text{Id} - \Phi(\cdot, t_0)\|_{\infty}} , q_h^0)_{\Omega^{h}(t_0)}| \\
&\lesssim (h^{qs} + \Delta t^{q_t+1}) \|u_0\|_{W^{1,\infty}(\Omega(t_0))} \|q_h\|_j.
\end{aligned}
\quad (A.1) $$

where we have used (3.6a) in the last step and note that $\|q_h\|_j$ contains the initial data as a summand in $\|q_h\|$. For the remaining terms we obtain

$$ 
\begin{aligned}
|(f^e - |\det D\Phi^{\text{st}}[f^l, q_h]_{Q^h}|) \lesssim & (\underbrace{(1 - |\det D\Phi^{\text{st}}|)}_{|\cdot| \lesssim h^{qs} + \Delta t^{q_t+1} \text{ with (3.6b)}} f^l, q_h)_{Q^h}| + |(|f^e - f^l|, q_h)_{Q^h}| \\
& \le \|f^e\|_{W^{1,\infty}(Q^h)} \cdot \|\text{Id} - \Phi^{\text{st}}\|_{\infty, Q^h} \lesssim (h^{qs+1} + \Delta t^{q_t+1}) \|f^e\|_{W^{1,\infty}(Q^h)}.
\end{aligned}
\quad (A.2) $$

with

$$ \|f^e - f^l\|_{\infty, Q^h} \lesssim \|f^e\|_{W^{1,\infty}(Q^h)} \cdot \|\text{Id} - \Phi^{\text{st}}\|_{\infty, Q^h} \lesssim (h^{qs+1} + \Delta t^{q_t+1}) \|f^e\|_{W^{1,\infty}(Q^h)}. $$

In total, we arrive at

$$ 
\begin{aligned}
|f_h(q_h) - f(q_h^{-l})| &\lesssim (h^{qs} + \Delta t^{q_t+1}) \|f\|_Q \cdot \|q_h\|_{Q^h} \\
&\quad + (h^{qs+1} + \Delta t^{q_t+1}) \|f^e\|_{W^{1,\infty}(Q)} \|q_h\|_{Q^h} + (h^{qs} + \Delta t^{q_t+1}) \|u_0\|_{W^{1,\infty}(\Omega(t_0))} \|q_h\|_j.
\end{aligned}
$$

Using $\|q_h\|_{Q^h} \lesssim \|q_h\|_j$ from (5.5b) we arrive at the result.

### A.4 Proof of lemma 6.3

Writing out $B_h(u^l, q_h)$ and $B(u, q_h^{-l})$ and exploiting $[u^l]^n = 0$ we have

$$ 
\begin{aligned}
B_h(u^l, q_h) &= \underbrace{(\partial_t(u^l), q_h)_{Q^h}}_{\textcircled{1}} + \underbrace{(\vec{w}\nabla(u^l), q_h)_{Q^h}}_{\textcircled{2}} + \underbrace{(\nabla(u^l), \nabla q_h)_{Q^h}}_{\textcircled{3}} + \underbrace{(u_0^l, (q_h)_+^0)_{\Omega^{h}(0)}}_{\textcircled{4}} \\
B(u, q_h^{-l}) &= \underbrace{(|\det D\Phi^{\text{st}}|(\partial_t u)^l, q_h)_{Q^h}}_{\textcircled{1}} + \underbrace{(|\det D\Phi^{\text{st}}|(\vec{w}\nabla u)^l, q_h)_{Q^h}}_{\textcircled{2}} \\
&\quad + \underbrace{(|\det D\Phi^{\text{st}}|(\nabla u)^l, (\nabla(q_h^{-l})^l)_{Q^h})}_{\textcircled{3}} + \underbrace{(|\det D_x\Phi^{\text{st}}(\cdot, 0)|u_0^l, (q_h)_+^0)_{\Omega^{h}(0)}}_{\textcircled{4}}.
\end{aligned}
$$
We can estimate the differences as follows:

$$
\begin{aligned}
|\mathring{\text{a}} - \mathring{\text{b}}| & \le \left( (\nabla u)^l \underbrace{\partial_t \Phi}_{\|\cdot\|_{\infty} \lesssim h^{qs+1} + \Delta t^{qt}} - \underbrace{|1 - \det D\Phi^{\text{st}}|}_{\lesssim h^{qs} + \Delta t^{qt+1}} (\partial_t u)^l, q_h \right)_{Q_h} \\
& \lesssim \|q_h\|_{Q^h} (h^{qs} + \Delta t^{qt+1}) \|\partial_t u\|_{Q} + (h^{qs+1} + \Delta t^{qt}) \|\nabla u\|_{Q}.
\end{aligned}
$$

Next, we continue with the second term

$$
\mathring{\text{a}} - \mathring{\text{b}} = \left( \left( \vec{w} D_x (\Phi^{\text{st}})^T - \vec{w}^l |\det D\Phi^{\text{st}}| \right) (\nabla u)^l, q_h \right)_{Q^h}.
$$

We observe that for the inner term

$$
\begin{aligned}
|\vec{w} D_x (\Phi^{\text{st}})^T - \vec{w}^l |\det D\Phi^{\text{st}}|| & \le |\vec{w} D_x (\Phi^{\text{st}})^T - \vec{w}^l| + |\vec{w}^l||1 - \det D\Phi^{\text{st}}| \\
& \lesssim (h^{qs+1} + \Delta t^{qt+1}) \|w\|_{W^{1,\infty}} \text{ by } (3.6a) \\
& \qquad + \underbrace{|(I - D_x (\Phi^{\text{st}})^T)|}_{\lesssim h^{qs} + \Delta t^{qt+1} \text{ by } (3.6a)} |\vec{w}| + \underbrace{|\vec{w}^l|}_{\lesssim h^{qs} + \Delta t^{qt+1} \text{ by } (3.6b)} |1 - \det D\Phi^{\text{st}}|.
\end{aligned}
$$

In total, with $\|\vec{w}\|_{W^{1,\infty}} \lesssim 1$ we arrive at

$$
| \mathring{\text{a}} - \mathring{\text{b}} | \lesssim (h^{qs} + \Delta t^{qt+1}) \| \nabla u \|_{Q} \| q_h \|_{Q^h}.
$$

Next, we estimate the third difference

$$
\mathring{\text{a}} - \mathring{\text{b}} = \left( \left( (D_x \Phi^{\text{st}})^T - |\det D\Phi^{\text{st}}|(D_x \Phi^{\text{st}})^{-1} \right) (\nabla u)^l, \nabla q_h \right)_{Q^h}.
$$

Bounding the inner term we find

$$
\begin{aligned}
& \| D_x (\Phi^{\text{st}})^T - |\det D\Phi^{\text{st}}|(D_x \Phi^{\text{st}})^{-1} \| \\
& \le \underbrace{\| I - D_x (\Phi^{\text{st}})^T \|}_{\lesssim h^{qs} + \Delta t^{qt+1} \text{ by } (3.6a)} + \underbrace{|1 - |\det D\Phi^{\text{st}}||}_{\lesssim h^{qs} + \Delta t^{qt+1} \text{ by } (3.6b)} \underbrace{\| (D_x \Phi^{\text{st}})^{-1} \|}_{\lesssim 1} + \underbrace{\| (D_x \Phi^{\text{st}})^{-1} - I \|}_{\lesssim h^{qs} + \Delta t^{qt+1} \text{ by } (3.6a)} \\
& \Rightarrow |\mathring{\text{a}} - \mathring{\text{b}}| \lesssim (h^{qs} + \Delta t^{qt+1}) \cdot \|\nabla u\|_{Q} \cdot \|q_h\|_{j}.
\end{aligned}
$$

For the last pair of summands we note

$$
|\mathring{\text{a}} - \mathring{\text{b}}| \lesssim \underbrace{(|1 - \det D_x \Phi^{\text{st}}(\cdot, 0)| u_0^l, (q_h)_+^0)_{\Omega^{h(0)}}}_{\lesssim h^{qs} + \Delta t^{qt+1} \text{ by } (3.6b)} \lesssim (h^{qs} + \Delta t^{qt+1}) \underbrace{\|u_0^e\|_{\Omega^{h(0)}} \|q_h\|_j}_{\lesssim \|u_0\|_{\Omega^{(0)}}}
$$

Using $\|u_0\|_{\Omega(0)} = \|u\|_{\Omega(0)} \lesssim \|u\|_{H^1(Q)}$ and collecting the subresults yields the claim.