The aim of this paper is to develop a numerical scheme to approximate evolving interface problems for parabolic equations based on the abstract evolving finite element framework proposed in [20]. An appropriate weak formulation of the problem is derived for the use of   evolving finite elements designed to accommodate for a moving  interface. Optimal order error bounds are proved for arbitrary order evolving  isoparametric finite elements. The paper concludes with numerical results for a model problem verifying orders of convergence.

§ INTRODUCTION

The model studied in the paper is the following, let $\Omega$ be a stationary domain with a moving interface $\Gamma(t)$ that encloses a subdomain $\Omega_1(t)$ and let $\Omega_2(t) = \Omega \setminus \overline{\Omega}_1(t)$. We are interested in well posedness and a suitable finite element scheme for the solutions $u_1, u_2$, scalar fields on the subdomains $\Omega_1(t)$ and $\Omega_2(t)$,  to:

\begin{align}
\partial_t u_i - \nabla\cdot (\mathcal{A}_i(t;x) \nabla u_i) +\mathcal{B}(t;x) \cdot \nabla u_i + \mathcal{C}(t;x)u_i &= f_i(t;x)&&\quad \text{in} \; \Omega_i(t),\\
u_2 &= 0&&\quad \text{on} \; \partial\Omega,\\
u_1 -u_2 &= 0 &&\quad \text{on} \; \Gamma(t), \\
\mathcal{A}_1(t;x)\frac{\partial}{\partial \nu_{\Gamma}} u_1 \bigg|_{\Gamma(t)}- \mathcal{A}_2(t;x)\frac{\partial}{\partial\nu_{\Gamma}}u_2\bigg|_{\Gamma(t)} &= G(t;x)&&\quad \text{on} \; \Gamma(t),\\
u_i(0) &= u_i^0 && \quad \text{on} \; \Omega_i(0).
\end{align}

Such equations can arise as subproblems when modelling the concentration of a dissolved chemical species transport around two species in space.
In particular, we mention applications in fluid dynamics [39, 11, 1], materials science [10, 24] and cell biology [40, 37, 25].

There are two main difficulties concerning this problem. The first of which is the evolution of the sub-domains and the second is the presence of a discontinuous jump across the interface. One common approach to moving domains is the ALE (Arbitrary Eulerian Lagrangian) method, see [38, 27, 30]. This involves having a parametrisation of the evolving region, the flow associated with this parametrisation could be physical or could be made to fit a specific purpose such as in [17]. The flow constructed in this paper only requires the knowledge of the surface velocity and by use of a harmonic extension with some interpolation, a global flow is then extracted to move the nodes of the mesh.
For the interface condition, common methods are either to use a discontinuous or immersed Galerkin method [41, 2, 32]. In this paper we propose  using an evolving fitted mesh allowing the use of isoparametric elements that more accurately approximate the boundary, resulting in higher order error estimates.

The key contributions of this work are:

* We provide a functional analytic setting to show well posedness of the continuous problem (<ref>).

* We provide a new isoparametric finite element method for moving interface problems. Our method for moving the mesh is to move the Lagrange nodes with a given, smooth velocity. The difficulty in achieving a higher order method arises from constructing a good initial mesh.

* We provide a robust error bound which demonstrates the error in an $L^2$ norm is bounded, up to a constant, by $h^{k+1}$, where $h$ represents the mesh size and $k$ is the degree of polynomials used both for the discretisation of the domain and the solution. This is the same order error as if we interpolated a known smooth solution.

* Numerical results and the simulation code are provided both to demonstrate the results and to allow others to use the implementation.

The assumption is made that we are given a global, smooth velocity field $\mathbf{w}$.  Furthermore, said velocity field is such that it preserves the regularity of the mesh over time. The velocity may be derived from physical considerations or otherwise an arbitrary velocity constructed in order to define a well behaved numerical scheme. We do not address how to achieve such a velocity in this work.
There are methods in the literature to prevent mesh deformation, which involve re-parametrising the flow responsible for the movement of the interface into a more suitable flow, see, for example, [9, 18, 19].

§.§ Outline

Section <ref> gives a well posedness analysis of the continuous equations along with the necessary functional analysis setting.
The finite element construction is in Section <ref> and the finite element scheme is in Section <ref>.
An optimal order error bound is shown in Section <ref> under smoothness assumptions on the domain and its evolution and the solution.
Section <ref> includes a time discretisation of the finite element scheme along with numerical experiences demonstrating the error bounds are tight.
The Appendix includes further details of the proof the well posedness of the continuous scheme.

§ EVOLVING SPACE FORMULATION AND WELL POSEDNESS

§.§ Evolving Sobolev Spaces

We set up the necessary tools from the theory of evolving Sobolev spaces which were introduced and developed in [5, 3]. Let $I = [0,T]$ be a closed time interval and let $\{X(t)\}_{t \in I}$ be a family of Banach spaces equipped with norm $\|\cdot\|_{X(t)}$. Assume that there exists a linear map $\phi_t:X(0) \to X(t)$ that satisfies the following properties:

The map $\phi_t$ is invertible for all $t \in I$ with inverse denoted by $\phi_{-t}$ and $\phi_0$ is the identity.
There exists a constant $C$ independent of time such that $\norm{\phi_t \eta}_{X(t)} \leq C\norm{\eta}_{X(0)}$, $\norm{\phi_{-t} \tilde{\eta}}_{X(0)} \leq C\norm{\tilde{\eta}}_{X(t)}$, for all $\eta \in X(0)$ and $\tilde{\eta} \in X(t)$, for all $t \in I$.
The map $t \mapsto \norm{\phi_t \eta}_{X(t)}$ is measurable  for all $\eta \in X(0)$.

Then if such a map exists, we call it the flow map and the pair $(X(t), \phi_t)|_{t \in I}$ a compatible pair. Given a compatible pair, define the Banach moving spaces as:

\begin{align*}
L^p_{X} :&= \bigg \{ \eta: I \to \bigcup_{t \in I} X(t) \times \{t\}, \quad t \to (\hat{\eta}(t),t)|\, \phi_{-t}\hat{\eta}(t) \in L^p(I; X(0)) \bigg\}.
\end{align*}

We identify $\eta(t) = (\hat{\eta}(t), t)$ with $\hat{\eta}(t)$ and equip $L^p_X$ with norm:

\[
\norm{\eta}_{L^p_X} :=
\begin{cases}
\left({\int_0^T \norm{\eta(t)}_{X(t)}^p}\right)^{\frac{1}{p}} &\text{for } p \in [1,\infty),\\
\esssup_{t \in [0,T]} \norm{\eta(t)}_{X(t)} &\text{for } p=\infty.
\end{cases}
\]

$L^p_X$ is indeed a Banach space, see Theorem 3.4 in [3].
The analogues of the spaces of continuous functions and of compactly supported smooth functions are defined as:

\begin{align*}
C^k_X &:= \bigg \{ \eta : I \to \bigcup_{t \in I} X(t) \times \{t\}, \quad t \to (\hat{\eta}(t),t)|\, \phi_{-t}\hat{\eta}(t) \in C^k(I; X(0)) \bigg\}, \\
\mathcal{D}_X &:= \bigg \{ \eta : I \to \bigcup_{t \in I} X(t) \times \{t\}, \quad t \to (\hat{\eta}(t),t)|\, \phi_{-t}\hat{\eta}(t) \in \mathcal{D}(I; X(0)) \bigg\}.
\end{align*}

The strong material derivative in the evolving Banach space setting is defined as follows:

\begin{align*}
\partial_t^\bullet \eta := \phi_t \partial_t \phi_{-t} \eta, \quad \eta \in C^1_X.
\end{align*}

Now assume $\{X(t)\}_{t \in I}$, $\{Y(t)\}_{t \in I}$ and $\{X^*(t)\}_{t \in I}$ are families of Banach spaces, with $Y(t)$ being Hilbert and $X^*(t)$ the dual of $X(t)$ for all $t \in I$. Assume further that for all $t \in I$,  $X(t) \subset Y(t) \subset X^*(t)$ constitutes a Banach triple. It is also assumed that there exists a map $\phi_t:Y(0) \to Y(t)$ with $\phi_{t}|_{X(0)}: X(0) \to X(t)$ with adjoint flow $\phi_{-t}^*:X^*(0) \to X^*(t)$,

\begin{align*}
\langle \phi_{-t}^* f, v \rangle_{X(t)} := \langle f, \phi_{-t} v\rangle_{X(0)}, \quad f \in X^*(0), \, v \in X(t),
\end{align*}

such that $(X(t), \phi_t|_{X(t)})|_{t \in I}$, $(Y(t), \phi_{t})|_{t \in I}$ and $(X^*(t), \phi_{-t}^*)|_{t \in I}$ all define compatible pairs and therefore we can define the spaces $L^p_X, \; L^2_Y, \; L^{p'}_{X^*}$ with their respective flows. In this case, just as for Bochner spaces, we have the have that $(L^p_X)^*$ is isometrically isomorphic to $L^{p'}_{X^*}$ where $p'$ is the Holder complement of $p$ (Theorem 3.7 in [3]). Moreover, the Banach triple structure is preserved: $L^p_X \subset L^2_Y \subset L^{p'}_{X^*}$ for $p \in [2,\infty]$. Note that $L^2_Y$ remains a Banach space with a natural inner product structure, see Remark 3.9 [3]. In order to generalise the concept of a “weak time derivative” to the evolving space, we first assume the following:

The map $t \mapsto \langle \phi_t w_0, \phi_t v_0 \rangle_{X(t)} = (\phi_t w_0, \phi_t v_0)_{Y(t)}$ is continuously differentiable for fixed $w_0, v_0 \in X_0$.
For all $t \in I$, the map:

\begin{align*}
(w_0, v_0) \mapsto \frac{\partial}{\partial t} (\phi_t w_0, \phi_t v_0)_{Y(t)},
\end{align*}

for $(w_0, v_0) \in X(0) \times X(0)$ is continuous
There exists a constant $C_\lambda$ independent of time such that, for almost all $t \in I$ and $w_0, v_0 \in X(0)$, we have:

\begin{align*}
\bigg| \frac{\partial}{\partial t} (\phi_t w_0, \phi_t v_0)_{Y(t)} \bigg| \leq C_\lambda \norm{w_0}_{Y(0)} \norm{v_0}_{Y(0)}.
\end{align*}

Let assumption <ref>–<ref> hold and label:

\begin{align*}
\lambda(t; w,v):= \bigg[\frac{\partial}{\partial t} (\phi_t w_0, \phi_t v_0)_{Y(t)}\bigg]\bigg|_{(w_0,v_0) = (\phi_{-t} w, \phi_{-t} v)}, \quad w,v \in X(t).
\end{align*}

Then $\lambda(t; \cdot, \cdot): Y(t) \times Y(t) \to \mathbb{R}$ is a continuous, symmetric, bounded and bilinear for almost all $t \in I$. We say $w \in L^1_X$ has a weak material derivative if there exists $v \in L^1_{X^*}$ such that:

\begin{align*}
\int_0^T ( w(t), \partial ^\bullet_t \eta)_{Y(t)}\, dt = \int_0^T \langle v(t), \eta \rangle_{X(t)} + \lambda(t; w, \eta)\,dt,
\end{align*}

for all $\eta \in \mathcal{D}_X$. We label $v = \partial_t^\bullet w$.
Note that this can be further generalised to the non Banach triple setting (see Section 4.2 in [3]).
This definition follows all properties we expect from a weak derivative, such as being equivalent to the strong material derivative if the function is regular enough.

This allows us to define the equivalent of the Bochner solution space.
We define the solution space $W^{p,p'}(X,X^*) := \{v \in L^p_X, \; \partial_t^\bullet v \in L^{p'}_{X^*}\}$ with norm:

\begin{align*}
\norm{v}_{W}: = \norm{v}_{L^p_X} + \norm{\partial_t^\bullet v}_{L^{p'}_{X^*}}.
\end{align*}

In the case where $p =p' = 2$, we label $W^{2,2}(X,X^*) = W(X,X^*)$ and use the equivalent norm:

\begin{align*}
\norm{v}^2_{W}: = \norm{v}^2_{L^2_X} + \norm{\partial_t^\bullet v}^2_{L^{2}_{X^*}}.
\end{align*}

The last assumption needed is the moving space equivalence:
The solution space is said to satisfy the moving space equivalence if:

\begin{align*}
v \in W^{p,p'}(X,X^*) \iff \phi_{-t} v \in \mathcal{W}^{p,p'}(X(0), X^*(0)),
\end{align*}

where $\mathcal{W}^{p,p'}(X(0),X^*(0)) = \{v_0 \in L^p(I;X(0)), \; \partial_t v \in L^{p'}(I; X^*(0))\}.$
(The Transport Theorem) Assume $v,w \in W^{p, p'}(X,X^*)$ and the moving space equivalence is satisfied, then ,the map $t \mapsto (v,w)_{Y(t)}$ is uniformly continuous and for almost all $t \in I$, and the following holds:

\begin{align*}
\frac{d}{dt} (v,w )_{Y(t)} = \langle \partial_t^\bullet v,w\rangle_{X(t)} +\langle \partial_t^\bullet w,v\rangle_{X(t)} + \lambda(t;u,v).
\end{align*}

Moreover, $C_Y \hookrightarrow  W^{p, p'}(X,X^*)$.
See Section 4.5 in [3] for proofs.

§.§ Setting up the Domain

An example configuration of the domain.

Let $\Omega$ be a stationary domain in $\mathbb{R}^d$, $d=2,3$, with piecewise linear boundary and let $\{\Gamma(t), \; t \in I\}$ be a family of closed compact connected $C^{2+k}$ ($k\geq 0$) hypersurfaces with $\Gamma(t) \subset \Omega$. Assume that there is a domain $\Omega_1(t)$ bounded by $\Gamma(t)$ such that $\partial \Omega_1(t) = \Gamma(t)$ for all $t \in I$. Let $\Omega_2(t): = \Omega \setminus \Omega_1(t)$ and assume that $\Gamma(t) \cap \partial \Omega = \varnothing$, then:

\begin{align*}
\overline{\Omega} = \overline{\Omega}_1(t) \cup \overline{\Omega}_2(t), \; \overline{\Omega}_1(t) \cap \overline{\Omega}_2(t) = \Gamma(t), \, \partial\Omega_2(t) = \Gamma(t) \cup \partial \Omega.
\end{align*}

A sketch of the domains is shown in Figure <ref>.

Remark. The assumption that the outer boundary is piecewise linear is made to avoid having to analysis perturbation of the domain for Dirichlet boundary conditions, however the method presented can easily be altered if one removes this assumption.

We label the outer normals of $\Omega_1(t)$ and $\Omega_2(t)$ by $\nu_{\Gamma(t)}$ and $\nu_{\partial \Omega_2(t)}$ respectively.

\begin{align*}
\mathcal{Q}_i := \bigcup_{t \in I} \Omega_i(t) \times \{t\}, \; \mathcal{Q}:= \Omega \times I.
\end{align*}

Furthermore there exists a given velocity field $\mathbf{w}$ transporting $\Omega_1(t)$ and $\Omega_2(t)$, i.e $ \mathbf{w} \cdot \nu_{\Gamma(t)}|_{\Gamma(t)}= V_{\Gamma}$ where $V_\Gamma$ is the normal velocity of $\Gamma(t)$, $\mathbf{w} \cdot \nu_{\partial \Omega_2(t)}|_{\partial \Gamma(t)} = -V_{\Gamma}$ and $\mathbf{w} \cdot \nu_{\partial \Omega_2(t)}|_{\partial \Gamma} = 0$. This velocity is assumed to be of regularity $\mathbf{w} \in C(I; C(\overline{\Omega}; \mathbb{R}^d))$ with $\mathbf{w}_i(t; \cdot) \in C^{2+k}(\overline{\Omega}_i(t); \mathbb{R}^d)$. Let $\mathbf{\Phi}_i(t;x): \overline{\Omega}_i(0) \to \overline{\Omega}_i(t)$ be the solution to:

\begin{align*}
\frac{d}{dt} \mathbf{\Phi}_i(t;x) =\mathbf{w}(t; \mathbf{\Phi}_i(t;x)), \qquad
\mathbf{\Phi}_i(0; x) = x.
\end{align*}

The solution exist and is of regularity $\mathbf{\Phi}_i \in C^{1}(\overline{\mathcal{Q}}_i; \mathbb{R}^d)$ with $\mathbf{\Phi}_i(t;\cdot): \overline{\Omega}_i(0) \to \overline{\Omega}_i(t)$ and $\mathbf{\Phi}_i(t;\cdot) \in C^{2+k}(\Omega_i(0); \mathbb{R}^d)$, see [26]. Furthermore, both $\mathbf{\Phi}_i(t; \cdot)$ are invertible diffeomorphisms for all $t \in I$ with $\mathrm{Im}(\mathbf{\Phi}_i(t;\cdot)) = \Omega_i(t)$. Since we assumed $\mathbf{w} \in C(I; C(\overline{\Omega}; \mathbb{R}^d))$, it follows that $\mathbf{\Phi}_1(t;x)|_{\Gamma(0)} = \mathbf{\Phi}_2(t;x)|_{\Gamma(0)}$.

Let $J^t_i$ denote the determinant of Jacobian matrix, $J^t_i = \text{det} [\nabla \mathbf{\Phi}_i(t;x)]$. The prior assumptions imply $J^t_i(\cdot) \in C^{1+k}(\overline{\Omega}_i(0); \mathbb{R})$ and there exists $C_\Omega$ independent of time and space such that:

\begin{align*}
\frac{1}{C_\Omega} \leq |J_i^t | \leq C_\Omega,
\end{align*}

and $J^{-t}_i$ denotes its inverse.

Remark If instead only the velocity $\mathbf{w}_\Gamma$ on $\Gamma(t)$ such that $\mathbf{w}_\Gamma \cdot \nu_{\Gamma(t)} = V_\Gamma$ was known, a global velocity field $\mathbf{w}$ such that $\mathbf{w}|_{\Gamma(t)} = \mathbf{w}_\Gamma$, $\mathbf{w} \cdot \nu_{\partial \Omega} = 0$ can still be constructed, see, for example, Section 2.5 in [33].

Let $d_\Gamma(t;x)$ be the signed distance function:

\begin{align*}
d_\Gamma(t;x) = \begin{cases}
-\inf\{|x - y|: y \in \Gamma(t)\}, \quad &\text{for} \, x \in \overline{\Omega}_1(t), \\
\inf\{|x - y|: y \in \Gamma(t)\}, \quad &\text{for} \, x \in \Omega_2(t) .
\end{cases}
\end{align*}

Then, since the interface is of class $C^{2 +k}$, there exists a constant  $\delta >0$ such that if $ x \in \mathcal{N}_{\Gamma(t)} := \{x \in \Omega. \; | d_{\Gamma}(t;x)| \leq \delta \}$, then it can be uniquely decomposed as:

\begin{align*}
x = d_{\Gamma}(t;x) \nu_{\Gamma(t)}(\Pi_{t}(x)) + \Pi_t(x),
\end{align*}

where $\Pi_t(\cdot)$ is the nearest point on $\Gamma(t)$, i.e, $\Pi_t(x) := \{y \in \Gamma(t), \; |y - x| =  \inf_{y' \in \Gamma(t)} |x - y'| \}$  (see Section 2.3 [33]). We refer to the set $\mathcal{N}_{\Gamma(t)}$ as the tubular neighbourhood of $\Gamma(t)$. Moreover, since the interface is $C^{2+k}$,  $d_{\Gamma}(t; \cdot) \in C^{2+k}(\mathcal{N}_{\Gamma(t)}; \mathbb{R})$ and $\Pi_t \in C^{1 +k}(\mathcal{N}_{\Gamma(t)}; \mathbb{R}^d)$ (see [21] and Lemma 14.16 in [23]). Note that $\delta$ can be chosen independent of time and, by the assumption that $\Gamma(t) \cap \partial \Omega = 0$ for all $t \in I$, such that $\mathcal{N}_{\Gamma(t)} \cap \partial \Omega = \varnothing$ for all $t \in I$.

§.§ Realisation

For a given function $v$ acting on $\Omega$, we decompose it as:

\begin{align*}
v_1  := v \chi_{\overline{\Omega}_1(t)}, \quad v_2 :=  v \chi_{\overline{\Omega}_2(t)},
\end{align*}

where $\chi_{\overline{\Omega}_i(t)} = 1$ if $x \in \overline{\Omega}_i(t)$ and zero otherwise. The jump of the function $v$ is defined as follows:

\begin{align*}
\llbracket v \rrbracket  = [v_1 - v_2]|_{\Gamma(t)}.
\end{align*}

A function $v$ on $\Omega$ will be decomposed as a pair $v = (v_1, v_2)$. Define the following spaces:

\begin{align*}
H(t) &= L^2(\Omega_1(t)) \times L^2(\Omega_2(t)), \quad \norm{v}_{H(t)}^2 = \sum_{i = 1}^2 \norm{v_i}^2_{L^2(\Omega_i(t))},  \\
V(t) &= \{v \in H^1(\Omega_1(t)) \times H^1(\Omega_2(t)),\, \llbracket v \rrbracket = 0, \; v_2|_{\partial \Omega} = 0\}, \quad \norm{v}_{V(t)}^2 = \sum_{i = 1}^2 \norm{v_i}^2_{H^1(\Omega_i(t))}.
\end{align*}

Note that due to the continuity of the trace operators $T_i(t)v := v|_{\partial \Omega_i(t)}$ on $H^1(\Omega_i(t))$, $V(t)$ defines a closed subset of $H^1(\Omega_1(t)) \times H^1(\Omega_2(t))$ and contains $H_0^1(\Omega_1(t)) \times H_0^1(\Omega_2(t))$, hence is dense within $H(t)$. For a $v = (v_1, v_2) \in V(t)$, we will identify $v|_{\Omega_1(t)} = v_1$ and $v|_{\Omega_2(t)} = v_2$.
We also define the interface space:

\begin{align*}
H^{1/2}(\Gamma(t)) &= \{ v \in L^2(\Gamma(t)), \, |v|_{H^{1/2}(\Gamma(t))} < \infty\}, \; |v|_{H^{1/2}(\Gamma(t))} := \int_{\Gamma(t)}\int_{\Gamma(t)} \frac{|v(x) - v(y)|^2}{|x - y|^d }  \, dxdy,  \\
\intertext{with norm given by:}
\norm{v}^2_{H^{1/2}(\Gamma(t))}:&= \norm{v}^2_{L^2(\Gamma(t))} + |v|^2_{H^{1/2}(\Gamma(t))}.
\end{align*}

Then, $H^{1/2}(\Gamma(t))$ is a Hilbert space and moreover is dense and compactly embedded in $L^2(\Gamma(t))$ (see Section 2 in [31]). For consistency of notation, let $\mathcal{V}_\Gamma(t) = H^{1/2}(\Gamma(t))$ and $\mathcal{H}_\Gamma(t) = L^2(\Gamma(t))$, and identify the Banach triple $\mathcal{V}_\Gamma(t) \subset \mathcal{H}_\Gamma(t) \subset \mathcal{V}_\Gamma^*(t)$.

Now for a function $v \in H(t)$ and $w \in \mathcal{H}_\Gamma(t)$, the respective flows are defined as:

\begin{align*}
\phi_t v := (v_1(t;\mathbf{\Phi}_1(-t;x)), v_2(t;\mathbf{\Phi}_2(-t;x))), \quad \phi_t w := w(t;\mathbf{\Phi}_1(-t;x)).
\end{align*}

The pairs $(V(t), \phi_t)|_{t \in I}$, $(H(t), \phi_t)|_{t \in I}$, $(V^*(t), \phi_{-t}^*)|_{t \in I}$, $(\mathcal{V}_\Gamma(t), \phi_t)|_{t \in I}$, $(\mathcal{H}_\Gamma(t), \phi_t)|_{t \in I}$ and $(\mathcal{V}_\Gamma^*(t), \phi_{-t}^*)|_{t \in I}$ are all compatible.

Assumptions <ref>–<ref> need to be checked. This will be checked only for $V(t)$ as a similar logic can be employed for the remaining spaces.
<ref> follows from both $\mathbf{\Phi}_1(-t;x)$, $\mathbf{\Phi}_2(-t;x)$ being invertible diffeomorphisms.
For <ref>, via simple manipulation:

\begin{align}\label{compatibility}
\norm{\phi_t v}^2_{V(t)}  &= \sum_{i = 1}^2 \int_{\Omega_i(t)} |v_i(t;\mathbf{\Phi}_i(-t;x))|^2 + |\nabla v_i(t;\mathbf{\Phi}_i(-t;x))|^2, \nonumber \\
&= \sum_{i = 1}^2 \int_{\Omega_i(0)}[ \,|v_i(t;x)|^2 + |[\nabla \mathbf{\Phi}_i(-t;y)]^T|_{y = \mathbf{\Phi}(t;x)} \nabla v_i(t;x)|^2]J^t_i, \\
&\leq C(|J^t_i|_{L^\infty(\Omega_i(0))}, |\nabla \mathbf{\Phi}_i(-t,x)|_{L^\infty(\Omega_i(t))}) \norm{v}^2_{V(0)}, \nonumber
\end{align}

The bound follows from the assumption on the regularity of the velocity field. The same method shows a similar bound for $\norm{\phi_{-t} \tilde{v}}_{V(0)} \leq C\norm{\tilde{v}}_{V(t)}$ for all $\tilde{v} \in V(t)$. To show measurability, <ref>, note that the second equality in (<ref>) is continuous.

Identifying both $V(t) \subset H(t) \subset V^*(t)$ and $\mathcal{V}_\Gamma(t) \subset \mathcal{H}_\Gamma(t) \subset \mathcal{V}_\Gamma^*(t)$ with the structure $X(t) \subset Y(t) \subset X^*(t)$ developed in Section <ref>.
The moving space equivalence is satisfied between $W(V,V^*)$ and $\mathcal{W}(V(0), V^*(0))$.

The proof follows identically from the one given in Proposition 7.4 [3] as by assumption the Jacobian determinate is at least of regularity $J^t_i \in C^1 (I; C^1(\overline{\Omega}_i(0); \mathbb{R}))$.

Remark. It does not matter which flow $\mathbf{\Phi}_i(t;x)$ is used to define $\mathcal{H}_\Gamma(t)$ as $\mathbf{\Phi}_1(t;x)|_{\Gamma(0)} = \mathbf{\Phi}_2(t;x)|_{\Gamma(0)}$. Moreover, it can be shown that $v = (v_1, v_2) \in V(t)$ if, and only if, $v_1\chi_{\Omega_1(t)} + v_2(1 - \chi_{\Omega_1(t)}) \in H^1_0(\Omega)$ with equivalent norms, hence the space $V(t)$ can be thought as an identification of the components of a functions in $H^1_0(\Omega)$.

Hence we may define both moving space triples $L^2_V \subset L^2_H \subset L^2_{V^*}$ and $L^2_{\mathcal{V}_\Gamma} \subset L^2_{\mathcal{H}_\Gamma} \subset L^2_{\mathcal{V}^*_\Gamma}$.
Let $g_i \in C^1(\Omega_i(t); \mathbb{R})$, then:

\begin{align*}
\frac{d}{dt} \int_{\Omega_i(t)} g_i = \int_{\Omega_i(t)} \partial_t g_i +  \mathbf{w} \cdot \nabla g_i  + g_i \nabla \cdot \mathbf{w}.
\end{align*}

Hence for a function $v \in C^1_V$, we define:

\begin{align*}
\partial_t^\bullet v = \phi_t \frac{d}{dt}(v_1(t;\mathbf{\Phi}_1(-t;x)), v_2(t;\mathbf{\Phi}_2(-t;x))) = ([\partial_t +\mathbf{w} \cdot \nabla]v_1(t;x), [\partial_t +\mathbf{w} \cdot \nabla]v_2(t;x)) =: (\partial_t^\bullet v_1, \partial_t^\bullet v_2).
\end{align*}

One can check that due to the regularity of the flow, the assumptions of <ref>–<ref> are satisfied on the triple $V(t) \subset H(t) \subset V^*(t)$, moreover, via Reynold's transport theorem, one can check that the bilinear form $\lambda$ introduced Section <ref> in this case becomes:

\begin{align*}
\lambda(t;v, \eta) = (\nabla \cdot \mathbf{w} v, \eta)_{H(t)} = \sum_{i = 1}^2 \int_{\Omega_i(t)}[\nabla \cdot \mathbf{w} ]v_i\,\eta_i .
\end{align*}

§.§ The Weak Formulation

Taking the strong problem (<ref>), assuming there exists a regular enough solution $u$, we can rewrite the partial differential equation as:

\begin{align}
\partial_t^{\bullet} u_i- \nabla \cdot (\mathcal{A}_i (t;x) \nabla u_i)  +[\mathcal{B}(t;x)- \mathbf{w}] \cdot \nabla u_i +\nabla \cdot \mathbf{w}u_i+[\mathcal{C}(t;x) - \nabla \cdot \mathbf{w}]u_i=  f_i.
\end{align}

Where the term $\nabla \cdot \mathbf{w}$, corresponding to the previously identified bilinear form $\lambda(t; \cdot,\cdot)$ is introduced to get the equation in a more convenient form. Then testing with a function $v \in L^2_V$ and using the interface condition, we arrive at the following variational problem:

\begin{align*}
\int_0^T \langle \partial_t^{\bullet} u, v \rangle_{V(t)}\,dt +\int_0^T \underbrace{\sum_{i = 1}^2\int_{\Omega_i(t)} \mathcal{A}_i \nabla u_i \cdot \nabla v_i + [\mathcal{B} - \mathbf{w}] \cdot \nabla u_i \,v_i + [\mathcal{C} - \nabla\cdot\mathbf{w}]u_i\, v_i}_{=:a(t;u,v)}\,dt  +\lambda(t;u,v),\\
= \int_0^T \underbrace{\langle f, v\rangle_{V(t)} +\langle G,v \rangle_{H^{1/2}(\Gamma(t))}}_{=:l(t;v)}\,dt,
\end{align*}

with initial condition $u(0) = u_0$. This gives us the weak formulation:

\begin{align}\label{variational}
\int_0^T \langle \partial_t^\bullet u, v \rangle_{V(t)} +a(t;u,v) +\lambda(t;u,v) = \int_0^T l(t;v),
\end{align}

for all $v \in L^2_V$. Moreover if, instead $v \in W(V,V^*)$, we get the equivalent formulation via the transport theorem (for notational convenience later on, we will label the inner product $(\cdot, \cdot)_{H(t)} =: m(t;\cdot,\cdot)$):

\begin{align*}
\frac{d}{dt} m(t; u, v) + a(t; u, v) = m(t; u, \partial^\bullet_t v) + l(t; v).
\end{align*}

If $\partial_t^\bullet u \in L^2_H$, then via identification of the Banach triple, we have:

\begin{align*}
\langle \partial_t^\bullet u ,v \rangle_{V(t)} = m(t; \partial_t^\bullet u, v),
\end{align*}

and the problem can be restated abstractly in this case as $u \in W(V,H)$ being the solution to:

\begin{align}\label{eq:abstract}
m(t; \partial_t^\bullet u, v) + a(t;u,v) + \lambda(t; u,v)= l(t;v),
\end{align}

for almost all $t \in I$, and all $v \in L^2_V$.

§.§ Well Posedness

Assume the following:

$A1$ the coefficients $\mathcal{A}_i \in C(\overline{\mathcal{Q}}_i; \mathbb{R})$, $\mathcal{B} \in C(I \times \overline{\Omega}; \mathbb{R}^d)$ and $\mathcal{C} \in C(I \times \overline{\Omega}; \mathbb{R})$;
$A2$ there exists a constant $\gamma >0$ such that:

\begin{align*}
\inf_{t \in I} \inf_{x \in \Omega_i(t)} \mathcal{A}_i(t;x) \geq \gamma;
\end{align*}

$A3$ $(u_0, f, G, \mathbf{w}_i, \mathbf{w}) \in H(0) \times L^2_{V^*} \times L^2_{V_\Gamma^*} \times C^1(\overline{\mathcal{Q}}_i, \mathbb{R}^d) \times C(I \times \overline{\Omega}; \mathbb{R}^d))$,

then there exists a unique solution $u \in W(V,V^*)$ with inequality:

\begin{align*}
\norm{u}_{W} \leq C \bigg(\norm{f}_{L^2_{V^*}} + \norm{G}_{L^2_{\mathcal{V}_\Gamma^*}} +\norm{u_0}_{H(0)}\bigg).
\end{align*}

Furthermore, if it holds that:

$A4$ $(u_0,f,G, \mathcal{A}_i) \in V(0) \times  L^2_H  \times W(\mathcal{V}_\Gamma, \mathcal{V}^*_\Gamma) \times C^1(\overline{\mathcal{Q}_i}; \mathbb{R})$,

then the solution is of additional regularity $u \in W(V,H)$ with bound:

\begin{align*}
\norm{u}_{W(V,H)} \leq C\bigg(\norm{f}_{L^2_H} +\norm{u_0}_{V(0)} +\norm{G}_{W(\mathcal{V}_\Gamma, \mathcal{V}^*_\Gamma)} \bigg).
\end{align*}

The existence and uniqueness follows from a standard application of the Babuska-Lax-Milgram theorem in conjunction with Poincaré's inequality, detailed in Theorem 3.6 in [4], the proof of additional regularity under assumption $A4$ is given in the appendix.

Furthermore, it is convenient to define the derivative of the bilinear form $a(t;\cdot, \cdot)$ as:

\begin{align}\label{eq;derive}
b(t;v,w) := \partial_t [a(t;v,w)] - a(t; \partial^\bullet_t v,w) - a(t;v, \partial^\bullet_t w) \quad \forall v,w \in W(V,V).
\end{align}

Then, assuming furthermore that $\mathcal{A}_i \in C^1(\mathcal{Q}_i; \mathbb{R})$, $\mathcal{B} \in C^1(I \times \Omega; \mathbb{R}^d)$ and $\mathcal{C} \in C^1(I \times \Omega; \mathbb{R})$, the bilinear form $b(t; \cdot, \cdot)$ exists and can be explicitly calculated as:

\begin{align}
b(t;v,w) &= \sum_{i = 1}^2 \int_{\Omega_i(t)} \mathcal{D}^\mathcal{A}_i(\mathbf{w}, \mathcal{A}_i, v_i,w_i) + \mathcal{D}^\mathcal{B}(\mathbf{w},B, v_i,w_i) \\
& \qquad\qquad\qquad\qquad +  v_i\,w_i  \partial^\bullet_t [\mathcal{C} - \nabla \cdot \mathbf{w}] +  \nabla \cdot \mathbf{w}[\mathcal{C} - \nabla \cdot \mathbf{w}] v_i\,w_i, \nonumber\\
\mathcal{D}^\mathcal{A}_i(\mathbf{w}, \mathcal{A}_i, v_i, w_i) &= (\partial^\bullet_t \mathcal{A}_i(t;x)+ \nabla \cdot \mathbf{w}\mathcal{A}_i(t;x))\nabla v_i \cdot \nabla w_i  - D_i(\mathbf{w}, v_i,w_i),\nonumber\\
\mathcal{D}^\mathcal{B}(\mathbf{w},B,v_i,w_i) &= \partial^\bullet_t [\mathcal{B}(t;x) - \mathbf{w}] \cdot \nabla v_i\, w_i + [\mathcal{B}(t;x) - \mathbf{w}]\cdot\nabla v_i\, w_i \nabla \cdot \mathbf{w} - \sum_{j,k = 1}^{d} [\mathcal{B} - \mathbf{w}]_j  (\nabla_j \mathbf{w}_k)\nabla_k v_i\,w_i,\nonumber \\
D_i(\mathbf{w}, v_i,w_i) &= \sum_{l,k = 1}^d \mathcal{A}_i(t;x)[\nabla_l \mathbf{w}_k +\nabla_k \mathbf{w}_l]\nabla_l v_i \nabla_k w_i. \nonumber
\end{align}

Note that the derivative of the bilinear form $m(t;\cdot,\cdot)$ is already assumed to exist and equals $\lambda(t;\cdot,\cdot)$ introduced in Section <ref>.

§ EVOLVING FINITE ELEMENTS

§.§ Construction of the Domain

We detail the initial triangulation of the domain:

§.§.§ Initial Mesh Construction/Assumption:

$M1$  We first perform a partition into $d$-dimensional simplices corresponding to a polyhedral approximation of the interior domain $\tilde{\Omega}^h_1 \subset \overline{\Omega}_1(0)$, $\tilde{\Omega}^h_1 = \cup_{j = 1}^{M_1} \tilde{K}^j_1$, $\tilde{\mathcal{J}}^h_1:= \{\tilde{K}^j_1\}_{j = 1}^{M_1}$, where $\tilde{K}^j_1$ are the triangular elements of maximum diameter $h >0$ bounded by some $\tilde{h}$ and $M_1$ the number of elements.
$M2$  Let $\tilde{\Omega}^h_2:= \overline{\Omega \setminus \tilde{\Omega}^h_1}$, and assume that $\tilde{h}$ is small enough such that there exists an exact triangular partition $\tilde{\mathcal{J}}^h_2$ such that $\tilde{\Omega}^h_2 = \cup_{j = 1}^{M_2} \tilde{K}^j_2$, $\tilde{\mathcal{J}}^h_2 = \{\tilde{K}^j_2\}_{j = 1}^{M_2}$ with maximum diameter $2\tilde{h}$. Let $\tilde{\mathcal{J}}^h = \cup_{i = 1}^2 \tilde{\mathcal{J}}^h_i$ assume that all partitions $\{\tilde{\mathcal{J}}_1^h, \tilde{\mathcal{J}}_2^h, \tilde{\mathcal{J}}^h\}$ are admissible, shape regular and uniform in $\{\tilde{\Omega}_1^h, \tilde{\Omega}^h_2, \overline{\Omega}\}$ respectively, see [12] Definition 5.1.
$M3$  Each element $\tilde{K}$ contains $d+1$ facets labelled $\{\tilde{E}^j\}_{j = 1}^{d+1} \subset \tilde{K}$. We refer to the set of all facets of all elements in $\tilde{\mathcal{J}}$ by $\tilde{\mathcal{G}}$.
$M4$  For $\tilde{E} \in \tilde{\mathcal{G}}$, if there exists $\tilde{K}_1 \in \tilde{\mathcal{J}}_1^h$ and $\tilde{K}_2 \in \tilde{\mathcal{J}}^h_2$ such that $\tilde{E} = \tilde{K}_1 \cap  \tilde{K}_2$, then we call $\tilde{E}$ an interface facet and label the collection of those facets $\tilde{\mathcal{G}}_\Gamma^h$ and the union of interface facets $\tilde{\Gamma}^h_0$. If for a given $\tilde{E}$, there is only one element $\tilde{K}_i \in \tilde{\mathcal{J}}^h_i$ such that $\tilde{E} \subset \tilde{K}_i$, then such a facet is called a boundary facet and label the set of boundary facets $\partial \tilde{\Omega}^h_i(0)$, $i = 1,2$.
$M5$  We restrict the vertices of interface facets to be on $\Gamma_0$, i.e, if $\tilde{E}$ is an interface facet, and $\{\tilde{a}^K_j\}_{j = 1}^{d-1}$ are the vertices of $\tilde{E}$, then $\{\tilde{a}^K_j\}_{j = 1}^{d-1} \subset \tilde{E} \cap \Gamma_0$. Conversely, we will assume that if a facet has all its vertices on the interface, then it is a boundary facet. See Figure <ref> for an example.
$M6$  Let $\hat{K}_{ref}$ be the reference element of $\tilde{\mathcal{J}}$ (i.e for all $\tilde{K} \in \tilde{\mathcal{J}}$, there exists an invertible diffeomorphism $F_{\tilde{K}}$ such that $F_{\tilde{K}}(\hat{K}_{ref}) = \tilde{K}$). The reference element is then equipped with the standard $k^{th}$ Lagrangian element triple $(\hat{K}_{ref}, \hat{P}^k, \hat{\Sigma}^k)$ (see [13] section 3.2) where $\hat{P}^k$ is the set of $k^{th}$ order Lagrange polynomials and $\hat{\Sigma}^k$ is the dual basis of $\hat{P}^k$, which in this case takes the form $\hat{\Sigma}= \{ \chi \to \chi(\hat{\alpha}), \; \hat{\alpha} \in N(\hat{K}_{ref})\}$, where $N(\hat{K}_{ref})$ is the set of Lagrangian nodes in $\hat{K}_{ref}$. Let $(\tilde{K}, \tilde{P}^k,\tilde{\Sigma}^k)$ and $(\tilde{K}', \tilde{P}'^k, \tilde{\Sigma}'^k)$ be two adjacent elements in $\mathcal{\tilde{J}}$, the following assumption is made

\begin{align*}
\bigg(\bigcup_{\alpha \in N(\tilde{K})} \tilde{\alpha} \bigg) \cap \tilde{K}' = \bigg(\bigcup_{\alpha' \in N(\tilde{K}')} \tilde{\alpha}' \bigg) \cap \tilde{K},
\end{align*}

i.e the Lagrangian nodes are shared between two elements.

Showing the difference between two interior meshes, one left one breaking <ref> whereas the right one following <ref>.

Note that via construction, $\tilde{E} \in \tilde{\mathcal{G}}_\Gamma^h$, if and only if there exits an element $\tilde{K}_1 \in \tilde{\mathcal{J}}^h_1$ and $\tilde{K}_2 \in\tilde{\mathcal{J}}^h_2$ with $\tilde{E} = \tilde{K}_1 \cap \tilde{K}_2$ and hence $\tilde{\Gamma}^h_0 = \tilde{\Omega}^h_1(0) \cap \tilde{\Omega}^h_2(0)$. Note that this construction defines Lagrangian triangulated bulk domains $(\tilde{\Omega}^h_1, \tilde{\Omega}^h_2, \tilde{\Omega}^h)$, and $\tilde{\Gamma}^h_0$ defines a triangulated hypersurface, see Definitions 4.14 and 6.14 in [20].

After the initial triangulation, we define the isoparametric version using the same method as in [20] Section 8.5. Let $\hat{K}_{ref}$ be the reference element of the partition $\tilde{\mathcal{J}}^h$, with reference map $F_{\tilde{K}}:\hat{K}_{ref} \to \tilde{K}$. For $\eta \in C(\overline{\Omega}; \mathbb{R})$ and for some $\tilde{K} \in \tilde{\mathcal{J}}^h$, we define the interpolation operator:

\begin{align*}
\tilde{I}^k(\eta)|_{\tilde{K}} = \sum_{j = 1}^{N(k)} \sigma^j_k(\eta) \chi_j = \sum_{j = 1}^{N(k)} \eta(\alpha_j) \chi_j,
\end{align*}

where $N(k) = |N(\tilde{K})|$, the number of elements in $\hat{\Sigma}^k$. Let $\{\tilde{a}^{\tilde{K}}_j\}_{j  =1}^{d+1}$ be the vertices of an element $\tilde{K} \in \tilde{\mathcal{J}}^h$. If two or more of said vertices are on the interface $\Gamma_0$, then the element is referred to as an interface element. Let $\tilde{\mathcal{F}}$ be the set of all interface elements and define the following function $\mathbf{\Psi}^h: \Omega \to \Omega$ element-wise as follows. If $\tilde{K} \notin \tilde{\mathcal{F}}$, then $\mathbf{\Psi}^h(x) = x$ for $x \in \tilde{K}$. If instead $\tilde{K} \in \tilde{\mathcal{F}}$, then expand $x \in \tilde{K}$ into barycentric coordinates:

\begin{align*}
x = \sum_{j = 1}^{d+1} \mu_j(x) \tilde{a}^{\tilde{K}}_j.
\end{align*}

Let $L_K$ be the number of vertices in $\tilde{K}$ that lie on $\Gamma_0$ ($L_K \geq 2$ by assumptions) and assume that the vertices are ordered so that the first $L_K$ lie on $\Gamma_0$. Let:

\begin{align*}
\tilde{\mu}_K(x):= \sum_{j = 1}^{L_K} \mu_j(x), \quad \sigma_K := \{ x \in \tilde{K}, \; \tilde{\mu}_K(x) = 0 \}.
\end{align*}

From the properties of barycentric coordinates, $\tilde{\mu}_K$ can be seen as the distance from the interface interpolant, with $\tilde{\mu}_K(x) = 1$ when $x$ is on a facet between vertices on the interface, and $\tilde{\mu}_K(x) = 0$ when $x$ is on facet between non-interface vertices.

\begin{align*}
y(x) = \sum_{j = 1}^{L_K} \frac{\mu_j(x)}{\tilde{\mu}_K(x)} \tilde{a}^{\tilde{K}}_j.
\end{align*}

Note that $y(x) \in \tilde{K}$ since $0 \leq \mu_j(x) \leq \tilde{\mu}_K(x)$. Hence define:

\begin{align*}
\mathbf{\Psi}^h|_{\tilde{K}}(x) := \begin{cases}
x + (\tilde{\mu}_K(x))^{k+2}(\Pi_0(y(x)) - y(x)) \quad &\text{if} \, x \notin \sigma, \\
x \quad &\text{otherwise}.
\end{cases}
\end{align*}

We summarise the properties of this map in the following Lemma, for the definition of triangulated bulk domain and k-bulk finite element, see Definitions 4.14 and 4.5 in [20].
For $\tilde{h}$ small enough, the map $ \mathbf{\Psi}^h|_{\tilde{K}} \in C^{k+1}(\tilde{K}; \mathbb{R}^d)$ and is invertible for each $\tilde{K} \in \tilde{\mathcal{J}}^h$ and $ \mathbf{\Psi}^h: \tilde{\Gamma}^h_0 \to \Gamma_0$. Moreover, define the following:

\begin{align*}
F_{K} &:= [\tilde{I}^h  \mathbf{\Psi}^h](F_{\tilde{K}}), \\
K&:= F_K(\tilde{K}), \\
P^k&:= \{ \tilde{\chi}_k \circ F^{-1}_K: \tilde{\chi}_k \in \tilde{P}^k \}, \\
\Sigma^k&:= \{\chi \mapsto \tilde{\sigma}(\chi \circ F_K): \, \tilde{\sigma} \in \tilde{\Sigma}^k\},
\end{align*}

then the triplet $(K, P^k, \Sigma^k)$ with reference map $F_{K}$ defines a $k$-bulk finite element triplet (Definition 4.5 in [20]). Let $\mathcal{J}^h_i = \{[\tilde{I}^h  \mathbf{\Psi}^h](\tilde{K}_i), \; \tilde{K}_i \in \tilde{\mathcal{J}}^h_i\}$, $\mathcal{G}_\Gamma^h = \{ [\tilde{I}^h\mathbf{\Psi}^h]_1(E), \, E \in \tilde{\mathcal{G}}^h_\Gamma\}$ (here $[\tilde{I}^h\mathbf{\Psi}^h]_1$ refers to taking the interpolation with the adjacent element in $\tilde{\mathcal{J}}_1$) and $\mathcal{J}^h = \mathcal{J}^h_1 \cup \mathcal{J}^h_2$, then $\{\mathcal{J}^h _1, \mathcal{J}^h _2, \mathcal{J}^h \}$ are conforming admissible sub-divisions. Furthermore, let:

\begin{align*}
\Omega^h_i(0) := \bigcup_{K_i \in \mathcal{J}^h_i} K_i, \, \Omega^h := \bigcup_{K \in \mathcal{J}^h} K, \, \Gamma^h(0) := \bigcup_{E\in \mathcal{G}_\Gamma^h} E =  \Omega^h_1 \cap  \Omega^h_2.
\end{align*}

Then $(\Omega^h_1(0),\Omega^h_2(0))$ define triangulated bulk domains approximating $(\overline{\Omega}_1(0), \overline{\Omega}_2(0))$, $\Gamma^h(0)$ a triangulated hypersurface approximating $\Gamma(0)$, and $\Omega^h$ defines a triangulated bulk domain that is an exact partition of $\overline{\Omega}$.

Most of the statements follow from Section 8.5 in [20]. What remains to check is whether the interface elements are mapped in such a way to preserve interface facets. For an interface facet $\tilde{E}$ with two adjacent element $\tilde{K}_i$ we require $\mathbf{\Psi}^h|_{\tilde{K}_1}(\tilde{E}) = \mathbf{\Psi}^h|_{\tilde{K}_2}(\tilde{E})$. By construction of the mesh, $x \in \tilde{E}$, $\mu_{K_1}(x) = \mu_{K_2}(x) = 1$, hence both maps $\mathbf{\Psi}^h|_{\tilde{K}_i}(x) = \Pi_0(x)$. Any Lagrangian node $\tilde{\alpha}_i$ on $\tilde{E}$ will be mapped by both maps to $\alpha_i := \Pi_0(\tilde{\alpha}_i)$. Since each interface facet contains the exact amount of nodes to uniquely define a polynomial on the facet, which must equal the restriction on the interface element of the Lagrangian polynomial on the full element (see [12] Remark 5.4), hence for $x \in \tilde{E}$:

\begin{align*}
[I^h \mathbf{\Psi}^h]|_{\tilde{K}_1}(x) &= \sum_{j = 1}^{N(k)} \mathbf{\Psi}^h(\tilde{\alpha}^{\tilde{K}_1}_j) \chi_j(x) = \sum_{\tilde{\alpha}_j^{\tilde{K}_1} \in \tilde{E}} \Pi_0(\tilde{\alpha}_j^{\tilde{K}_1}) \chi_j(x),\\
&= \sum_{\tilde{\alpha}_j^{\tilde{K}_2} \in \tilde{E}} \Pi_0(\tilde{\alpha}_j^{\tilde{K}_2}) \chi^j(x) = \sum_{j = 1}^{N(k)} \mathbf{\Psi}^h(\tilde{\alpha}^{\tilde{K}_2}_j) \chi_j(x) =  [I^h \mathbf{\Psi}^h]|_{\tilde{K}_2}(x).
\end{align*}

For $\tilde{E} \in \tilde{\mathcal{G}}^h_{\Gamma}$, let $\tilde{K}_i \in \tilde{\mathcal{J}}^h_i$ be the adjacent elements and $K_i = I^h \mathbf{\Psi}^h(\tilde{K}_i)$. Since the map $I^h \mathbf{\Psi}^h$ is invertible onto its image for each $\tilde{K} \in \tilde{\mathcal{J}}^h$ and is continuous across the intersection $\tilde{E}$, it holds that $I^h \mathbf{\Psi}^h$ is invertible on $\tilde{K}_1 \cup \tilde{K}_2$, since both elements are closed, hence:

\begin{align*}
E := [I^h \mathbf{\Psi}^h]|_{\tilde{K}_1}(\tilde{E}) = I^h \mathbf{\Psi}^h(\tilde{K}_1 \cap \tilde{K}_2) = I^h \mathbf{\Psi}^h(\tilde{K}_1) \cap I^h \mathbf{\Psi}^h(\tilde{K}_2) = K_1 \cap K_2.
\end{align*}

Therefore the image of an interface facet remains an interface facet.

Figure <ref> gives an example of the initial parametrisation $I^h \mathbf{\Psi}^h$ for two tetrahedra near an interface. The initial decomposition interface elements may have curved edges along the interface but edges away from the interface remain straight. This means that each face of the tetrahedron may be curved.

The intersection of two interface elements (teal and blue respectively) of different domains. The shared interface facet (prescribed in the tile coloured triangle) is then pushed by the map $\mathbf{\Psi}^h$ to become a piece of $\Gamma_0$, the map $\tilde{I}^h \mathbf{\Psi}^h$ maps the original mesh to a curved mesh approximating the interface.

Let $\alpha^K_j|_{j = 1}^{N(k)}$ be the Lagrangian nodes on an element $K \in \mathcal{J}^h$, which by construction are defined as $\alpha^K_j = F_K(\hat{\alpha}^{ref}_j)$, the corresponding interpolation operator for a given function $\eta \in C(\overline{\Omega}; \mathbb{R})$ is given by:

\begin{align*}
[I^h \eta]|_{K} = \sum_{j = 1}^{N(k)} \eta(\alpha^K_j) \chi_j,
\end{align*}

for $\chi_j \in P^k$.

§.§ Time Dependent Mesh

To induce a flow on the mesh such that it “preserves” the shape of the mesh, since the map $[\tilde{I}^h\mathbf{\Psi}^h]$ is element-wise invertible, define the flow:

\begin{align*}
\mathbf{\Phi}^h_i(t; \cdot)|_{K_i(0)} :=
I^h_{K_i(0)} [ \mathbf{\Phi}_i(t; \mathbf{\Psi}^h \circ (\tilde{I}^h \mathbf{\Psi}^h)^{-1}(\cdot))],
\end{align*}

and note that $\mathbf{\Phi}^h_i(0; \cdot)$ is the identity map and Lagrange points move under the smooth flow $\mathbf{\Phi}_i$.

For $\tilde{h}$ small enough, this map is an invertible diffeomorphism on each element $K_0 \in \mathcal{J}^h$. As before, we summarise the construction in the following lemma:
For $\tilde{h}$ small enough, map $\mathbf{\Phi}^h_t|_{\tilde{K}} \in C^{k+1}(\tilde{K}; \mathbb{R}^d)$ and is invertible onto its image for each $\tilde{K} \in \mathcal{J}^h$. Moreover, define the following:

\begin{align*}
F_{K(t)} &:=  \mathbf{\Phi}^h_t(F_{\tilde{K}}), \\
K(t)&:= F_{K(t)}(\tilde{K}), \\
P^k(t)&:= \{ \chi_k \circ F^{-1}_{K(t)}: \chi_k \in P \}, \\
\Sigma^k(t)&:= \{\chi \mapsto \sigma(\chi \circ F_{K(t)}): \, \sigma \in \Sigma^k\},
\end{align*}

then the triplet $(\hat{K}_{ref}, \hat{P}^k, \hat{\Sigma}^k)$ with reference map $ F_{K(t)}$ defines a bulk evolving finite element triplet. Let $\mathcal{J}^h_i(t) = \{\mathbf{\Phi}^h_i(t;\tilde{K}_i), \; \tilde{K}_i \in\mathcal{J}^h_i\}$, $\mathcal{G}_\Gamma^h(t) = \{\mathbf{\Phi}^h_1(t; E),\, E \in \mathcal{G}_\Gamma^h\}$ and $\mathcal{J}^h(t) = \mathcal{J}^h_1(t) \cup \mathcal{J}^h_2(t)$, then $\{\mathcal{J}^h _1(t), \mathcal{J}^h_2(t), \mathcal{J}^h(t) \}$ are evolving conforming admissible sub-divisions (see Definition 4.32 in [20]). Furthermore, let:

\begin{align*}
\Omega^h_i(t) = \bigcup_{K_i(t) \in \mathcal{J}^h_i(t)} K_i(t), \, \Omega^h(t) = \bigcup_{K \in \mathcal{J}^h} K(t), \, \Gamma^h(t) = \bigcup_{E \in \mathcal{G}_\Gamma^h} E= \Omega^h_1(t) \cap \Omega^h_2(t).
\end{align*}

Then $(\Omega^h_1(t),\Omega^h_2(t))$ define triangulated bulk domains approximating $(\overline{\Omega}_1(t), \overline{\Omega}_2(t))$, $\Gamma^h(t)$ a triangulated hypersurface approximating $\Gamma^h(t)$, and $\Omega^h(t)$ defines a triangulated bulk domain that is an exact partition of $\overline{\Omega}$

The proof follows the same way as in Lemma <ref>.

It will be assumed that the mesh remains uniformly quasi-uniform in time, see Definition 4.35 in [20]. Moreover, this allows us to move the Lagrangian nodes via $\alpha^{K_i(t)}_j(t) = \mathbf{\Phi}_i^h(t; \alpha_j^{K_i})$. The Lagrangian interpolation operator, $I^h|_{K(t)}$, is then defined in the canonical way. Moreover for $x \in K_i(0)$:

\begin{align*}
\partial_t  \mathbf{\Phi}_i^h(t; x) = \sum_{\substack{j = 1 \\ \tilde{\chi}_j \in \tilde{P}^k}}^{N(k)}\partial_t  \mathbf{\Phi}_i(t; \mathbf{\Psi}^h(\alpha_j^{\tilde{K}_i})) \tilde{\chi}_j(x)&= \sum_{\substack{j = 1 \\ \tilde{\chi}_j \in \tilde{P}^k}}^{N(k)}\mathbf{w}(t;  \mathbf{\Phi}_i\circ\mathbf{\Psi}^h(\alpha_j^{\tilde{K}_i}))\tilde{\chi}_j(x), \\
&= \sum_{\substack{j = 1 \\ \chi_j \in P^k(t)}}^{N(k)}\mathbf{w}(t; \alpha_j^{K_i(t)}(t))\chi_j(t;\mathbf{\Phi}_i^h(t; x)) =: \mathbf{w}^h(t;  \mathbf{\Phi}_i^h(t; x)),
\end{align*}

where one sees that $\mathbf{w}^h$ is the interpolated velocity with respect to the moving nodes:

\begin{align*}
\mathbf{w}^h(t;\cdot)|_{K_i(t)} = I^h|_{K_i(t)}[\mathbf{w}(t;\cdot)].
\end{align*}

An example of the temporal deformation of an evolving element is shown in Figure <ref>.
Despite interior elements of the initial partition being linear, since the velocity used to displace the elements is a polynomial interpolant of the velocity, the resulting element might not remain linear and can be deformed. However, as stated earlier, we assume that the mesh maintains uniform in time quasi-regularity.

Example of the temporal deformation of an interior element in three space dimensions.

§.§ The Lift

Finally, the last mesh related object we need is the lift map (see Section 8.6 in [20]). Fix $t \in I$, if $K_i(t) \in \mathcal{J}^h_i(t)$ is an interior element, then define the lift $\Lambda^h_i(t;\cdot)$ as:

\begin{align*}
\Lambda^h_i(t;x) = x, \; \text{for} \; x \in K_i(t).
\end{align*}

If instead $K_i(t)$ is an interface element, we first pull-back the reference map to $\hat{x} \in \hat{K}^{ref}$ such that $x = F_{K_i(t)}(\hat{x})$, then decomposing $\hat{x}$ into barycentric coordinates with respect to the vertices $\hat{a}^{ref}_j$ of $\hat{K}^{ref}$, we have:

\begin{align*}
\hat{x} = \sum_{j = 1}^{d+1} \mu_j(\hat{x}) \hat{a}^{ref}_j,
\end{align*}

and once again, let $L_K$ be the number of vertices on the interface and assume the vertices are ordered so that $\hat{a}_j^{ref}$, $j=1, \ldots, L_K$, get mapped on to $\Gamma(t)$, then we introduce the interface distance and the singular set analogously:

\begin{align*}
\tilde{\mu}(\hat{x}) = \sum_{j = 1}^{L_K}\mu_j (\hat{x}), \; \sigma = \{\hat{x} \in \hat{K}^{ref}|\; \tilde{\mu}(\hat{x}) = 0\}.
\end{align*}

The projection is now defined on the reference element:

\begin{align*}
\tilde{y}(\hat{x}) = \sum_{j = 1}^{L_K} \frac{\mu_i(\hat{x})}{\tilde{\mu}(\hat{x})}, \; y(t;x):= F_{K_i}(t; \tilde{y}(\hat{x})).
\end{align*}

Hence the lift operator can now be defined on interface elements as:

\begin{align*}
\Lambda^h_i(t; x)|_{K_i(t)} = \begin{cases}
x + (\tilde{\mu}(\hat{x}))^{k+2}(\Pi_t(y(t;x)) - y(t;x)) \quad &\text{if} \, x \notin \sigma,\\
x \quad &\text{otherwise}.
\end{cases}
\end{align*}

This gives us this following Lemma:
For $h$ small enough, the map $ \Lambda^h_i(t; \cdot)|_{K_i(t)}$ is a $ C^{k+1}(K_i(t); \mathbb{R}^d)$ element-wise diffeomorphism with image $ \Lambda^h_i(t; \Omega^h_i(t)) = \overline{\Omega}_i(t)$. Moreover, define the following:

\begin{align*}
\mathcal{J}_i^l(t) := \{ \Lambda^l_i(t; K_i(t))|\; K_i(t) \in \mathcal{J}^h_i(t) \}, \quad \mathcal{J}^l(t) := \mathcal{J}_1^l(t) \cup \mathcal{J}_2^l(t).
\end{align*}

Then the triplet $(\mathcal{J}_1^l(t), \mathcal{J}_2^l(t), \mathcal{J}^l(t))$ respectively define a uniform $k$-regular  evolving subdivision of $(\overline{\Omega}_1(t), \overline{\Omega}_2(t), \overline{\Omega})$.
This follows from Lemma 8.12 in [20] and the same method as lemma <ref>.

§.§ Finite Element Spaces

The Broken Sobolev spaces are defined as follows:

\begin{align*}
W^{m,p}_T(\mathcal{J}^h_i(t)) = \{&\eta \in L^1(\Omega^h_i(t)),\, \eta|_{K_i(t)} \in W^{m,p}(K_i(t)) \quad \forall K_i(t) \in \mathcal{J}^h_i(t), \\
&\; \eta|_{\partial K_i(t)} = \eta|_{\partial K'_i(t)} \quad \forall K'_i(t) \in \mathcal{J}^h_i(t) \; s.t \; K'_i(t) \cap K_i(t) \neq \emptyset \},
\end{align*}

equipped with the Broken Sobolev norm:

\begin{align*}
\norm{\eta^h}^p_{W^{m,p}(t)} := \sum_{K_i(t) \in \mathcal{J}^h_i(t)} \norm{\eta^h}^p_{W^{m,p}(K_i(t))}.
\end{align*}

Remark. This does in fact define a Banach space see [6], Section 3.

The discrete spaces are then defined as:

\begin{align*}
H^h(t)& :=   L^2(\Omega^h_1(t))  \times  L^2(\Omega^h_2(t)), \\
V^h(t)&:= \{ \eta^h \in W^{1,2}_T(\mathcal{J}^h_1(t)) \times W^{1,2}_T(\mathcal{J}^h_2(t)), \eta_1^h - \eta_2^h|_{\Gamma^h(t)} = 0 ,\; \text{and} \; \eta^h_2|_{\partial \Omega^h} = 0\},
\end{align*}

Let $\alpha(t)$ be a Lagrangian node and $\mathcal{J}_i(\alpha(t))$ be the set of elements in $K_i(t) \in \mathcal{J}_i(t)$ such that $\alpha(t) \in K_i(t)$, and let $\mathcal{N}^h_i(t)$ be the global set of all Lagrangian nodes in $\mathcal{J}_i^h(t)$. Introduce the finite dimensional subspace:

\begin{multline*}
\mathcal{S}_i^h(t)
:= \biggl\{
\chi^h_i = (\chi^h_i)_{K_i(t) \in \mathcal{J}^h_i(t)} \in \prod_{K_i(t) \in \mathcal{J}^h_i(t)}
\{
\hat{\chi} \circ F^{-1}_{K_i(t)} : \, \hat{\chi} \in P_k(\hat{K})
\} : \\
\chi_i^h|_{K}(\alpha(t)) = \chi^h_i|_{K'}(\alpha(t)) \quad \text{for all} \, K_i(t), K_i'(t) \in \mathcal{J}_i(\alpha(t)), \forall \alpha(t) \in \mathcal{N}^h_i(t)
\biggr\}.
\end{multline*}

Combining two copies of the space yields the adequate solution space:

\begin{multline*}
\mathcal{S}^h (t):= \bigl\{\eta^h = (\eta^h_1, \eta^h_2) \in \mathcal{S}^h_1(t) \times \mathcal{S}^h_2(t) |\; \chi_1^h(\alpha(t)) = \chi_2^h(\alpha(t)) \\ \text{for all } \alpha(t) \in \Gamma^h(t) \cap \mathcal{N}^h_1(t) \; \text{and} \; \chi_2^h (\alpha(t)) = 0 \quad \forall \alpha(t) \in \partial \Omega^h(t) \cap \mathcal{N}^h_2(t)\bigr\},
\end{multline*}

and define the $\phi^h_t:H^h(0) \to H^h(t)$ as:

\begin{align*}
\phi^h_t v^h := (v^h_1(\mathbf{\Phi}^h_1(-t;x)), v^h_2(\mathbf{\Phi}^h_2(-t;x)).
\end{align*}

By construction, all pairs $(V^h(t), \phi^h_t)|_{t \in I}, \, (H^h(t), \phi^h_t)|_{t \in I}, \, (\mathcal{S}^h(t), \phi^h_t)|_{t \in I}$ are compatible. Hence the moving spaces $L^2_{V^h}, \, L^2_{H^h}, \, L^2_{\mathcal{S}^h}$ are well defined. Denote the discrete material derivative by:

\begin{align*}
\partial^h_t \eta := \phi^h_{t} \frac{\partial}{\partial t} \phi_{-t}^h \eta,
\end{align*}

for $\eta \in C^1_{V^h}$. The bilinear form $\lambda^h(t; \cdot, \cdot)$ of Definition <ref> associated with that material derivative is:

\begin{align*}
\lambda^h(t;\eta^h,v^h) = (\nabla \cdot \mathbf{w}^h \eta^h, v^h)_{H^h(t)},
\end{align*}

where: $\mathbf{w}^h$ is the previously defined discrete velocity from Section <ref>.

For a function $v^h \in H^h(t)$, the lift is denoted by $(\cdot)^l:H^h(t) \to H(t)$ and defined as follows:

\begin{align*}
v^{h,l} := (v^h_1(t;\Lambda_1^h(t;x)),v^h_2(t;\Lambda_2^h(t;x))).
\end{align*}

Its inverse will be labelled by $(\cdot)^{-l}$, i.e $(v^{h,l})^{-l}  =v^h$. Since $\Lambda^h_i(t;\cdot) \in W^{k+1, \infty}_T(\mathcal{J}^h_i(t); \mathbb{R}^d)$ and invertible, via a similar change of variable method as Lemma <ref>:

\begin{align*}
&c_1 \norm{v^{h,l}}_{H(t)} \leq \norm{v^h}_{H^h(t)} \leq c_2 \norm{v^{h,l}}_{H(t)},\\
&c_1 \norm{v^{h,l}}_{V(t)} \leq \norm{v^h}_{V^h(t)} \leq c_2 \norm{v^{h,l}}_{V(t)} , \quad \text{if} \; v^h \in V^h(t).
\end{align*}

The lifted solution space can now be defined as:

\begin{align*}
\mathcal{S}^l(t) := \{ \chi^{h,l}|\; \chi^h \in \mathcal{S}^h(t) \}.
\end{align*}

The interpolation operator onto $\mathcal{S}^l(t)$, $I^l:C(\Omega) \to \mathcal{S}^l(t) $ can also be defined in a similar way:

\begin{align*}
I^l(\eta)|_{K(t)} := \sum_{j = 1}^{N(k)} \eta(\alpha^l_j(t)) \chi^{h,l}_j.
\end{align*}

where $\{\alpha^l_j(t)\}_{j = 1}^{N(k)}$ are the lifted Lagrangian Nodes.
We equip the lifted spaces with an analogous flow $\mathbf{\Phi}^l_t: \overline{\Omega}_i(0) \to \overline{\Omega}_i(t)$ defined via the equation: $\mathbf{\Phi}^l_i(t;\Lambda^h_i(0;x)) = \Lambda^h_i(t;\mathbf{\Phi}^h_i(t;x))$. By the invertability of $\Lambda^l_i$, this defines a flow, for which we can associate a push-forward map $\phi^l_t$ and inverse $\phi^l_{-t}$ as before. Note that this flow satisfies all properties <ref>–<ref> and <ref>–<ref> on the triplet $V(t) \subset H(t) \subset V^*(t)$ and therefore can be equipped with its own material derivative $\partial^l_t$ and bilinear form $\lambda^l(t; \cdot, \cdot)$ (Definition <ref>). Moreover, note that:

\begin{align*}
\partial^l_t \eta^{h,l} & = \phi^l_t \frac{d}{dt} \phi^{l}_{-t} (\eta_1^h(t; \Lambda_1^h(t;x)),\eta_2^h(t; \Lambda_2^h(t;x))) \\
&= \phi^l_t\frac{d}{dt} (\eta_1^h(t;\mathbf{\Phi}^h_1(t;( (\Lambda_1^{h})^{-1}(0;x)))),\eta_2^h(t;( (\Lambda_2^{h})^{-1}(0;x))))
= \bigg(\phi^h_t \frac{d}{dt} \phi^h_{-t} \eta^h\bigg)^l = (\partial^h_t \eta^h)^l,
\end{align*}

for $\eta^h \in C^1_{V^h}$. Let:

\begin{align*}
Z_k(t) := \{w \in V(t)|\; w_i \in H^{1+k}(\Omega_i(t))\}, \; \norm{w}^2_{Z_k(t)} := \sum_{i = 1}^2 \norm{w_i}^2_{H^{1+k}(\Omega_i(t))}.
\end{align*}

Assume that the pair $(Z_k(t), \phi_t)|_{t \in I}$ is compatible (which follows from the flow being of added regularity $\mathbf{\Phi}_i(t; \cdot) \in C^{1 +k}(\Omega_i(0); \mathbb{R}^d))$. Then the following variant of the approximation lemma holds:

We have the estimates:

\begin{align*}
\norm{w - I^l w}_{H(t)} +h \norm{w - I^lw}_{V(t)} & \leq c h^{k+1} \norm{w}_{Z_{k}(t)} && \text{for } w \in Z_k(t) \\
\norm{w - I^l w}_{H(t)} + h \norm{w - I^lw}_{V(t)} & \leq c h^2 \norm{w}_{Z_1(t)} && \text{for } w \in Z_1(t).
\end{align*}

See Section 7 in [20].

§ EVOLVING FINITE ELEMENT METHOD

§.§ Scheme

For $U^h, \, \zeta^h \in V^h(t)$, let:

\begin{align*}
m^h(t; U^h, \zeta^h) &:= \sum_{i = 1}^2 \int_{\Omega^h_i (t)} U^h_i \zeta^h_i,\\
a^h(t; U^h, \zeta^h) &:= \sum_{i = 1}^2 \int_{\Omega^h_i (t)} \mathcal{A}_i(t;\Lambda^{h}_l(t;x))\nabla U^h_i \cdot \nabla  \zeta^h_i + [\mathcal{B}(t;\Lambda^{h}_l(t;x)) - \mathbf{w}^h]\cdot \nabla U^h\, \zeta^h\\
& \qquad\qquad+ [\mathcal{C}(t; \Lambda^{h}_l(t;x)) - \nabla \cdot \mathbf{w}^h(t; x)] U^h \, \zeta^h,  \\
l^h(t; \zeta^h) & := (f^{-l}J^h, \zeta^h)_{H^h(t)} + (G^{-l}\mu^h, \zeta^h)_{L^2(\Gamma^h(t))},
\end{align*}

and set the initial condition:

\begin{align*}
U^h_0 = \sum_{j = 1}^{dim(\mathcal{S}^h)} (u_0, \chi_j^h)_{H(0)} \chi_j^h.
\end{align*}

Then, $U^h(t) \in \mathcal{S}^h(t)$ is a solution to the discrete variational problem if it solves:

\begin{align}\label{variational-eq}
m^h(t; \partial^h_t U^h, \zeta^h) + a^h(t; U^h, \zeta^h) +\lambda^h(t; U^h, \zeta^h)&= l^h(t; \zeta^h) \quad \forall \zeta^h \in L^2_{\mathcal{S}^h}, \forall t \in I, \\
U^h(0) &=U^h_0.
\end{align}

Here $J^h$, $\mu^h$ are the discrete Jacobians with respect to the lift maps $\Lambda^h_i(t;x)$, $\Lambda^h_1(t;x)|_{\Gamma^h(t)}$, and by regularity of $\Lambda^h$, are of class $C^k(K_i(t); \mathbb{R}^d)$, $C^k(E(t); \mathbb{R}^{d-1}), \; \forall K_i(t)\in \mathcal{J}^h_i(t), \; \forall E(t) \in \mathcal{G}^h_\Gamma(t)$ respectively.

Remark. It might not be practical to calculate $l^h(t;\cdot)$ for an arbitrary pair $(f,G) \in L^2_H \times L^2_{\mathcal{V}_\Gamma}$ as it would have to be calculated via numerical integration. We will assume for the rest of the paper that it is possible to calculate these integrals exactly, see [14] for numerical integration on curved domains.

This formulation can be rearranged to a more useful form via the transport theorem with respect the form $m^h(t; \cdot, \cdot)$:

\begin{equation}
\label{eq:weak-form}
\frac{d}{dt} m^h(t; U^h, \zeta^h) - m^h(t; U^h, \partial^h_t \zeta^h) +a^h(t; U^h, \zeta^h)= l^h(t; \zeta^h),
\end{equation}

for $\zeta^h \in C^1_{\mathcal{S}^h}$. Moreover, by construction of the $l^h(t;\cdot)$ term, for a function $\eta^h \in H^h(t)$:

\begin{align*}
l(t; \eta^{h,l}) -  l^h(t; \eta^{h} ) = 0.
\end{align*}

§.§ Well posedness of the finite element scheme

There exists a unique solution to Problem <ref> with continuous bound:

\begin{align*}
\sup_{t \in I}\norm{U^h}^2_{H^h(t)} +\int_0^T\norm{U^h}^2_{V^h(t)}\leq C(T)\bigg(\norm{U^{h}_0}^2_{H^h(0)} + \norm{f}_{L^2_H} +\norm{G}_{L^2_{\mathcal{V}_\Gamma^*}} \bigg).
\end{align*}

Substituting the Ansatz:

\begin{align*}
U^h(t;x) = \sum_{j= 1}^{N(k)} \alpha_j(t) \chi^h_j(t;x),
\end{align*}

where $\{\chi^h_j(t;x)\}^{N(k)}_{j = 1}$ are the basis functions of the evolving solution space $L^2_{\mathcal{S}^h}$. We refer to [20] Lemma 3.1 for proof of the transport property:

\begin{align*}
\partial^h_t \chi^h_j = 0 \quad \forall j \in {1,...., N(k)}.
\end{align*}

Then the problem can be restated as the finite dimensional problem:

\begin{align*}
\frac{d}{dt} \bigl(\mathbf{M}(t) \mathbf{\alpha}(t)\bigr) +\mathbf{A}(t)\mathbf{\alpha}(t) & = \mathbf{L}(t) \\
\mathbf{\alpha}(0) & = \mathbf{\alpha}_0
\end{align*}

\begin{align*}
\mathbf{\alpha} = (\alpha_1,....., \alpha_{N(k)}), \; [ \mathbf{M}]_{j,k} = m^h(t;\chi^h_j,\chi^h_k), \; [\mathbf{A}(t)]_{j,k} = a^h(t; \chi_j^h, \chi_k^k), \; \mathbf{L} = (l^h(t; \chi^h_1),...., l^h(t; \chi^h_{N(k)})).
\end{align*}

Note that $\mathbf{M}(t)$ is a Gram matrix (and hence invertible). Hence, by use of standard ODE theory (see Section 1.6 in [36]), there exists a solution $\alpha(t) \in \mathcal{W}^{1,1}(\mathbb{R}; \mathbb{R}^{N(k)})$. The uniform bound  and uniqueness follows from testing with  $U^h$ and using the transport theorem.

§ ERROR BOUND

The main result of this article is the following optimal order error bound.

If the solution to (<ref>) is of regularity $u \in W(Z_k, Z_k)$ with uniform bound:

\begin{align*}
\norm{u}_{W(Z_k, Z_k)} \leq C_u,
\end{align*}

then there exists a constant $\mathcal{C}$ depending on $C_u$ such that the following holds:

\begin{align*}
\sup_{t \in I} \norm{u - U^{h,l}}^2_{H(t)} + h^2\int_0^T \norm{u - U^{h,l} }^2_{V(t)} \leq \norm{u_0 - u_0^{h,l}}_{H(0)}^2 +h^{2k+2} \mathcal{C}(C_u).
\end{align*}

In the next section we set out preliminary approximation results and then prove the error bound in the subsequent section.

§.§ Preliminaries

We introduce only the necessary tools from Section 3.3 in [20] in order to obtain suitable orders of convergence. Since we have two possible flows on $V(t)$, $\phi$ and $\phi^l$, we need to guarantee some compatibility condition that are reflected on the problem, to that end, note that we can define the derivative of the bilinear forms with respect to $\mathbf{\Phi}^l$ instead. Let:

\begin{align*}
b^l(t; u,v) := \frac{d}{dt}[a(t;u,v)] - a(t; \partial_t^l u,v) - a(t;u, \partial_t^l v),
\end{align*}

and define $\lambda^l(t;\cdot , \cdot)$ to be the bilinear form of Definition <ref> with respect to the flow $\mathbf{\Phi}^l$, which can be calculated to be:

\begin{align*}
\lambda^l(t;v, w) = \sum_{i = 1}^2 \int_{\Omega_i(t)} \nabla \cdot \mathbf{w}^l v\, w,
\end{align*}

where $\mathbf{w}^l$ is defined as:

\begin{align*}
\mathbf{w}^l(t;x) = \frac{d}{dt} \mathbf{\Phi}^l(t;y)|_{y = \mathbf{\Phi}^l(-t;x)}.
\end{align*}

Then the following holds for the bilinear forms introduced in (<ref>) and (<ref>):

There exists a constant $c >0$ such that for almost all $t \in I$ and for all $w^h,v^h \in V^h(t)$, $w^{h,l}, v^{h,l} \in V(t)$ the following error bounds hold:

\begin{align}
\label{eq:m-error} \tag{P1}
|m( t; w^{h,l}, v^{h,l} ) - m^h( t; w^h, v^h ) |& \leq c h^{k+1} \norm{w^{h,l}}_{V(t)} \norm{v^{h,l}}_{V(t)}, \\
\label{eq:c-error} \tag{P2}
|\lambda( t; w^{h,l}, v^{h,l} ) - \lambda^h( t; w^h, v^h )| &\leq c h^{k+1} \norm{w^{h,l}}_{V(t)}\norm{v^{h,l}}_{V(t)}, \\
\label{eq:ct-error} \tag{P3}
|\lambda^l( t; w^{h,l}, v^{h,l} ) - \lambda( t; w^{h,l}, v^{h,l} )| &\leq c h^{k} \norm{w^{h,l}}_{V(t)}\norm{v^{h,l}}_{V(t)}, \\
\label{eq:a-error} \tag{P4}
|a( t; w^{h,l}, v^{h,l} ) - a^{h}( t; w^h, v^h )|& \leq c h^{k} \norm{w^{h,l}}_{V(t)} \norm{v^{h,l}}_{V(t)}, \\
\label{eq:b-error} \tag{P5}
|b^l( t; w^{h,l}, v^{h,l}) - b^h(t;w^h,v^h)|&\leq ch^k \norm{w^{h,l}}_{V(t)} \norm{v^{h,l}}_{V(t)}, \\
\label{eq:bt-error} \tag{P6}
|b^l( t; w^{h,l}, v^{h,l} ) - b( t; w^{h,l}, v^{h,l} )|& \leq c h^{k} \norm{w^{h,l}}_{V(t)} \norm{v^{h,l}}_{V(t)}.
\end{align}

For $\eta, \zeta \in Z_1(t)$ with inverse lifts $\eta^{-l}, \zeta^{-l}$:

\begin{align}
\label{eq:a-error2} \tag{\ref*{eq:a-error}'}
|a( t; \eta, \zeta ) - a^{h}( t; \eta^{-l}, \zeta^{-l} ) |& \leq ch^{k+1}\norm{\eta}_{Z_1(t)}\norm{\zeta}_{Z_1(t)}, \\
\label{eq:b-error2} \tag{\ref*{eq:b-error}'}
|b^l( t; \eta, \zeta) - b^{h}( t; \eta^{-l}, \zeta^{-l})|& \leq c h^{k+1} \norm{\eta}_{Z_1(t)} \norm{\zeta}_{Z_1(t)}.
\end{align}

For $\eta \in C^1_{Z_1}$ and $\zeta \in Z_1(t)$, with inverse lifts $\eta^{-l}$ and $\zeta^{-l}$:

\begin{align}
\label{eq:amd-error2} \tag{P7}
|a( t; \partial_t^l \eta, \zeta ) - a^{h}( t; \partial_t^h  \eta^{-l}, \zeta^{-l} )|
&\leq c h^{k+1} ( \norm{\eta}_{Z_1(t)} +\norm{\partial_t^\bullet\eta}_{Z_1(t)} )\norm{\zeta}_{Z_1(t)}.
\end{align}

The material derivatives satisfy

∂_t^l ζ- ∂_t^∙ζ_H(t)
≤c h^k+1 ζ_V(t)

∂^l_t ζ- ∂^∙_t ζ_V(t)
≤c h^k ζ_Z_1(t)

Those statement follow by applying the estimates from [20] Section 8.

Before moving on, we slightly modify the problem, for  a test function $v \in W(V,H)$, we can rewrite the weak formulation of our problem as:

\begin{align*}
\frac{d}{dt} m(t;u,v) - m(t;u,\partial^\bullet_t v) +a(t;u,v) = l(t;v).
\end{align*}

Doing the standard parabolic rescaling $\check{u} = e^{-\kappa t}u$, the problem becomes:

\begin{align}\label{eq:error}
\frac{d}{dt}m(t;\check{u}, v) - m(t;u,\partial^\bullet_t v)+ \underbrace{\kappa m(t; \check{u}, v)+a(t;\check{u}, v)}_{=:a_\kappa(t;\check{u},v)} = \underbrace{e^{-\kappa t}l(v)}_{ =:\tilde{l}(v)}.
\end{align}

Testing with $v = \check{u}$, we arrive at (using Young's inequality):

\begin{align*}
a_\kappa(t; \check{u}, \check{u})
& = \sum_{i=1}^2 \int_{\Omega_i(t)} \mathcal{A} \nabla \check{u} \cdot \nabla \check{u}
+ (\mathcal{B} - \mathbf{w}) \nabla \check{u} \cdot \check{u}
+ (\mathcal{C} - \nabla \cdot \mathbf{w}) \check{u}^2 \\
& \ge \gamma \norm{\check{u}}^2_V
- \norm{\mathcal{B} - \mathbf{w}}_{L^\infty} \norm{\check{u}}_V \norm{\check{u}}_{H}
+ (\kappa - \norm{\mathcal{C}} - \norm{\nabla \cdot \mathbf{w}}_{L^\infty}) \norm{u}_{H}^2 \\
& \ge (\gamma - \frac{\epsilon}{2} \norm{\mathcal{B} - \mathbf{w}}_{L^\infty}) \norm{ \check{u}}_V^2
+ (\kappa - \norm{\mathcal{C}} - \norm{\nabla \cdot \mathbf{w}} - \frac{\norm{\mathcal{B} - \mathbf{w}}}{2\epsilon}).
\end{align*}

\begin{multline*}
a_\kappa(t;\check{u}, \check{u}) \geq \bigg(\gamma  - \frac{\epsilon \norm{\mathcal{B} - \mathbf{w}}_{L^\infty(\Omega)}}{2}\bigg)\norm{\nabla \check{u}}^2_{H(t)} \\
+\bigg(\kappa - \norm{\nabla \cdot \mathbf{w}}_{L^\infty(\Omega)} - \norm{\mathcal{C}}_{L^\infty(\Omega)} - \frac{\norm{\mathcal{B} - \mathbf{w}}_{L^\infty(\Omega)}}{2 \epsilon} \bigg) \norm{\check{u}}^2_{H(t)}.
\end{multline*}

Taking $\epsilon < 2 \gamma\norm{\mathcal{B - \mathbf{w}}}_{L^\infty(\Omega)}^{-1}$ and $\kappa >\norm{\nabla \cdot \mathbf{w}}_{L^\infty(\Omega)} +\norm{\mathcal{C}}_{L^\infty(\Omega)} + (4\gamma)^{-1}\norm{\mathcal{B - \mathbf{w}}}^2_{L^\infty(\Omega)}$ guarantees the coercivity of the bilinear form $a_\kappa$.
Performing the same transformation to the discrete analogue: define $\check{U}_h = e^{-\kappa t} U_h$ which satisfies

\begin{align}\label{eq:errordis}
\frac{d}{dt}m^h(t;\check{U}^h, v^h) - m^h(t;\check{U},\partial^h_t v^h)+ \underbrace{\kappa m^h(t; \check{U}^h, v^h)+a^h(t;\check{U}^h, v^h)}_{=:a^h_\kappa(t;\check{U}^h,v^h)} = \underbrace{e^{-\kappa t}l^h(v^h)}_{ =:\check{l}^h(v^h)},
\end{align}

which is also coercive provided $\kappa$ is large enough by the same argument. The Ritz projection $\Pi^h :V(t) \to V^h(t)$ with respect to those bilinear forms is defined as the solution to:

\begin{align*}
a^h_\kappa(t; \Pi^h(u), v^h) = a_\kappa(t; u, v^{h,l}) \quad \forall v^h \in V^h(t),
\end{align*}

by the coercivity and boundedness, this gives us a bounded and linear operator $\Pi^h:V(t) \to V^h(t)$.

(Lemma 3.8 - 3.10 in [20]).
Let $\pi^h z := (\Pi^h z)^l$, then, for $w \in C^1_{Z_k}$, if (<ref>)–(<ref>), Lemma <ref> and (<ref>) hold, there exists a $c >0$ such that:

\begin{align*}
\norm{\partial^h_t \Pi^h w}_{V^h(t)} &\leq \mathcal{C}(\norm{w} + \norm{\partial^\bullet_t w}_{V(t)}),\\
\norm{w - \pi^h w}_{H(t)} +h \norm{w - \pi^hw}_{V(t)} &\leq c h^{k+1} \norm{w}_{Z_k(t)}, \\
\norm{\partial^l_t ( w - \pi^h w)}_{H(t)} + h \norm{\partial^l_t(w - \pi^hw)}_{V(t)} &\leq c h^{k+1}(\norm{w}_{Z_k(t)} +\norm{\partial_t^\bullet w}_{Z_k(t)}).
\end{align*}

The same duality argument as in [20] and [16] is used. Set $\mathcal{T}:\mathcal{H}_\Gamma(t) \to \mathcal{V}_\Gamma(t)$ as:

\begin{align*}
(\mathcal{T}w, v)_{\mathcal{V}_\Gamma(t)} := \langle w, v \rangle_{\mathcal{V}_\Gamma(t)} = (w,v)_{\mathcal{H}_\Gamma(t)}
\text{ for all } w, v \in \mathcal{V}_\Gamma(t),
\end{align*}

i.e $\mathcal{T}$ acts a Riesz map mapping to the element in $\mathcal{V}_\Gamma(t)$ that corresponds to the functionals in $\mathcal{H}_\Gamma(t) \subset \mathcal{V}_\Gamma^*(t)$. Notice that:

\begin{align*}
\norm{\mathcal{T}w}^2_{\mathcal{V}_\Gamma(t)} = \norm{w}^2_{\mathcal{V}^*_\Gamma(t)} = \int_{\Gamma(t)} w \mathcal{T}w,
\end{align*}

for any $w \in \mathcal{H}_\Gamma(t)$. Define the solution operator $\mathcal{R}: \mathcal{V}_\Gamma(t) \to Z_1(t)$ to be the solution to:

\begin{align*}
a_\kappa(t; \mathcal{R}(v), w) = (v,w)_{\mathcal{H}_\Gamma(t)} \quad \forall w \in V(t).
\end{align*}

Since the trace of a function $w \in V(t)$ is a function in $\mathcal{H}_\Gamma(t)$, the problem is well posed and the solution satisfies the regularity result (as in [29] Theorem 1):

\begin{align*}
\norm{\mathcal{R}(v)}_{Z_1(t)} \leq C\norm{v}_{\mathcal{V}_\Gamma(t)},
\end{align*}

where the regularity constant is bounded on $[0, T]$.
Next, we substitute $v =\mathcal{T}(\xi)$:

\begin{align}
\label{eq:trace-reg}
\norm{\mathcal{R} \circ \mathcal{T}(\xi)}_{Z_1(t)} \leq C\norm{\mathcal{T}(\xi)}_{\mathcal{V}_\Gamma(t)} = C\norm{\xi}_{\mathcal{V}^*_\Gamma(t)} = C(\xi, \mathcal{T}\xi)_{\mathcal{H}_\Gamma(t)}^{1/2}.
\end{align}

This allows us to prove the following:

For $w \in Z_k(t)$ and $\eta  := w - \pi^h w$ it holds that

\begin{align}\label{eq:bil}
|b(t; \eta,v)| \leq c (\norm{\eta}_{H(t)} +h\norm{\eta}_{V(t)} +h^{k+1}\norm{w}_{Z_k(t)}) \norm{v}_{Z_1(t)}.
\end{align}

By construction:

\begin{align*}
\norm{v}^2_{\mathcal{V}_\Gamma^*(t)} = \int_{\Gamma(t)} v \cdot \mathcal{T}(v) = a_\kappa(t;v, \mathcal{R}\circ \mathcal{T}(v)) = a_\kappa(t; v, \mathcal{R}\circ \mathcal{T}(v) - I^l [\mathcal{R}\circ \mathcal{T}(v)]) + a_\kappa(t; v, I^l[ \mathcal{R}\circ \mathcal{T}(v)]).
\end{align*}

Then, for the first part, we have:

\begin{align*}
|a_\kappa(t; v, \mathcal{R}\circ \mathcal{T}(v) - I^l [\mathcal{R}\circ \mathcal{T}(v))]| &\leq \norm{v}_{V(t)} \norm{\mathcal{R}\circ \mathcal{T}(v) - I^l [\mathcal{R}\circ \mathcal{T}(v)]}_{V(t)}, \\
& \leq ch\norm{v}_{V(t)} \norm{\mathcal{R}\circ \mathcal{T}(v)}_{Z_1(t)} \leq ch \norm{v}_{V(t)}\norm{v}_{\mathcal{V}_\Gamma^*(t)}.
\end{align*}

Hence we have:

\begin{align}\label{b-ineq}
\norm{v}^2_{\mathcal{V}_\Gamma^*(t)} \leq ch \norm{v}_{V(t)}\norm{v}_{\mathcal{V}_\Gamma^*(t)} + a_\kappa(t; v, I^l [\mathcal{R}\circ \mathcal{T}(v)]).
\end{align}

Going back to the bilinear form $b(t;\cdot, \cdot)$, using it's explicit form <ref>, for $\eta:= u - \pi^h u$:

\begin{align*}
|b(t; \eta,v)| & \leq C\norm{\eta}_{H(t)} \norm{v}_{H(t)} + \bigg|\sum_{i = 1}^2  \int_{\Omega_i(t)} \mathcal{D}^B (\mathbf{w}, \mathcal{B}, \eta, v) + \mathcal{D}^A(\mathbf{w}, \mathcal{A}_i, \eta_i, v_i) \bigg|.
\end{align*}

For $\mathcal{D}^A$, integrating by parts yields:

\begin{align*}
\bigg| \sum_{i = 1}^2 \int_{\Omega_i(t)} &\mathcal{D}^A(\mathbf{w}, \mathcal{A}_i, \eta_i, v_i) \bigg| \leq \bigg| \int_{\Gamma(t)} \bigg \llbracket (\partial_t^\bullet \mathcal{A}_i+ \nabla\cdot \mathbf{w} \mathcal{A}_i) \frac{\partial v}{\partial \nu_\Gamma} \bigg\rrbracket \eta - \sum_{i = 1}^2 \int_{\Omega_i(t)} \eta_i \nabla\cdot ([\partial^\bullet_t \mathcal{A}_i + \nabla \cdot\mathbf{w} \mathcal{A}_i]\nabla v_i) \\
& +\sum_{i = 1}^2 \int_{\Omega_i(t)} \eta_i \sum_{l,k = 1}^d \nabla_l(\mathcal{A}_i[\nabla_l \mathbf{w}_k +\nabla_k \mathbf{w}_l]\nabla_k v_i)  - \int_{\Gamma(t)}\eta \sum_{l,k = 1}^d \llbracket \nu_l \mathcal{A}_i[\nabla_l \mathbf{w}_k +\nabla_k \mathbf{w}_l]\nabla_k v_i\rrbracket \bigg|, \\
&\leq C(|\mathcal{A}_i|_{C^1(\overline{\mathcal{Q}}_i; \mathbb{R})}, |\nabla \mathbf{w}|_{C^1(\overline{\Omega}; \mathbb{R}^d)}) \norm{v}_{Z_1(t)} (\norm{\eta}_{H(t)} + \norm{\eta}_{\mathcal{V}^*_{\Gamma}(t)}).
\end{align*}

In the last line we have used both the generalised trace inequality and the Banach triple identification for the boundary terms. Similarly for $\mathcal{D}^B$:

\begin{align*}
\bigg| \sum_{i = 1}^2 \int_{\Omega_i(t)} &\mathcal{D}^\mathcal{B}(\mathbf{w}, \mathcal{B}, \eta_i, v_i) \bigg| \leq C(|B|_{C^1(\overline{\Omega}; \mathbb{R}^d)}, |\nabla\cdot \mathbf{w}|_{C^1(\overline{\Omega}; \mathbb{R}^d)} )\norm{v}_{Z_1(t)} (\norm{\eta}_{H(t)} + \norm{\eta}_{\mathcal{V}^*_{\Gamma}(t)}).
\end{align*}

Using (<ref>):

\begin{align} \label{b-ineq2}
|b(t; \eta,v)| & \leq C\bigg( \norm{\eta}_{H(t)}  +  h\norm{\eta}_{V(t)} + \frac{1}{\norm{\eta}_{\mathcal{V}_\Gamma^*(t)}}a_\kappa(t;\eta, I^l\mathcal{R} \circ \mathcal{T}(\eta))\bigg)\norm{v}_{Z_1(t)}.
\end{align}

Finally, this yields:

\begin{align*}
|a_\kappa(t;\eta, I^l[\mathcal{R} \circ \mathcal{T}(\eta)])|
&= |a_\kappa(t; \pi^h w, I^l [\mathcal{R}\circ \mathcal{T}(\eta)]) - a_\kappa^h(\Pi^hw, I^h [\mathcal{R} \circ \mathcal{T}(\eta)]) |\\
& \leq |a_\kappa(\pi^h w - w, I^l [\mathcal{R} \circ \mathcal{T}(\eta)]) - a_\kappa^h (\Pi^h w - w^{-l}, I^h [\mathcal{R}\circ \mathcal{T}(\eta)]) | \\
& +|a_\kappa(w, I^l [\mathcal{R} \circ \mathcal{T}(\eta)] - \mathcal{R}\circ \mathcal{T}(\eta)) - a^h_\kappa(w^{-l}, I^h [\mathcal{R} \circ \mathcal{T}(\eta)^{-l} ]- \mathcal{R} \circ\mathcal{T}(\eta)^{-l})| \\
& + |a_\kappa(w, \mathcal{R} \circ \mathcal{T}(\eta)) - a_\kappa^h(w^{-l}, \mathcal{R} \circ \mathcal{T}(\eta)^{-l}) | \\
%& \leq c^{h^k} \norm{\pi^h w - w}_{V(t)} \norm{I^l[\mathcal{R} \circ \mathcal{T}(\eta)]}_{V(t)} \\
%& + c h^{k+1} \norm{w}_{V(t)} \norm{\mathcal{R} \circ \mathcal{T}(\eta)}_{Z_1(t)} \\
%& + c h^{2k+1} \norm{w}_{Z_1(t)} \norm{\mathcal{R} \circ \mathcal{T}(\eta)}_{Z_1(t)} \\
& \leq c h^{k + 1} \norm{w}_{Z_1(t)} \norm{\mathcal{R} \circ \mathcal{T}(\eta)}_{V(t)}
+ c h^{k+1} \norm{w}_{V(t)} \norm{\eta}_{\mathcal{V}_\Gamma^*(t)}
+ c h^{k+1} \norm{w}_{Z_1(t)} \norm{\eta}_{\mathcal{V}_\Gamma^*(t)} \\
&\leq ch^{k+1} \norm{w}_{Z_1(t)} \norm{\eta}_{\mathcal{V}_\Gamma^*(t)},
\end{align*}

using Lemmas <ref> and <ref>, the regularity estimate above (<ref>) and (<ref>), (<ref>) and (<ref>). Substituting the final inequality here in (<ref>) yields (<ref>).

§.§ Proof of error bound

We have an additional right-hand side functional term that is not present in the  original proof by Elliott $\&$ Ranner. This requires  a modification of the proof of Theorem 3.11 in [20].

Set $\theta :=  \check{U} - \Pi^h \check{u}$, then using (<ref>) and using the fact that $l^h(t;\cdot)$ equals $l(t; (\cdot)^l)$ for functions in $H^h(t)$, we arrive at:

\begin{align*}
&\frac{d}{dt}m^h(t; \Pi^h \check{u}, v^h) +a^h_\kappa(t; \Pi^h\check{u}, v^h) - m^h(t; \Pi^h\check{u}, \partial^h_t v^h) - \tilde{l}^h(v^h), \\
&= \frac{d}{dt}m^h(t; \Pi^h \check{u}, v^h) +a_\kappa(t; \check{u}, v^{h,l}) - m^h(t; \Pi^h\check{u}, \partial^h_t v^h) - \tilde{l}^h(v^h) ,\\
&= \frac{d}{dt}\bigg[m^h(t; \Pi^h \check{u}, v^h) - m(t; \check{u}, v^{h,l}) \bigg] - \bigg[m^h(t; \Pi^h \check{u}, \partial^h_t v^h) - m(t; \check{u}, \partial^\bullet_t v^{h,l})\bigg],\\
&= m^h(t; \partial_t^h \Pi^h \check{u}, v^h) - m(t; \partial^l_t \check{u}, v^{h,l}) +\lambda^h(t; \Pi^h\check{u}, v^{h,l}) - \lambda^l(t;\check{u}, v^{h,l}) +m(t; \check{u}, \partial^\bullet_t v^{h,l} - \partial^l_t v^{h,l}).
\end{align*}

Now subtracting this equation from (<ref>), and rearranging yields:

\begin{align}\label{eq:erroru}
&\frac{d}{dt}m^h(t;\theta, v^h) +a^h_\kappa(t; \theta, v^h) - m^h(t; \theta, \partial^h_t v^h) \\ \nonumber
&= - \bigg[m^h(t; \partial_t^h \Pi^h \check{u}, v^h)  - m(t; \partial^l_t \pi^h\check{u}, v^{h,l}) + m(t; \partial^l_t [\pi^h\check{u}-\check{u}], v^{h,l}) +\lambda^h(t; \Pi^h\check{u}, v^{h,l}) - \lambda^l(t;\pi^h\check{u}, v^{h,l})\\ \nonumber
&\qquad +\lambda^l(t;[\pi^h\check{u} - \check{u}], v^{h,l}) +m(t; \check{u}, \partial^\bullet_t v^{h,l} - \partial^l_t v^{h,l})\bigg]=:-\Xi^h(\check{u}, v^h).
\end{align}

Using the identity $\partial^l_t (v^{l,h}) = (\partial^h_t v^h)^l$ and looking at $\Xi(\cdot, \cdot)$ term by term:

\begin{align*}
& |m^h(t; \partial_t^h \Pi^h \check{u}, v^h)  - m(t; \partial^l_t \pi^h\check{u}, v^{h,l})| \\
&= |m^h(t; \partial_t^h \Pi^h \check{u}, v^h) -m(t;(\partial_t^h \Pi^h \check{u})^l, v^{h,l} ) +  m(t; (\partial_t^h \Pi^h\check{u})^l  - \partial^l_t \pi^h\check{u}, v^{h,l})|, \\
&\leq ch^{k+1} (\norm{\check{u}}_{Z_k(t)} + \norm{\partial^\bullet_t \check{u}}_{Z_k(t)})\norm{v^h}_{V^h(t)},
\end{align*}

by (<ref>). Similar rearrangement and the use of Lemma <ref> with (<ref>), (<ref>) ,(<ref>) yields:

\begin{align*}
|\Xi^h(\check{u}, v^h)| \leq ch^{k+1} (\norm{\check{u}}_{Z_k(t)} + \norm{\partial^\bullet_t \check{u}}_{Z_k(t)} ) \norm{v^h}_{V^h(t)}.
\end{align*}

Using (<ref>) and substituting $v^h = \theta$, we obtain:

\begin{align*}
\frac{d}{dt}m^h(t;\theta, \theta) +a^h_\kappa(t; \theta, \theta) - m^h(t; \theta, \partial^h_t \theta) = \Xi^h(\theta,\theta).
\end{align*}

Using the transport formula and the bound on $\Xi$:

\begin{align*}
\frac{1}{2} \frac{d}{dt}m^h(t; \theta, \theta) +a^h_\kappa(t; \theta, \theta) \leq -\frac{1}{2} \lambda^h(t;\theta,\theta) +ch^{k+1} (\norm{\check{u}}_{Z_k(t)} + \norm{\partial^\bullet_t \check{u}}_{Z_k(t)} ) \norm{\theta}_{V^h(t)}.
\end{align*}

Integrating over time and using Young's and Gröwal's inequality:

\begin{align*}
\sup_{t \in I} \norm{\theta}^2_{H^h(t)} +\int_0^T \norm{\theta}^2_{V^h(t)} \leq \norm{\theta}^2_{H^h(0)} +C\int_0^T \norm{\theta}^2_{H^h(t)} +ch^{2k+2} \int_0^T (\norm{\check{u}}^2_{Z_k(t)} + \norm{\partial^\bullet_t \check{u}}^2_{Z_k(t)}).
\end{align*}

Finally, using the decomposition:

\begin{align*}
\check{u} - {\check{U}}^{h,l} = \check{u} - \pi^h\check{u} +\pi^h\check{u} - {\check{U}}^{h,l}.
\end{align*}

Using the previous bound, the fact that the lift is a diffeomorphism and the bound on the Ritz map, we finally obtain:

\begin{align*}
\sup_{t \in I} \norm{\check{u} - {\check{U}}^{h,l} }^2_{H(t)} +h^2\int_0^T \norm{\check{u} - {\check{U}}^{h,l} }^2_{V(s)}\,ds &= \sup_{t \in I} \norm{\check{u} - \pi^h\check{u} +\theta^l}^2_{H(t)} +h^2\int_0^T \norm{\check{u}- \pi^h\check{u} +\theta^{l} }^2_{V(s)}\,ds, \\
&\leq \norm{u_0 - u_0^{h,l}}_{H(0)} +h^{2k+2}c(C_u).
\end{align*}

Undoing the scaling $u = e^{\kappa t}\check{u}$ gives us the desired error bound.

§ NUMERICAL RESULTS

All numerical results are computed using the  package[7, 8, 15, 35].
Simulation code is available in [34].
Results are computed on a sequence of meshes generated using [22] rather than successive refinement.

The main challenges in implementing the numerical scheme are:

* computing the initial geometry:
We start with a piecewise linear geometry given by .
The initial domain is computed through an explicit parametrisation applying directly the method from Section <ref> efficiently using custom written  code.
The domain evolution is carried out simply by moving the initial Lagrange nodes according to the smooth, given velocity field.

* labelling and tracking different parts of the domain:
Along side the geometry and topology of the mesh we must track labels which say which elements are in domain $\Omega^h_1(t)$ or $\Omega^h_2(t)$ and which facets are on $\Gamma^h(t)$. Once this is fixed for the initial domains $\tilde{\Omega}^h_1$ and $\tilde{\Omega}^h_2$, this information is passed between different times.
provides physical tags to each element and facets which can be used to identify the different domains.

Efficient and accurate quadrature rules are used to perform element-wise integrals.
Note that system matrices must be reassembled at each time step due to the evolution of the domain.

§.§ Time discretisation of advection-diffusion problem

We start from the spatial discretisation from Section <ref>. We will apply a backward difference formula (BDF) time discretisation of order $q$, see [28] for more including analysis of a similar surface only problem.
We take a partition of the time interval $0 = t_0 < t_1 < \cdots < t_M = T$.
For simplicity we assume that each time interval is of the same length: $\tau := t_j - t_{j-1}$ for $j = 1, 2, \ldots, M$.

We use temporal interpolations of each domain at each time step to construct a sequence of triangulations  $\mathcal{J}^h(t_j)$ each equipped with finite element spaces $\mathcal{S}^h(t_j)$ for $j=0, 1, \ldots, M$.
We define the discrete velocity $W^j \in S_h(t_j)^d$ by

\begin{equation}
W^j = \frac{1}{\tau} \sum_{l=0}^q \delta_l X^{j-l},
\end{equation}

where $X^{j}$ are the positions of the Lagrange nodes of the triangulation at time $t_j$ and $\{ \delta_l \}_{l=0}^q$ are the backward difference formula weights determined from the relation:

\begin{equation}
\label{eq:bdf-weights}
\delta(\zeta) = \sum_{l=0}^q \delta_j \zeta^q
= \sum_{\lambda=1}^q \frac{1}{\lambda} (1 - \zeta)^\lambda.
\end{equation}

The fully discrete problem is the time discritisation of (<ref>):
Given starting values $U^0 \in S^h(t_0)$, $\ldots$, $U^{q-1} \in S^h(t_0)$, and data $\mathcal{A}, \mathcal{B}, \mathcal{C}$ and $l_h$, for $j = q, \ldots, M$, we wish to find $U^j \in S^h(t_j)$ as the solution of

\begin{multline}
\label{eq:full-discrete}
\frac{1}{\tau} \sum_{l=0}^q \delta_{l} m^h(t^{j-l}; U^{j-l}, \chi_i^{j-l})
+ a^h(t^j; U^j, \chi^j))
=  l_h(t^j; \chi^i_j)
\qquad \mbox{ for all basis function } \chi_i^j \in V_h^j,
\end{multline}

where again $\delta_l$ are the BDF weights (<ref>). Note that the first term on the left hand side is computed by summing over $q$ different meshes to approximate the time derivative.

§.§ Numerical examples of advection-diffusion problem

For $d = 2, 3$, let $\Omega = [-1,1]^d$, for $t \in [0,T]$, we define the evolution of the domain through the flow map $\Phi_t$ given by:

\begin{align*}
\Phi_t(x) = x + \frac{|x|^{1/3}\prod_{i=1}^d(1-x_i^2)}{0.5 \prod_{i=1}^d (1-4 x_i^2 / |x|)}
\begin{cases}
((\alpha(t)-1) x_1, (\beta(t)-1) x_2) & \mbox{ if } d = 2, \\
((\alpha(t)-1) x_1, (\beta(t)-1) x_2, 0) & \mbox{ if } d = 3,
\end{cases}
\end{align*}

for $\alpha(t) = 1 + 0.25 \sin(t)$ and $\beta = 1 + 0.25 \cos(t)$.
This is a special motion which ensures that notes initially on $\partial \Omega$ do not move and the surface $\Gamma(t)$ is described by the level set function $\phi(\cdot, t)$ given by

\begin{align*}
\phi(\cdot, t) = \begin{cases}
\frac{x_1^2}{\alpha(t)^2} + \frac{x_2^2}{\beta(t)^2} - \frac{1}{2}
& \mbox{ if } d = 2 \\
\frac{x_1^2}{\alpha(t)^2} + \frac{x_2^2}{\beta(t)^2} + x_3^2 - \frac{1}{2}
& \mbox{ if } d = 3.
\end{cases}
\end{align*}

We define $\Omega_1(t)$ as the interior of $\Gamma(t)$ and $\Omega_2(t) = \Omega \setminus \Omega_1(t)$.

We set $\mathcal{A}|_{\Omega_1(t)} = 10$, $\mathcal{A}|_{\Omega_2(t)} = 1$, $\mathcal{B} = 0$ and $\mathcal{C} = 0$.
We set data, $l_h$, such that the exact solution $u$ is given by

\[
u(x, t) =
\sin(t) |\Phi(x)| \prod_{i=1}^d \sin(2 \pi x_i).
\]

In order to simplify the implementation the right hand side data ($l_h$) is computed by taking interpolations of smooth data.
We compute using isoparametric elements of order 1, 2, 3 on a sequence of given meshes.
For order $k$ discretisation in space we use BDF order $k+1$ in time.
The initial solution $U^0 = 0$ matches the exact solution at $t = 0$. The other starting values are computed using lower order BDF methods.
For elements of order $k$ we expect convergence of order $k+1$ for the error at the final time, $u(T) - U^M$, in the $L^2(\Omega)$ norm and order $k$ is the $H^1(\Omega)$ semi-norm.
The results are shown in Figure <ref> for the cases $d=2,3$ respectively.
The precise numerical values are shown in Table <ref>.
We see that the numerical results support for analytical convergence results.

$L^2$ error for advection-diffusion problem for $d=2$ (top) and $d=3$ (bottom).

$h$                       $\tau$                  $L^2$ error                  $H^1$ error     eoc($L^2$ error)     eoc($H^1$ error)

$1.00000$                    $1.00000$     $3.61161 \times 10^{-1}$                    $3.33340$                  —                  —

$6.66667 \times 10^{-1}$     $5.00000 \times 10^{-1}$                    $1.61100$                    $5.31814$            -3.687830            -1.152088

$4.37747 \times 10^{-1}$     $2.50000 \times 10^{-1}$     $2.03904 \times 10^{-1}$                    $2.34493$             4.913737             1.946674

$2.40008 \times 10^{-1}$     $1.25000 \times 10^{-1}$     $9.43824 \times 10^{-2}$                    $1.45816$             1.281758             0.790531

$1.34513 \times 10^{-1}$     $6.25000 \times 10^{-2}$     $2.92209 \times 10^{-2}$     $7.82215 \times 10^{-1}$             2.024952             1.075622

Order 1, $d=2$

$h$                       $\tau$                  $L^2$ error                  $H^1$ error     eoc($L^2$ error)     eoc($H^1$ error)

$1.00000$                    $1.00000$     $3.04065 \times 10^{-1}$                    $3.34129$                  —                  —

$6.66667 \times 10^{-1}$     $5.00000 \times 10^{-1}$     $2.93554 \times 10^{-1}$                    $2.04598$             0.086769             1.209675

$4.37747 \times 10^{-1}$     $2.50000 \times 10^{-1}$     $5.19037 \times 10^{-2}$     $8.64156 \times 10^{-1}$             4.119035             2.048919

$2.40008 \times 10^{-1}$     $1.25000 \times 10^{-1}$     $6.67721 \times 10^{-3}$     $2.89035 \times 10^{-1}$             3.412317             1.822405

$1.34513 \times 10^{-1}$     $6.25000 \times 10^{-2}$     $9.32995 \times 10^{-4}$     $8.42882 \times 10^{-2}$             3.398988             2.128289

Order 2, $d=2$

$h$                       $\tau$                  $L^2$ error                  $H^1$ error     eoc($L^2$ error)     eoc($H^1$ error)

$1.00000$                    $1.00000$     $1.96574 \times 10^{-1}$                    $2.56225$                  —                  —

$6.66667 \times 10^{-1}$     $5.00000 \times 10^{-1}$     $1.07806 \times 10^{-1}$     $9.17747 \times 10^{-1}$             1.481515             2.532199

$4.37747 \times 10^{-1}$     $2.50000 \times 10^{-1}$     $1.50610 \times 10^{-2}$     $1.99583 \times 10^{-1}$             4.679020             3.626986

$2.40008 \times 10^{-1}$     $1.25000 \times 10^{-1}$     $6.00060 \times 10^{-4}$     $4.39335 \times 10^{-2}$             5.362743             2.518530

$1.34513 \times 10^{-1}$     $6.25000 \times 10^{-2}$     $5.81492 \times 10^{-5}$     $2.73137 \times 10^{-2}$             4.031031             0.820862

Order 3, $d=2$

$h$                       $\tau$                  $L^2$ error     $H^1$ error     eoc($L^2$ error)     eoc($H^1$ error)

$1.25000$                    $1.00000$     $6.90652 \times 10^{-1}$       $8.44874$                  —                  —

$8.66599 \times 10^{-1}$     $5.00000 \times 10^{-1}$     $6.57805 \times 10^{-1}$       $7.69627$             0.133018             0.254644

$6.31590 \times 10^{-1}$     $2.50000 \times 10^{-1}$     $5.80704 \times 10^{-1}$       $6.69577$             0.394097             0.440229

$3.33531 \times 10^{-1}$     $1.25000 \times 10^{-1}$     $3.16133 \times 10^{-1}$       $4.46142$             0.952350             0.635873

$1.75870 \times 10^{-1}$     $6.25000 \times 10^{-2}$     $1.01279 \times 10^{-1}$       $2.34846$             1.778601             1.002685

Order 1, $d=3$

$h$                       $\tau$                  $L^2$ error                  $H^1$ error     eoc($L^2$ error)     eoc($H^1$ error)

$8.66599 \times 10^{-1}$     $5.00000 \times 10^{-1}$     $3.87129 \times 10^{-1}$                    $5.41237$                  —                  —

$6.31590 \times 10^{-1}$     $2.50000 \times 10^{-1}$     $1.72322 \times 10^{-1}$                    $2.90077$             2.558644             1.971663

$3.33531 \times 10^{-1}$     $1.25000 \times 10^{-1}$     $2.21292 \times 10^{-2}$     $9.11823 \times 10^{-1}$             3.214491             1.812495

$1.75870 \times 10^{-1}$     $6.25000 \times 10^{-2}$     $2.87214 \times 10^{-3}$     $2.57098 \times 10^{-1}$             3.190428             1.978141

Order 2, $d=3$

$h$                       $\tau$                  $L^2$ error                  $H^1$ error     eoc($L^2$ error)     eoc($H^1$ error)

$1.25000$                    $1.00000$     $3.39847 \times 10^{-1}$                    $4.82017$                  —                  —

$8.66599 \times 10^{-1}$     $5.00000 \times 10^{-1}$     $1.31777 \times 10^{-1}$                    $2.72351$             2.586204             1.558424

$6.31590 \times 10^{-1}$     $2.50000 \times 10^{-1}$     $3.56234 \times 10^{-2}$     $9.04974 \times 10^{-1}$             4.135188             3.482911

$3.33531 \times 10^{-1}$     $1.25000 \times 10^{-1}$     $2.55744 \times 10^{-3}$     $1.41605 \times 10^{-1}$             4.125257             2.905010

$1.75870 \times 10^{-1}$     $6.25000 \times 10^{-2}$     $1.69910 \times 10^{-4}$     $2.10189 \times 10^{-2}$             4.236777             2.980704

Order 3, $d=3$

Results for advection-diffusion problem for $d=2, 3$.

§ PROOF OF REGULARITY

Under the assumptions of $A4$ in Theorem <ref>, we get the additional regularity, the solution $u \in W(Z_0,H)$

It follows from the regularity assumptions on the flow $\mathbf{\Phi}_t$ that the pair $(Z_1(t), \phi_t|_{Z_1(t)})|_{t \in I}$ is compatible. We first show that there exists a $u_g \in L^2_{Z_1}$ such that:

\begin{align*}
\sum_{i = 1}^2 \int_{\Omega_i(t)}  \mathcal{A}_i(t;x) \nabla u_g \cdot \nabla v= \int_{\Gamma(t)} G v\quad \text{for almost all} \; t \in I,\; \forall v \in L^2_V,\; \norm{u_g}_{L^2_{Z_1}} \leq C\norm{G}_{L^2_{\mathcal{V}_\Gamma}}
\end{align*}

Fix $t \in I$, set $u^t_g \in Z_1(t)$ to be the solution to:

\begin{align} \label{eq:regular}
\sum_{i = 1}^2 \int_{\Omega_i(t)}  \mathcal{A}(t;x) \nabla u^t_g \cdot \nabla v  = \int_{\Gamma(t)} G v
\end{align}

For all $v \in V(t)$. Then, by Theorem 1 in [29] we have that for all $t \in I$, there exists a unique solution $u^t_g( \cdot) \in Z_1(t)$ to (<ref>).

Set $u_g(t; \cdot) := u^t_g(\cdot)$, we show that this solution is in-fact in $L^2_{Z_1}$. It suffices to show that the map $t \to (u_g, w)_{Z_1(t)}$ is measurable for all $w \in L^2_{Z_1}$ and $\norm{u_g}_{Z_1(t)}$ uniformly bounded in time (see Lemma 2.14 in [4]). We first show via this method that $u_g \in L^2_{V}$, testing the differential equation (<ref>) with $v = u_g$ and integrating in time, we have:

\begin{align*}
\norm{u_g}^2_{L^2_V} \leq C(\gamma) \norm{G}^2_{L^2_{V^*}},
\end{align*}

via Young's and Poincare's inequalities. To show the measurability, since $\mathcal{A}(t;x)$ is bounded from both below and above independent of time, we can induce the equivalent inner product $(u,v)_{\tilde{V}(t)} := (\mathcal{A}(t;x)\nabla u, \nabla v)_{H(t)}$. Showing measurability then follows as:

\begin{align*}
(u_g,v)_{\tilde{V}(t)} = \sum_{i = 1}^2 \int_{\Omega_i(t)}  \mathcal{A}(t;x) \nabla u_g \cdot \nabla v  = \int_{\Gamma(t)} G v = \int_{\Gamma(t)}  \langle G, v\rangle_{\mathcal{V}_\Gamma(t)}
\end{align*}

and since $(v,G) \in L^2_V \times L^2_{\mathcal{V}_\Gamma^*}$, by Lemma 2.13 [4], the map $t\mapsto \langle G, v\rangle_{\mathcal{V}_\Gamma(t)}$ is measurable and hence so is $(u_g,v)_{\tilde{V}(t)}$, so $u_g \in L^2_{V}$. Before moving on, note that for any fixed $t \in I$, since $u_g(t; \cdot) \in Z_1(t)$, we can integrate by parts, obtaining:

\begin{align*}
\int_{\Gamma(t)}\bigg\llbracket \mathcal{A}(t;x) \frac{\partial u}{\partial \nu_\Gamma} \bigg\rrbracket  v-\sum_{i = 1}^2 \int_{\Omega_i(t)} \nabla \cdot (\mathcal{A}_i(t;x) \nabla u_g) v = \int_{\Gamma(t)} G v
\end{align*}

Since the space $C^\infty_0(\Omega_1(t)) \times C^\infty_0(\Omega_2(t))$ is dense in $H(t)$, we get that:

\begin{align*}
\nabla \cdot (\mathcal{A}_i(t;x) \nabla u_g)  = 0 \; \text{a.e}
\end{align*}

Since $\mathcal{A}(t;x)$ is assumed to be differentiable, we introduce the equivalent inner product on $Z_1(t)$:

\begin{multline*}
(\eta,v)_{\tilde{Z}_1(t)} := \\
\sum_{i = 1}^2 \int_{\Omega_i(t)} \nabla \cdot (\mathcal{A}_i(t;x) \nabla \eta_i)  \nabla \cdot (\mathcal{A}_i(t;x) \nabla v_i) + (\nabla \mathcal{A}_i(t;x))\cdot \nabla \eta_i(\nabla \mathcal{A}_i(t;x))\cdot \nabla v_i + \nabla \eta_i \cdot \nabla v_i.
\end{multline*}

To show the equivalence of inner products $(\cdot, \cdot)_{\tilde{Z}_1(t)}$ and $(\cdot, \cdot)_{Z_1(t)}$, note:

\begin{align*}
\norm{\eta}^2_{\tilde{Z}_1(t)} &= \sum_{i = 1}^2 \int_{\Omega_i(t)} |\nabla \cdot (\mathcal{A}_i(t;x) \nabla \eta_i)|^2 + |(\nabla \mathcal{A}_i(t;x))\cdot \nabla \eta_i|^2 +|\nabla \eta_i|^2\\
&= \sum_{i = 1}^2 \int_{\Omega_i(t)}|\mathcal{A}_i(t;x) \Delta \eta_i|^2 +  2|(\nabla  \mathcal{A}_i(t;x) )\cdot \nabla \eta_i|^2  +|\nabla \eta_i|^2 \geq \norm{\eta}^2_{Z_1(t)}
\end{align*}

\begin{align*}
\norm{\eta}^2_{\tilde{Z}_1(t)} \leq C(|\mathcal{A}_i|_{C^1(\Omega_i(t))})(\norm{\Delta \eta_i}^2_{H(t)} +\norm{\nabla \eta_i}^2_{H(t)} +\norm{\eta_i}^2_{H(t)}) \leq C\norm{\eta}^2_{Z_1(t)}.
\end{align*}

Substituting in $\eta = u_g$, we arrive at:

\begin{align*}
(u_g,v)_{\tilde{Z}_1(t)}  = \sum_{i = 1}^2 \int_{\Omega_i(t)} (\nabla \mathcal{A}_i(t;x))\cdot \nabla u_i(\nabla \mathcal{A}_i(t;x))\cdot \nabla v_i +\nabla u_i \cdot \nabla v_i.
\end{align*}

Since we already know that $u_g \in L^2_V$, $u_g$ is measurable and hence the map $t \to (u_g,v)_{\tilde{Z}_1(t)}$ is measurable for all $v \in Z_1(t)$. For the bound, Theorem 1 in [29] gives us a constant $C_t$ (that depends on time) such that:

\begin{align*}
\norm{u_g}_{Z_1(t)} \leq C_t \norm{G}_{\mathcal{V}_\Gamma(t)}
\end{align*}

Using a similar method as Lemma <ref> and changing the variables, it follows that there exists $C_T> C_t$, $C_T <\infty$. Hence:

\begin{align*}
\int_0^T  \norm{u_g}_{Z_1(t)}^2 \leq C_T^2 \int_0^T \norm{G}_{\mathcal{V}_\Gamma(t)}^2 <\infty.
\end{align*}

Hence $u_g \in L^2_{Z_1}$.

We now show that $u_g \in W(Z_1, H)$. We first calculate the following commutator, for a $v \in W(V,V)$:

\begin{align*}
[\nabla, \partial_t^\bullet ]_l v = \bigg[\frac{\partial}{\partial x^l}\partial_t^\bullet - \partial_t^\bullet  \frac{\partial}{\partial x^l}\bigg] v &=\bigg[\frac{\partial}{\partial x^l}(\partial_t + \mathbf{w}\cdot \nabla)- (\partial_t + \mathbf{w}\cdot \nabla) \frac{\partial}{\partial x^l}\bigg] v, \\
&= \bigg[\frac{\partial}{\partial x^l} \mathbf{w}\cdot \nabla - \mathbf{w}\cdot \nabla \frac{\partial}{\partial x^l} \bigg]v =\sum_{j = 1}^d (\nabla_l \mathbf{w}_j)\nabla_j v.
\end{align*}

Then, if $u_g$ were differentiable, testing with some $\eta \in \mathcal{D}_V$ and taking the derivative on both sides yields:

\begin{multline*}
\lambda(t; \mathcal{A} \nabla u_g, \nabla \eta) + \sum_{i = 1}^2 \int_{\Omega_i(t)} \partial^\bullet_t (\mathcal{A}_i) \nabla u_g \cdot \nabla \eta_i + \mathcal{A}_i \nabla \partial_t^\bullet  u_g \cdot \nabla \eta_i  +D_i (\mathbf{w}, u_g, \eta)
= \lambda^\Gamma (t; G, \eta) +\langle  \partial_t^\bullet G,  \eta\rangle_{\mathcal{V}_\Gamma^*(t)},
\end{multline*}

where $D_i$ was introduced in equation <ref> and $\lambda^\Gamma(t; \cdot, \cdot)$ is the bilinear form from definition <ref> with respect to the triple $L^2_{\mathcal{V}_\Gamma} \subset L^2_{\mathcal{H}_\Gamma} \subset L^2_{\mathcal{V}^*_{\Gamma}}$ (see section 7.2 in [3] for the explicit form of $\lambda^\Gamma(t; \cdot,\cdot)$). Integrating over time:

\begin{multline*}
\sum_{i = 1}^2 \int_0^T\int_{\Omega_i(t)} \mathcal{A}_i \nabla \partial_t^\bullet  u_g \cdot \nabla \eta_i  = \int_0^T  \langle  \partial_t^\bullet G,  \eta\rangle_{\mathcal{V}_\Gamma^*(t)} - \bigg[\sum_{i = 1}^2 \int_{\Omega_i(t)}\partial^\bullet_t (\mathcal{A}_i) \nabla u_g \cdot \nabla \eta_i +   D_i (\mathbf{w}, u_g, \eta)\bigg] + \tilde{\lambda}(t; G, \eta),
\end{multline*}

\begin{align*}
\tilde{\lambda}(t;u_g,G,\eta) := \lambda^{\Gamma}(t;G,\eta)- \lambda(t; \mathcal{A}\nabla u_g, \nabla \eta),
\end{align*}

Then, set $\tilde{u}_g$ to be the solution to

\begin{align}\label{eq:regularity}
\sum_{i = 1}^2 \int_0^T\int_{\Omega_i(t)} \nabla \tilde{u}_g\cdot \nabla \eta  = \int_0^T  \int_{\Gamma(t)} \langle  \partial_t^\bullet G,  \eta\rangle_{\mathcal{V}_\Gamma^*(t)} - \bigg[\sum_{i = 1}^2 \int_{\Omega_i(t)}\partial^\bullet_t (\mathcal{A}_i) \nabla u_g \cdot \nabla \eta_i +   D_i (\mathbf{w}, u_g, \eta)\bigg] + \tilde{\lambda}(t; G, \eta).
\end{align}

Using the same method as before (showing a solution for a fixed $t \in I$ by use of [29]) we have that if $\partial^\bullet_tG \in \mathcal{V}_\Gamma^*$, there exists a unique $\tilde{u}_g \in L^2_V$ solving equation (<ref>). Set:

\begin{align*}
w := \phi_t \int_0^t \phi_{-t} \tilde{u}_g.
\end{align*}

Then $w \in W(V,V)$ and $w$ solves:

\begin{multline*}
\sum_{i = 1}^2 \int_0^T\int_{\Omega_i(t)} \mathcal{A}_i\nabla \partial_t^\bullet w\cdot \nabla \eta  \\
=\int_0^T \int_{\Gamma(t)} \langle  \partial_t^\bullet G,  \eta\rangle_{\mathcal{V}_\Gamma^*(t)} - \bigg[\sum_{i = 1}^2 \int_{\Omega_i(t)}\partial^\bullet_t (\mathcal{A}_i) \nabla u_g \cdot \nabla \eta_i +   D_i (\mathbf{w}, u_g, \eta)\bigg] + \tilde{\lambda}(t; G, \eta).
\end{multline*}

Then, subtracting the original solution $u_g$ and rearranging, we obtain:

\begin{align*}
\sum_{i = 1}^2 \int_0^T \int_{\Omega_i(t)} \mathcal{A}_i\nabla  (w - u_g)\cdot \nabla \partial_t^\bullet \eta  = \int_0^T\sum_{i = 1}^2 \int_{\Omega_i(t)} &D_i(\mathbf{w}, u_g- w, \eta)\\
&- \bigg[
(\partial_t^\bullet \mathcal{A}_i)\nabla (u_g - w)\cdot \nabla \eta +\lambda(t;\mathcal{A}\nabla(u_g-w),\nabla \eta)\bigg].
\end{align*}

Testing with $\eta = \phi_{t}\int_0^t \phi_{-s} [w-u_g](s)$, we arrive at:

\begin{multline*}
\sum_{i =1}^2 \int_0^T \int_{\Omega_i(t)}\mathcal{A}_i |\nabla (w - u_g)|^2  \\ =  \int_0^T\sum_{i = 1}^2 \int_{\Omega_i(t)} D_i(\mathbf{w}, u_g - w, \eta)- \bigg[
(\partial_t^\bullet \mathcal{A}_i)\nabla (u_g - w)\cdot \nabla \eta +\lambda(t;\mathcal{A}\nabla(u_g-w),\nabla \eta)\bigg] =: T_1 + T_2 + T_3.
\end{multline*}

\begin{align*}
T_1 = \sum_{i = 1}^2 \int_{\Omega_i(t)}    D_i(\mathbf{w}, u_g - w,\eta) = \sum_{i = 1}^2 \int_{\Omega_i(t)} \sum_{l,k = 1}^d \mathcal{A}_i(t;x)[\nabla_l \mathbf{w}_k +\nabla_k \mathbf{w}_l]\nabla_l (u_g - w) \nabla_k \eta.
\end{align*}

Let $M_\mathbf{w} := [\nabla_l \mathbf{w}_k +\nabla_k \mathbf{w}_l]_{l,k}$ be the symmetric matrix of partial derivatives of the velocity field, then:

\begin{align*}
T_1 &= \sum_{i = 1}^2 \int_{\Omega_i(t)} \mathcal{A}_i (M_{\mathbf{w}}\nabla (u_g - w))\cdot \nabla \eta, \\
&= \sum_{i = 1}^2 \int_{\Omega_i(t)}\mathcal{A}_i(M_{\mathbf{w}}\nabla (u_g - w))\cdot [\nabla \mathbf{\Phi}_i(-t;y)]^T\phi_t\nabla \int_0^t  \phi_{-s} [w-u_g](s),\\
&=\sum_{i = 1}^2\int_{\Omega_i(t)}\mathcal{A}_i(M_{\mathbf{w}}\nabla (u_g - w)) \cdot [\nabla\mathbf{\Phi}_i(-t;y)]^T\phi_t\int_0^t [\nabla \mathbf{\Phi}_i(s;y) ]^T \phi_{-s}( \nabla[w-u_g])(s),\\
&=\sum_{i = 1}^2\int_{\Omega_i(t)}\int_0^t\mathcal{A}_i(M_{\mathbf{w}}\nabla (u_g - w))\cdot [\nabla\mathbf{\Phi}_i(-t;y)]^T\phi_t [\nabla \mathbf{\Phi}_i(s;y) ]^T \phi_{-s}( \nabla[w-u_g])(s)\,dsdy.\\
\end{align*}

Using Fubini's theorem to interchange the order of integration and Holder's inequality:

\begin{align*}
T_1 &\leq C(|A|,|\nabla \mathbf{\Phi}|, |\nabla \mathbf{w}| )\norm{\nabla  (u_g - w)}_{H(t)} \\
& \qquad\qquad \cdot \int_0^t \bigg(\int_{\Omega(t)}|\nabla [\mathbf{\Phi}_i(-t;y)]^T|  \phi_t([\nabla\mathbf{\Phi}_i(s;y) ]^T \phi_{-s}(\nabla[w-u_g])(s))|^2\bigg)^{1/2}\\
&\leq C(|A_i|_{C(\mathcal{Q}_i)}, |\nabla \mathbf{\Phi}|_{C(\mathcal{Q})}, |\nabla \mathbf{w}|_{C(\mathcal{Q})} )\norm{\nabla  (u_g - w)}_{H(t)} \int_0^t \norm{\nabla(u_g - w)}_{H(s)}.
\end{align*}

Since the two remaining terms $T_2$ and $T_3$ both also depend on $\nabla \eta$, and performing a similar argument, one arrives at:

\begin{align*}
|T_2| &\leq C(|\partial^\bullet_t \mathcal{A}_i|_{C(\mathcal{Q}_i)},|\nabla \mathbf{\Phi}|_{C(\mathcal{Q})}) \norm{\nabla  (u_g - w)}_{H(t)} \int_0^t \norm{\nabla(u_g - w)}_{H(s)}, \\
|T_3| &\leq C(| \mathcal{A}_i|_{C(\mathcal{Q}_i)},|\nabla \mathbf{\Phi}|_{C(\mathcal{Q})}) \norm{\nabla  (u_g - w)}_{H(t)} \int_0^t \norm{\nabla(u_g - w)}_{H(s)},
\end{align*}

We therefore arrive at (by use of a Poincaré inequality):

\begin{align*}
\norm{u_g - w}_{V(t)} \leq C\int_0^t \norm{u_g - w}_{V(s)}.
\end{align*}

By a standard Grönwall argument, we infer that $u_g = w$ for almost all $t \in I$ and hence $u_g \in W(Z_1,V)$.

Let $z = u-u_g$, where $u$ is the weak solution from the Problem <ref>, then:

\begin{align*}
\int_0^T \langle \partial_t^\bullet z, v\rangle_{V(t)} +a(t;z, v) +\lambda(t;z,v) = \int_0^T\tilde{l}(t;v),
\end{align*}

\begin{align*}
\tilde{l}(t;v) = (f,v)_{H(t)} - (\partial^\bullet_t u_g, v)_{H(t)} - \lambda(t; u_g,v)-([\mathcal{B} - \mathbf{w}]\cdot\nabla u_g - [\mathcal{C} - \nabla\cdot\mathbf{w}]u_g, v)_{H(t)}.
\end{align*}

Since this is a functional in $L^2_H$, by Theorem 3.13 [4], $z \in W(V,H)$ and hence $u \in W(V,H)$.

[1]
Abels, H., Garcke, H., Lam, K. F., and Weber, J.
Two-phase flow with surfactants: Diffuse interface models and their
In Transport Processes at Fluidic Interfaces. Springer
International Publishing, 2017, pp. 255–270.

[2]
Adjerid, S., Chaabane, N., and Lin, T.
An immersed discontinuous finite element method for Stokes
interface problems.
Comput. Methods Appl. Mech. Engrg. 293 (2015), 170–190.

[3]
Alphonse, A., Caetano, D., Djurdjevac, A., and Elliott, C. M.
Function spaces, time derivatives and compactness for evolving
families of Banach spaces with applications to PDEs, 2021.

[4]
Alphonse, A., Elliott, C., and Stinner, B.
An abstract framework for parabolic PDEs on evolving spaces.
Portugal. Math. 72, 1 (2015), 1–46.

[5]
Alphonse, A., Elliott, C., and Stinner, B.
On some linear parabolic PDEs on moving hypersurfaces.
Interface. Free Bound. 17, 2 (2015), 157–187.

[6]
Arnold, D. N., Brezzi, F., Cockburn, B., and Marini, L. D.
Unified analysis of discontinuous Galerkin methods for elliptic
SIAM J. Numer. Anal. 39, 5 (2001/02), 1749–1779.

[7]
Balay, S., Abhyankar, S., Adams, M. F., Brown, J., Brune, P., Buschelman,
K., Dalcin, L., Eijkhout, V., Gropp, W. D., Karpeyev, D., Kaushik, D.,
Knepley, M. G., May, D. A., McInnes, L. C., Mills, R. T., Munson, T., Rupp,
K., Sanan, P., Smith, B. F., Zampini, S., Zhang, H., and Zhang, H.
PETSc users manual.
Tech. Rep. ANL-95/11 - Revision 3.11, Argonne National Laboratory,

[8]
Balay, S., Gropp, W. D., McInnes, L. C., and Smith, B. F.
Efficient management of parallelism in object oriented numerical
software libraries.
In Modern Software Tools in Scientific Computing (1997),
E. Arge, A. M. Bruaset, and H. P. Langtangen, Eds., Birkhäuser Press,
pp. 163–202.

[9]
Barrett, J. W., Garcke, H., and Nürnberg, R.
On the parametric finite element approximation of evolving
hypersurfaces in $\Bbb R^3$.
J. Comput. Phys. 227, 9 (2008), 4281–4307.

[10]
Barrett, J. W., Garcke, H., and Nürnberg, R.
Stable phase field approximations of anisotropic solidification.
IMA J. Numer. Anal. 34, 4 (2014), 1289–1327.

[11]
Barrett, J. W., Garcke, H., and Nürnberg, R.
On the stable numerical approximation of two-phase flow with
insoluble surfactant.
ESAIM Math. Model. Numer. Anal. 49, 2 (2015), 421–458.

[12]
Braess, D.
Finite elements, third ed.
Cambridge University Press, Cambridge, 2007.
Theory, fast solvers, and applications in elasticity theory,
Translated from the German by Larry L. Schumaker.

[13]
Brenner, S. C., and Scott, L. R.
The mathematical theory of finite element methods, third ed.,
vol. 15 of Texts in Applied Mathematics.
Springer, New York, 2008.

[14]
Ciarlet, P. G., and Raviart, P.-A.
The combined effect of curved boundaries and numerical integration in
isoparametric finite element methods.
In The mathematical foundations of the finite element method
with applications to partial differential equations (Proc. Sympos.,
Univ. Maryland, Baltimore, Md., 1972) (1972), Academic Press, New
York, pp. 409–474.

[15]
Dalcin, L. D., Paz, R. R., Kler, P. A., and Cosimo, A.
Parallel distributed computing using python.
Adv. Water Resour. 34, 9 (Sept. 2011), 1124–1139.
New Computational Methods and Software Tools.

[16]
Douglas, J., and Dupont, T.
Galerkin methods for parabolic equations with nonlinear boundary
Numer. Math. 20, 3 (June 1973), 213–237.

[17]
Edelmann, D.
Finite element analysis for a diffusion equation on a harmonically
evolving domain.
IMA J. Numer. Anal. 42, 2 (May 2021), 1866–1901.

[18]
Elliott, C. M., and Fritz, H.
On algorithms with good mesh properties for problems with moving
boundaries based on the harmonic map heat flow and the DeTurck trick.
SMAI Journal of Computational Mathematics 2 (2016), 141–176.

[19]
Elliott, C. M., and Fritz, H.
On approximations of the curve shortening flow and of the mean
curvature flow based on the DeTurck trick.
IMA J. Numer. Anal. 37, 2 (2017), 543–603.

[20]
Elliott, C. M., and Ranner, T.
A unified theory for continuous-in-time evolving finite element space
approximations to partial differential equations in evolving domains.
IMA J. Numer. Anal. 41, 3 (July 2021), 1696–1845.

[21]
Foote, R. L.
Regularity of the distance function.
Proc. Amer. Math. Soc. 92, 1 (1984), 153–155.

[22]
Geuzaine, C., and Remacle, J.-F.
Gmsh: A 3-d finite element mesh generator with built-in pre- and
post-processing facilities.
Int. J. Numer. Meth. Eng. 79, 11 (May 2009), 1309–1331.

[23]
Gilbarg, D., and Trudinger, N. S.
Elliptic partial differential equations of second order.
Classics in Mathematics. Springer-Verlag, Berlin, 2001.
Reprint of the 1998 edition.

[24]
Gurtin, M. E.
Toward a nonequilibrium thermodynamics of two-phase materials.
Arch. Ration. Mech. An. 100, 3 (Sept. 1988), 275–312.

[25]
Häkkinen, T. J., Sova, S. S., Corfe, I. J., Tjäderhane, L.,
Hannukainen, A., and Jernvall, J.
Modeling enamel matrix secretion in mammalian teeth.
PLOS Computational Biology 15, 5 (05 2019), 1–12.

[26]
Hartman, P.
Ordinary differential equations, vol. 38 of Classics in
Applied Mathematics.
Society for Industrial and Applied Mathematics (SIAM), Philadelphia,
PA, 2002.
Corrected reprint of the second (1982) edition [Birkhäuser,
Boston, MA; MR0658490 (83e:34002)], With a foreword by Peter Bates.

[27]
Ivančić, F., Sheu, T. W.-H., and Solovchuk, M.
Arbitrary Lagrangian {Eulerian}-type finite element methods
formulation for PDEs on time-dependent domains with vanishing discrete
space conservation law.
SIAM J. Sci. Comput. 41, 3 (Jan. 2019), A1548–A1573.

[28]
Kovács, B., and Power Guerra, C. A.
Error analysis for full discretizations of quasilinear parabolic
problems on evolving surfaces.
Numer. Meth. Part. D. E. 32, 4 (Feb. 2016), 1200–1231.

[29]
Li, H., Nistor, V., and Qiao, Y.
Uniform shift estimates for transmission problems and optimal rates
of convergence for the parametric finite element method.
Lect. Notes. Comput. Sc. (2013), 12–23.

[30]
MacDonald, G., Mackenzie, J. A., Nolan, M., and Insall, R. H.
A computational method for the coupled solution of reaction-diffusion
equations on evolving domains and manifolds: application to a model of cell
migration and chemotaxis.
J. Comput. Phys. 309 (2016), 207–226.

[31]
Mikhailov, S. E.
Traces, extensions and co-normal derivatives for elliptic systems on
Lipschitz domains.
J. Math. Anal. Appl. 378, 1 (2011), 324–342.

[32]
Mu, L., and Zhang, X.
An immersed weak galerkin method for elliptic interface problems.
J. Comput. Appl. Math. 362 (Dec. 2019), 471–483.

[33]
Prüss, J., and Simonett, G.
Moving interfaces and quasilinear parabolic evolution
equations, vol. 105 of Monographs in Mathematics.
Birkhäuser/Springer, [Cham], 2016.

[34]
Ranner, T.
firedrake moving interfaces, 2022.

[35]
Rathgeber, F., Ham, D. A., Mitchell, L., Lange, M., Luporini, F., Mcrae,
A. T. T., Bercea, G.-T., Markall, G. R., and Kelly, P. H. J.
ACM Trans. Math. Software 43, 3 (Jan. 2017), 1–27.

[36]
Roubicek, T.
Nonlinear Partial Differential Equations with Applications,
vol. 153 of International seies of numerical mathematics.
Birkhäuser-Verlag, Basel Boston Berlin, 2005.

[37]
Ryder, L. S., Dagdas, Y. F., Kershaw, M. J., Venkataraman, C., Madzvamuse,
A., Yan, X., Cruz-Mireles, N., Soanes, D. M., Oses-Ruiz, M., Styles, V.,
Sklenar, J., Menke, F. L. H., and Talbot, N. J.
A sensor kinase controls turgor-driven plant infection by the rice
blast fungus.
Nature 574, 7778 (Oct. 2019), 423–427.

[38]
San Martín, J., Smaranda, L., and Takahashi, T.
Convergence of a finite element/ALE method for the stokes equations
in a domain depending on time.
J. Comput. Appl. Math. 230, 2 (Aug. 2009), 521–545.

[39]
Schramm, L. L., Stasiuk, E. N., and Marangoni, D. G.
2 surfactants and their applications.
Annu. Rep. Prog. Chem., Sect. C: Phys. Chem. 99 (2003),

[40]
Werner, P., Burger, M., Frank, F., and Garcke, H.
A diffuse interface model for cell blebbing including membrane-cortex
coupling with linker dynamics.
SIAM J. Appl. Math. 82, 3 (June 2022), 1091–1112.

[41]
Yang, Q., and Zhang, X.
Discontinuous galerkin immersed finite element methods for parabolic
interface problems.
J. Comput. Appl. Math. 299 (June 2016), 127–139.
Recent Advances in Numerical Methods for Systems of Partial
Differential Equations.