1]Physical Sciences Division, Pacific Northwest National Laboratory, Richland,
Washington, USA, 99354

# Active Learning with Fully Bayesian Neural Networks for Discontinuous and
Nonstationary Data

Maxim Ziatdinov  [

###### Abstract

Active learning optimizes the exploration of large parameter spaces by
strategically selecting which experiments or simulations to conduct, thus
reducing resource consumption and potentially accelerating scientific
discovery. A key component of this approach is a probabilistic surrogate
model, typically a Gaussian Process (GP), which approximates an unknown
functional relationship between control parameters and a target property.
However, conventional GPs often struggle when applied to systems with
discontinuities and non-stationarities, prompting the exploration of
alternative models. This limitation becomes particularly relevant in physical
science problems, which are often characterized by abrupt transitions between
different system states and rapid changes in physical property behavior. Fully
Bayesian Neural Networks (FBNNs) serve as a promising substitute, treating all
neural network weights probabilistically and leveraging advanced Markov Chain
Monte Carlo techniques for direct sampling from the posterior distribution.
This approach enables FBNNs to provide reliable predictive distributions,
crucial for making informed decisions under uncertainty in the active learning
setting. Although traditionally considered too computationally expensive for
’big data’ applications, many physical sciences problems involve small amounts
of data in relatively low-dimensional parameter spaces. Here, we assess the
suitability and performance of FBNNs with the No-U-Turn Sampler for active
learning tasks in the ’small data’ regime, highlighting their potential to
enhance predictive accuracy and reliability on test functions relevant to
problems in physical sciences.

###### keywords:

Active learning, Bayesian neural networks, Gaussian process, Uncertainty
quantification

##  1 Introduction

Active learning [1, 2] is a machine learning technique for efficiently
exploring large parameter spaces. By strategically choosing which experiments
or simulations to perform next, active learning algorithms can significantly
reduce the number of observations needed to navigate these spaces effectively.
This not only has the potential to speed up the discovery process but also to
minimize resource use, making it a valuable approach in fields where
experimental costs are high or when the experimental setup is particularly
complex.

The underlying principle is straightforward. We identify a space of
experimental or simulation parameters that we can control. For each
combination of these parameters, we can measure an associated physical
property. We assume the existence of a function that links these parameters to
the target property, but its exact form is unknown — if it were known, machine
learning would be unnecessary. Instead, we approximate it with a probabilistic
surrogate model. Given that it’s usually impractical to explore every possible
combination of parameters due to resource constraints, and random exploration
isn’t efficient, we rely on the probabilistic predictions from this surrogate
model to intelligently guide our selection of next measurements.

Gaussian Processes (GP) are often the default choice for the surrogate model.
They are prized for their capacity to provide uncertainty quantification and
smooth interpolations, which are essential in active learning contexts where
the goal is to reconstruct unknown (”black box”) functions across a parameter
space with limited data. Standard GP assumes the target function is smooth and
continuous, and their effectiveness relies on this assumption. As a result,
they struggle with discontinuous or non-stationary functions, where abrupt
changes or variable statistical properties across the input space challenge
the model’s basic premises.

Several techniques have been developed to mitigate this limitation. For
instance, incorporating a change point model [3, 4] can effectively address
the shortcomings of GPs in modeling functions with abrupt transitions between
different states. However, this approach depends on sophisticated heuristics
and assumptions about the specific type of processes involved, making it
challenging to apply universally to arbitrary physical systems, especially in
active learning scenarios, where the initial dataset is minimal and the
underlying data behavior is largely unknown.

Introducing a prior mean function, whether probabilistic [5] or deterministic
[6], into a GP is another approach to address the above limitations. This
prior mean function acts as a baseline expectation of the system’s behavior,
guiding the GP towards more plausible values based on prior knowledge or
theoretical insights about the function being modeled. For instance, if some
basic physical laws or trends are known, these can be encoded within the prior
mean to inform the GP’s predictions, potentially improving both the efficiency
and accuracy of the learning process. However, challenges arise when the
physical model of the system is unknown or when the relationships within the
data cannot be easily expressed in a learnable (i.e. differentiable) form
suitable for integration into a GP. In such cases, specifying an appropriate
prior mean function becomes non-trivial. If the system’s behavior is highly
complex or the underlying physics is not well-understood, the use of an
incorrect prior mean can lead to biased or misleading predictions.

Figure 1: Fully Bayesian neural networks can outperform their Gaussian process
counterparts on active learning tasks, particularly for systems exhibiting
rapid changes in behavior and/or discontinuous transitions between different
states. These types of behavior are often encountered in various physical
science domains.

In the absence of prior knowledge, deep kernel learning (DKL), which combines
a neural network with a GP, emerged as a promising strategy for overcoming
some of the limitations of traditional GPs [7]. In DKL, the neural network
transforms the input data into a new feature space that is more amenable to
GP. Thus far, most DKL work has focused on problems in high-dimensional input
spaces using either Maximum a Posteriori or Maximum Likelihood Estimation
approximations [8, 9]. However, Ober et al. highlighted that standard DKL
might be prone to overfitting issues [10] and proposed that a fully Bayesian
approach — encompassing probabilistic treatments of both neural network
weights and GP hyperparameters — could mitigate these issues. However,
considering that DKL was initially conceptualized as an approximation of a
Bayesian deep learning, it raises the question of whether directly employing a
fully Bayesian neural network might be more effective.

A Fully Bayesian Neural Network (FBNN) employs a probabilistic treatment of
all neural network weights, treating them as random variables with specified
prior distributions [11, 12]. It utilizes advanced Markov Chain Monte Carlo
techniques to sample directly from the posterior distribution of its weights,
allowing the FBNN to account for all plausible weight configurations. This
approach enables the network to make probabilistic predictions, not just
single-point estimates but entire distributions of possible outcomes,
quantifying the inherent uncertainty. Despite their theoretical advantages in
providing robust uncertainty quantification, FBNNs have not been widely
adopted due to their significant computational demands [13]. Generally,
research in Bayesian deep learning has been geared towards finding suitable
approximations to FBNNs for handling large, high-dimensional data problems.
Yet, in the physical sciences, especially in the experimental research, where
parameter spaces are often relatively lower-dimensional, a fully Bayesian
approach could be feasible. Employing FBNNs for active learning in these
settings could provide reliable uncertainty estimates and robust predictive
modeling, essential for scientific applications. Despite these potential
advantages, there has been limited effort to systematically evaluate the
performance of FBNNs against benchmarks relevant to problems in the physical
sciences.

Our objective is to assess and compare the performance of FBNNs across various
active learning tasks relevant to autonomous scientific discovery. The central
question we explore is: ”Should you, as a domain scientist, consider a FBNN
for your experimental or modeling needs?” To investigate this, we have
developed test functions that are pertinent to common challenges in the
physical sciences, such as phase transitions, and used them to rigorously
assess the comparative strengths and weaknesses of FBNNs against traditional
GPs. Our results show that FBNNs either outperform or at the very least show
similar performance to traditional GPs, making them a compelling option for
practitioners exploring active learning methods in experimental or modeling
applications.

##  2 Related work

To date, active learning with BNNs typically involved various approaches to
approximate inference of BNN weights, aimed at balancing computational
efficiency with the accuracy of uncertainty estimates. These approaches
include Monte Carlo Dropout, Deep Ensembles, and Variational Inference.

Monte Carlo (MC) dropout approach uses dropout at test time to approximate
posterior sampling, providing a practical yet rough approximation of Bayesian
inference and enabling crude uncertainty estimation in neural network
predictions [14]. Y. Gal et al. [15] first showed its application for active
learning on image classification tasks. J.O. Woo further employed the MC
dropout technique with a balanced entropy acquisition method for improved
uncertainty estimates on active learning for categorical prediction tasks
[16]. In related work, Fronzi et al. demonstrated the use of MC dropout-based
active learning and density functional theory to efficiently predict bandgaps
in novel Van der Waals heterostructures [17]. However, careful calibration of
dropout rates and neural network architecture design are critical to achieve
reliable uncertainty estimates with this method [18].

Deep ensembles involve training multiple models independently and using the
variance in their predictions to estimate uncertainty [19, 20]. While not
inherently Bayesian, this method provides a practical way to gauge model
confidence and can capture model disagreement effectively. Recent studies have
shown that deep ensembles outperform MC Dropout in terms of providing more
reliable and well-calibrated uncertainty estimates, thereby improving the
active learning process for classification tasks [21, 22]. Further extending
this approach, the work by Jesús Carrete et al. [23] illustrates how deep
ensembles can be adapted to improve uncertainty estimation in computational
chemistry, enhancing active learning workflows for molecular dynamics
simulations.

Variational Inference (VI) aims at finding the best parameters for the assumed
distribution that approximates the true posterior [24]. This method is
computationally efficient but may not capture the full complexity of the
posterior distribution. A notable application of VI in neural networks is
Bayes by Backprop (BBB), which maintains a distribution (usually Gaussian)
over each weight and using backpropagation to update the parameters of this
distribution, allowing in principle learning the uncertainty in the weights.
Hafner et al. [25] showed that integrating BBB with noise contrastive priors
enhances uncertainty estimation for active learning tasks by addressing out-
of-distribution data, improving the model’s ability to select informative
samples and preventing overfitting.

The advanced MCMC methods such as Hamiltonian Monte Carlo (HMC) are considered
fully Bayesian because they attempt to sample directly from the posterior
distribution. These methods are computationally intensive but provide the most
accurate representation of uncertainty by thoroughly integrating over all
model uncertainties. Li et al. has recently investigated FBNNs utilizing HMC
sampling for Bayesian optimization demonstrating promising results across
various benchmark datasets [26]. Our focus, however, is on evaluating FBNNs
for active learning, rather than Bayesian optimization, as exploring and
understanding the overall parameter space is a more prevalent objective in
scientific research aimed at discerning how physical properties or
functionalities behave under varying conditions.

##  3 Methods

###  3.1 Fully Bayesian Neural Networks

In Bayesian neural networks, the constant weights are replaced with prior
probability distributions (Figure 2). These priors are then updated through
Bayes’ rule to derive the posterior distributions after observing data. This
modification addresses some of the limitations inherent in conventional neural
networks, such as overfitting and the inability to quantify uncertainty in
predictions. Generally, given the initial dataset $\mathcal{D}$ comprised of
inputs $x$ and targets $y$, an FBNN can be defined as a probabilistic model of
the form:

Weights: $\displaystyle\quad w\sim p(w)$ (1) Noise Variance:
$\displaystyle\quad\sigma^{2}\sim\text{Half-Normal}(\sigma_{0})$ Predictions:
$\displaystyle\quad y\sim\mathcal{N}(g(x;w),\sigma^{2})$

where $g$ is a neural network parameterized by its weights $w$, and
$\sigma^{2}$ is the observational noise variance. The idea of using prior
distributions is to represent our initial beliefs about the model parameters
before observing any data. However, in practice, defining meaningful priors
for FBNN weights is highly non-trivial due to the large number of weights
involved and the fact that these weights do not have an easily interpretable
relationship with observable features of the data. Typically, generic priors
such as normal distributions are employed, which do not necessarily convey
specific prior knowledge but are chosen for their mathematical tractability
and convenience. These priors serve as regularizers to prevent overfitting and
allow for the incorporation of uncertainty into the model.

The HMC sampling technique [27] or its variants is then employed to
approximate the posterior distribution of the weights. The predictive mean and
uncertainty at new data points $x^{*}$ are then computed as:

$\displaystyle\mu(x^{*})$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}g(x^{*},w_{\text{post}}^{(i)})$ (2)
$\displaystyle U(x^{*})$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}(y_{*}^{(i)}-\mu(x^{*}))^{2}$ (3)

where
$y_{*}^{(i)}\sim\mathcal{N}(g(x^{*},w_{\text{post}}^{(i)}),\sigma_{\text{post}}^{2(i)})$
are samples drawn from the posterior predictive distribution and N is the
total number of HMC samples.

Figure 2: Fully Bayesian neural network (FBNN). The deterministic weights of
traditional neural networks (a) are replaced with probabilistic distributions
in FBNNs (b). (c) FBNN as probabilistic graphical model: colored circles
depict observed variables while unobserved variables are in empty circles. (d)
Schematic of active learning process.

###  3.2 Gaussian Process

A Gaussian Process is considered the ”gold standard” in active learning due to
its flexibility and capacity for predictive uncertainty estimation [28]. As
such, it serves as the benchmark against which we compare FBNNs. The GP models
a distribution over functions, making it highly adaptable to a wide range of
applications.

The standard GP can be written as a probabilistic model of the form:

Kernel Hyperparameters: $\displaystyle\quad\theta\sim p(\theta)$ (4) Noise
Variance: $\displaystyle\quad\sigma^{2}\sim\text{Half-Normal}(\sigma_{0})$
Latent Function: $\displaystyle\quad
f(x)\sim\mathcal{GP}(0,k(x,x^{\prime};\theta))$ Observations:
$\displaystyle\quad
y=f(x)+\epsilon,\quad\epsilon\sim\mathcal{N}(0,\sigma^{2})$

After training, the overall predictive mean at new points $x^{*}$ is given by:

$\mu(x^{*})=\frac{1}{N}\sum_{i=1}^{N}\mu_{\text{post}}^{(i)}(x^{*})$ (5)

where $\mu_{\text{post}}$ is the predictive mean for a single HMC sample with
kernel hyperparameters $\theta$ and noise variance $\sigma^{2}$, defined using
a textbook GP formula:

$\mu_{\text{post}}(x^{*})=k_{*}^{T}(K+\sigma^{2}I)^{-1}y$ (6)

and where $k_{*}$ is the vector of covariances between $x^{*}$ and all
training points $x$, $K$ is the covariance matrix of the training points, and
$I$ is the identity matrix.

Given that each individual GP produces its own uncertainty estimation, it is
instructive to express the overall predictive uncertainty using the average of
individual predictive variances and the variance across the predictive means
of different HMC samples:

$U(x^{*})=\frac{1}{N}\sum_{i=1}^{N}\left(\Sigma_{\text{post}}^{(i)}(x^{*})+(\mu_{\text{post}}^{(i)}(x^{*})-\mu(x^{*}))^{2}\right)$
(7)

where

$\Sigma_{\text{post}}(x^{*})=k_{**}-k_{*}^{T}(K+\sigma^{2}I)^{-1}k_{*}$ (8)

and where $k_{**}$ represents the variance of the process at the new point
$x^{*}$. Note that when using the Maximum A Posteriori (MAP) estimate, which
effectively corresponds to $N=1$, the predictive mean $\mu(x^{*})$ reduces to
$\mu_{\text{post}}(x^{*})$ and the predictive uncertainty $U(x^{*})$
simplifies to $\Sigma_{\text{post}}(x^{*})$.

###  3.3 Active Learning Strategy

Active learning aims to efficiently utilize limited experimental or
computational resources by iteratively selecting the most informative data
points from a predefined parameter space. Each selected points corresponds to
a parameter (1D case) or a combination of parameters (n-D case) that the model
predicts will yield the most significant insights about the system under
study. The process begins with the selection of the overall parameter space,
which is defined based on parameters that can be directly controlled and
manipulated in experiments or simulations to affect behaviour of a system
under study.

Starting with an initial dataset, usually represented by a small number of
randomly or uniformly distributed measurements, the algorithm then selects
data points from the parameter space based on the highest uncertainty in the
prediction of surrogate model. The model is updated continuously as new data
are acquired. This process is repeated until certain conditions are met, such
as reaching a set number of measurements. This strategy ensures that resources
are focused on exploring the most potentially valuable areas of the parameter
space.

Algorithm 1 Active Learning

1:Define the overall parameter space $\mathcal{X}_{\text{domain}}$

2:Select the surrogate model type (FBNN or GP)

3:Initialize dataset $\mathcal{D}=\\{(x_{i},y_{i})\\}_{i=1}^{N}$

4:Train the initial surrogate model on $\mathcal{D}$

5:repeat

6: Compute uncertainty $U(x^{*})$ for each
$x^{*}\in\mathcal{X}_{\text{domain}}$

7:
$x_{\text{next}}\leftarrow\arg\max_{x^{*}\in\mathcal{X}_{\text{domain}}}U(x^{*})$

8: Conduct an experiment or simulation at $x_{\text{next}}$ to obtain
$y_{\text{next}}$

9: Update $\mathcal{D}$ by adding $(x_{\text{next}},y_{\text{next}})$

10: Re-train the surrogate model on updated $\mathcal{D}$

11:until a predefined stopping criterion is met (e.g. a number of
measurements)

The goal of active learning is to optimize the efficiency and effectiveness of
data acquisition in experimental or computational settings, conserving
resources and potentially accelerating the discovery process.

###  3.4 Choice of hyperparameters

We standardized the hyperparameters across all tests for GP and FBNN models to
ensure consistency in our comparisons. Specifically, for the Gaussian Process,
we utilized a Matérn 5/2 kernel with a weakly informative prior on the
lengthscale, $l\sim\text{Log-Normal}(0,1)$. For the Bayesian Neural Network,
unless explicitly specified otherwise, our architecture consisted of a three-
layer fully-connected neural network with 32, 16, and 8 neurons in each
respective layer, each activated by the hyperbolic tangent function, with
normal priors over their weights, $w\sim\text{Normal}(0,1)$. In all models,
unless explicitly specified otherwise, we adopted an observational noise prior
of $\sigma\sim\text{Half-Normal}(1)$, which is typical for scenarios where the
variance of noise is expected to be moderate but should not dominate the
signal.

###  3.5 Model training

To infer the posterior distributions over model parameters, we use the No-U-
Turn Sampler (NUTS) [29]. It builds upon the HMC approach by automatically
tuning the step size and dynamically setting the number of steps to take in
each iteration. This automated tuning aspect of NUTS is particularly useful in
autonomous science applications where manual tuning is infeasible or
suboptimal. Here, we utilized the version of NUTS from NumPyro probabilistic
programming library [30] due to its compatibility with JAX’s [31] automatic
differentiation and just-in-time compiling capabilities for significantly
faster computations. All computations were performed with double (64bit)
precision on AWS EC2 t3.2xlarge instance.

###  3.6 Test functions and datasets

We have designed our test functions (Figure 3a) to reflect common challenges
in applying active learning to physical systems that are not always addressed
by available benchmark functions.

  * •

Discontinuous Transition 1: This function models a system undergoing a
discontinuous phase transition, starting with one power-law behavior and then
abruptly switching to another at a specific point, representing a change from
one phase to another.

  * •

Discontinuous Transition 2: This function represents another potential example
of phase transition dynamics. It transitions from a sinusoidal form to a
sinusoidal with an added constant shift. This changes the baseline level of
the function, leading to a discontinuity at that transition point.

  * •

Discontinuous Transition 3: This function models a system undergoing multiple
discontinuous transitions, starting with a power-law increase, followed by
moderated growth described by a polynomial of a lower degree and a negative
offset, and concluding with exponential progression. It showcases how complex
systems can exhibit multiple distinct states as external conditions change.

  * •

Nonstationary 1: This function is modified version of the higdon function [32]
that transitions from a slowly varying oscillatory pattern to a linear trend,
marking a clear shift in both the type and the nature of the function’s
output.

  * •

Nonstationary 2: Designed to test the model’s capability to handle abrupt
changes in frequency and amplitude, this function defines a sinusoidal wave
that switches to a high-frequency oscillation modulated by an exponential
decay within a specified range.

  * •

Nonstationary 3: Combining a predictable sinusoidal base with sporadic
Gaussian spikes, this function mimics real-world phenomena where sudden,
localized changes occur within a broader steady pattern, challenging the
model’s ability to predict under varying local conditions.

  * •

Phases2D: This function generates a two-dimensional domain mimicking a complex
phase diagram. It tests the model’s ability to discern and adapt to variations
in physical and structural properties under different conditions which is
essential in fields such as materials science.

  * •

Rays2D: This function simulates the intensity of rays over a 2D domain,
mimicking the study of magnetotransport phenomena [33] in condensed matter
physics, which can produce ray-like patterns.

  * •

Ising2D: This function generates magnetization and heat capacity maps for 2D
Ising model based on [34], providing a test bed for active learning
applications in physical simulations.

A small, zero-centered normally distributed noise is added to simulate
observational noise. We start with 4 uniformly spaced points for 1D test
functions and 10 randomly distributed points for three different random seeds
for 2D test functions. For this study, we make a simplifying assumption that
the system can be measured everywhere in the parameter space at any point in
time, which may not be the case for non-reversible processes. We leave the
exploration of models under constraints of irreversible changes as potential
avenues for future research.

##  4 Results and Discussion

Figure 3: Active learning results for 1D test functions. (a) Ground truth
function values with uniformly initialized starting (’seed’) points. (b, c)
Results of active learning with Gaussian process (GP) and Bayesian neural
network (FBNN). (d) Mean squared error (MSE) as a function of active learning
steps.

The comparison of active learning results for 1D test functions is shown in
Figure 3. One can see that for the discontinuous functions, the GP tends to
smooth over the transitions, lacking the ability to capture the abrupt shifts
in the function. In contrast, FBNNs exhibit a higher flexibility, representing
the sharp transitions with greater accuracy. For the non-stationary functions,
the traditional GP, constrained by its global length scale assumption,
struggles with specific behaviors such as straight lines or rapid high-
frequency oscillations in the first two functions. Here, FBNN demonstrates
greater adaptability. For the third nonstationary function, FBNN and GP
provide comparable performance. Overall, while FBNNs may not consistently
surpass GPs in the early stages of exploration, they tend to outperform GPs in
later phases. We stress that all the FBNN results were obtained for the exact
same network architecture, the same priors, and the same NUTS parameters. In
other words, we did not tune neural network hyperparameters or the parameters
of the training process for different test functions.

Figure 4: Perfomance evaluation of GP and FBNN in terms of negative log
predictive density (NLPD) over active learning steps for the 1D test
functions.

Figure 4 examines the performance of FBNN and GP in terms of Negative Log
Predictive Density (NLPD) over active learning steps. NLPD assesses the
predictive accuracy of probabilistic models by measuring how well they predict
unseen data, specifically focusing on their confidence in those predictions
[35]. A model with a lower NLPD not only predicts more accurately but also
indicates a more precise estimation of its uncertainty. FBNNs generally
achieve lower NLPD values for both discontinuous and nonstationary function
types, indicating better predictive accuracy and confidence. Notably, the FBNN
curves tend to exhibit significant volatility, as their predictive confidence
can fluctuate while adjusting to new information during the learning process.
Despite these fluctuations, FBNNs generally stabilize at lower NLPD values, or
even if they continue to fluctuate, their running averages tend to be smaller
than those of GPs.

Figure 5: Performance evaluation in terms of MSE and NLPD at the final active
learning step across (a) different FBNN architectures, (b) Different FBNN
noise priors. The GP ’baseline’ is shown by a green bar.

The obvious question is how sensitive FBNN results are to the choice of neural
network architecture. In Figure 5a, we compared FBNN predictions (red bars) in
terms of MSE and NLPD across three different architectures, as well as against
the GP ’baseline’ (green bar) at the last active learning step. Specifically,
we assess the following FBNN architectures: a three-layer FBNN with 64, 32,
and 16 neurons in successive layers (64-32-16), a three-layer FBNN with 32,
16, and 8 neurons (32-16-8), and a two-layer FBNN with 32 neurons in each
layer (32-32). In terms of MSE, FBNNs outperform GPs in 5 out of 6 cases,
regardless of the architecture. The only exception is the nonstationary-3, for
which a two-layer FBNN underperforms compared to the GP. For NLPD, FBNNs with
different architectures either exceed or closely match the performance of GPs,
except for the discontinuous-3, where none of the models provide reliable
performance (NLPD is positive in all cases).

In addition, for a Bayesian approach, the choice of observational noise prior
can be important [25]. In Figure 5b, we analyze FBNN predictions using three
different noise priors: $\sigma\sim\text{Half-Normal}(0.1)$,
$\sigma\sim\text{Half-Normal}(1)$, and $\sigma\sim\text{Log-Normal}(0,1)$.
Despite slight variations in MSE and NLPD values, FBNNs consistently
outperform GP across all tested priors, underscoring the effectiveness of
FBNNs for the discontinuous and nonstationary functions even under different
noise assumptions. We didn’t compare different weights priors since the normal
distribution with zero mean and a standard deviation of one is already
sufficiently generic. Overall, this capability of FBNNs to deliver competitive
results even with potentially suboptimal configurations suggests that they
offer a reliable alternative to traditional GPs.

Figure 6: Active learning results for 2D test functions. (a) Ground truth
function values. (b, c) Results of active learning with GP and FBNN (d) Mean
squared error as a function of active learning steps. The reconstruction error
was averaged over 3 runs with different initialization of starting points

We next compare GP and FBNN based active learning performance for 2D parameter
spaces (Figure 6). For Phases2D, FBNN focuses on the phase boundary and avoids
spending excessive time exploring regions of the phase space where values are
nearly constant, demonstrating a more targeted and efficient exploration
strategy. GP captures the overall landscape reasonably well but fails to
precisely detect sharp phase boundaries, resulting in inefficient use of
resources by focusing too much on areas with minimal variation. For Rays2D,
the GP significantly underperforms, failing to adequately model the underlying
data structure. By contrast, FBNN more accurately represents the pattern,
clearly delineating the rays and capturing most of their intrinsic structure.

Similarly, in the 2D Ising model, FBNN excels at delineating the phase
boundaries for magnetization and effectively focuses on regions of the
parameter space where the heat capacity is maximized. This clearly
demonstrates its ability to identify key areas based purely on predictive
uncertainty, without the need to introduce a tailored acquisition function as
it was done for GPs in [34]. We emphasize that we used the same FBNN
hyperparameters across all 2D test functions, and they were the same as those
used for the 1D test functions in Figure 3.

Given that FBNNs are at least several times more computationally expensive
than GPs, and that GPs often capture the overall landscape quite well, one
might consider whether it makes sense to perform active learning with GPs to
obtain the overall landscape and then use a ”one-shot” learning approach with
FBNN to reconstruct finer details. We note, however, that a better delineation
of transition and nonstationary regions by FBNNs is due to their focused
active learning along those regions, where FBNN predictions tend to exhibit
high uncertainty. This prompts the algorithm to allocate more resources to
these areas. In contrast, the sampling strategies employed by GPs (or those
based on random sampling) typically do not focus sufficiently on these
regions. Consequently, a ”one-shot” FBNN applied to data points collected by
GP is unlikely to achieve the level of detail and accuracy that an FBNN-driven
active learning process can provide. However, it would be interesting to
explore if partially stochastic strategies [36] for FBNN training could reduce
computational cost and time without sacrificing predictive accuracy.

##  5 Conclusions

To summarize, we explored the performance of Fully Bayesian Neural Networks
and Gaussian Processes in active learning contexts relevant to physical
sciences problems. Our findings reveal that FBNNs outperform GPs in handling
abrupt transitions and non-stationary behaviors, demonstrating robust
adaptability and more accurate uncertainty quantification. Furthermore, FBNNs
consistently demonstrated a competitive edge over GPs across multiple tested
architectures and priors, highlighting their potential as a practical
replacement for GPs, as they did not require extensive hyperparameter tuning.
Given the computational feasibility of FBNNs in lower-dimensional spaces and
the demonstrated superior performance, they emerge as a compelling choice for
active learning frameworks in autonomous science within ”small data” regimes.
However, further exploration and application in experimental and modeling
tasks are warranted, particularly to assess their performance in higher-
dimensional spaces.

## Acknowledgment

This work was supported by the Laboratory Directed Research and Development
Program at Pacific Northwest National Laboratory, a multiprogram national
laboratory operated by Battelle for the U.S. Department of Energy.

## Code and Data Availability

The source code for BNN, GP, and test functions is available at
https://github.com/ziatdinovmax/NeuroBayes.

## References

  *   * Cohn et al. [1996] Cohn, D.A., Ghahramani, Z., Jordan, M.I.: Active learning with statistical models. Journal of artificial intelligence research 4, 129–145 (1996) https://doi.org/10.1613/jair.295
  * Settles [2009] Settles, B.: Active learning literature survey. University of Wisconsin-Madison Department of Computer Sciences (2009). http://digital.library.wisc.edu/1793/60660
  * Saatçi et al. [2010] Saatçi, Y., Turner, R.D., Rasmussen, C.E.: Gaussian process change point models. In: Proceedings of the 27th International Conference on Machine Learning (ICML-10), pp. 927–934 (2010)
  * Caldarelli et al. [2022] Caldarelli, E., Wenk, P., Bauer, S., Krause, A.: Adaptive Gaussian Process Change Point Detection. PMLR (2022). https://proceedings.mlr.press/v162/caldarelli22a.html
  * Ziatdinov et al. [2022] Ziatdinov, M.A., Liu, Y., Morozovska, A.N., Eliseev, E.A., Zhang, X., Takeuchi, I., Kalinin, S.V.: Hypothesis learning in automated experiment: Application to combinatorial materials libraries. Advanced Materials 34(20), 2201345 (2022) https://doi.org/10.1002/adma.202201345
  * Noack et al. [2021] Noack, M.M., Zwart, P.H., Ushizima, D.M., Fukuto, M., Yager, K.G., Elbert, K.C., Murray, C.B., Stein, A., Doerk, G.S., Tsai, E.H., et al.: Gaussian processes for autonomous data acquisition at large-scale synchrotron and neutron facilities. Nature Reviews Physics 3(10), 685–697 (2021) https://doi.org/10.1038/s42254-021-00345-y
  * Wilson et al. [2015] Wilson, A.G., Hu, Z., Salakhutdinov, R., Xing, E.P.: Deep kernel learning (2015) arXiv:1511.02222 [cs.LG]
  * Liu et al. [2022] Liu, Y., Kelley, K.P., Vasudevan, R.K., Funakubo, H., Ziatdinov, M.A., Kalinin, S.V.: Experimental discovery of structure–property relationships in ferroelectric materials via active learning. Nature Machine Intelligence 4(4), 341–350 (2022) https://doi.org/10.1038/s42256-022-00460-0
  * Sivaraman and Jackson [2022] Sivaraman, G., Jackson, N.E.: Coarse-grained density functional theory predictions via deep kernel learning. Journal of Chemical Theory and Computation 18(2), 1129–1141 (2022) https://doi.org/10.1021/acs.jctc.1c01001
  * Ober et al. [2021] Ober, S.W., Rasmussen, C.E., Wilk, M.: The promises and pitfalls of deep kernel learning (2021) arXiv:2102.12108 [stat.ML]
  * Neal [1996] Neal, R.M.: Bayesian Learning for Neural Networks. Springer, Berlin, Heidelberg (1996)
  * Andrieu et al. [1999] Andrieu, C., Freitas, J.a., Doucet, A.: Robust full bayesian methods for neural networks. In: Solla, S., Leen, T., Müller, K. (eds.) Advances in Neural Information Processing Systems, vol. 12. MIT Press, Cambridge, MA (1999)
  * Izmailov et al. [2021] Izmailov, P., Vikram, S., Hoffman, M.D., Wilson, A.G.G.: What are bayesian neural network posteriors really like? In: International Conference on Machine Learning, pp. 4629–4640 (2021). https://doi.org/10.48550/arXiv.2104.14421
  * Gal and Ghahramani [2016] Gal, Y., Ghahramani, Z.: Dropout as a bayesian approximation: Representing model uncertainty in deep learning. International Conference on Machine Learning, 1050–1059 (2016). PMLR
  * Gal et al. [2017] Gal, Y., Islam, R., Ghahramani, Z.: Deep bayesian active learning with image data (2017) arXiv:1703.02910 [cs.LG]
  * Woo [2023] Woo, J.O.: Active learning in bayesian neural networks with balanced entropy learning principle (2023) arXiv:2105.14559 [cs.LG]
  * Fronzi et al. [2021] Fronzi, M., Isayev, O., Winkler, D.A., Shapter, J.G., Ellis, A.V., Sherrell, P.C., Shepelin, N.A., Corletto, A., Ford, M.J.: Active learning in bayesian neural networks for bandgap predictions of novel van der waals heterostructures. Advanced Intelligent Systems 3(11), 2100080 (2021) https://doi.org/10.1002/aisy.202100080
  * Verdoja and Kyrki [2020] Verdoja, F., Kyrki, V.: Notes on the behavior of mc dropout. arXiv preprint arXiv:2008.02627 (2020) https://doi.org/10.48550/arXiv.2008.02627
  * Lakshminarayanan et al. [2017] Lakshminarayanan, B., Pritzel, A., Blundell, C.: Simple and scalable predictive uncertainty estimation using deep ensembles. Advances in Neural Information Processing Systems 30 (2017) https://doi.org/10.48550/arXiv.1612.01474
  * Wilson and Izmailov [2020] Wilson, A.G., Izmailov, P.: Bayesian deep learning and a probabilistic perspective of generalization. In: Proceedings of the 34th International Conference on Neural Information Processing Systems. NIPS ’20. Curran Associates Inc., Red Hook, NY, USA (2020)
  * Beluch et al. [2018] Beluch, W.H., Genewein, T., Nurnberger, A., Kohler, J.M.: The power of ensembles for active learning in image classification. In: 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 9368–9377 (2018). https://doi.org/10.1109/CVPR.2018.00976
  * Pop and Fulop [2018] Pop, R., Fulop, P.: Deep ensemble bayesian active learning : Addressing the mode collapse issue in monte carlo dropout via ensembles (2018) 1811.03897 [cs.LG]
  * Carrete et al. [2023] Carrete, J., Montes-Campos, H., Wanzenböck, R., Heid, E., Madsen, G.K.H.: Deep ensembles vs committees for uncertainty estimation in neural-network force fields: Comparison and application to active learning. The Journal of Chemical Physics 158(20), 204801 (2023) https://doi.org/10.1063/5.0146905
  * David M. Blei and McAuliffe [2017] David M. Blei, A.K., McAuliffe, J.D.: Variational inference: A review for statisticians. Journal of the American Statistical Association 112(518), 859–877 (2017) https://doi.org/10.1080/01621459.2017.1285773
  * Hafner et al. [2019] Hafner, D., Tran, D., Lillicrap, T., Irpan, A., Davidson, J.: Noise contrastive priors for functional uncertainty (2019) arXiv:1807.09289 [stat.ML]
  * Li et al. [2023] Li, Y.L., Rudner, T.G.J., Wilson, A.G.: A study of bayesian neural network surrogates for bayesian optimization (2023) arXiv:2305.20028 [cs.LG]
  * Betancourt [2018] Betancourt, M.: A conceptual introduction to hamiltonian monte carlo (2018) arXiv:1701.02434 [stat.ME]
  * Rasmussen and Williams [2005] Rasmussen, C.E., Williams, C.K.I.: Gaussian Processes for Machine Learning. The MIT Press, Cambridge, MA, USA (2005). https://doi.org/10.7551/mitpress/3206.001.0001
  * Hoffman and Gelman [2011] Hoffman, M.D., Gelman, A.: The no-u-turn sampler: Adaptively setting path lengths in hamiltonian monte carlo (2011) arXiv:1111.4246 [stat.CO]
  * Phan et al. [2019] Phan, D., Pradhan, N., Jankowiak, M.: Composable effects for flexible and accelerated probabilistic programming in numpyro. arXiv preprint arXiv:1912.11554 (2019) 1912.11554 [stat.ML]
  * Bradbury et al. [2018] Bradbury, J., Frostig, R., Hawkins, P., Johnson, M.J., Leary, C., Maclaurin, D., Necula, G., Paszke, A., VanderPlas, J., Wanderman-Milne, S., Zhang, Q.: JAX: Composable Transformations of Python+NumPy programs. http://github.com/google/jax
  * Gramacy and Lee [2009] Gramacy, R.B., Lee, H.K.H.: Adaptive design and analysis of supercomputer experiments (2009) arXiv:0805.4359 [stat.AP]
  * Qiu et al. [2020] Qiu, G., Niu, C., Wang, Y., Si, M., Zhang, Z., Wu, W., Ye, P.D.: Quantum hall effect of weyl fermions in n-type semiconducting tellurene. Nature Nanotechnology 15(7), 585–591 (2020) https://doi.org/10.1038/s41565-020-0715-4
  * Kalinin et al. [2020] Kalinin, S.V., Valleti, M., Vasudevan, R.K., Ziatdinov, M.: Exploration of lattice Hamiltonians for functional and structural discovery via Gaussian process-based exploration–exploitation. Journal of Applied Physics 128(16), 164304 (2020) https://doi.org/10.1063/5.0021762
  * Candela et al. [2003] Candela, J.Q., Girard, A., Larsen, J., Rasmussen, C.E.: Propagation of uncertainty in bayesian kernel models - application to multiple-step ahead forecasting. 2003 IEEE International Conference on Acoustics, Speech, and Signal Processing 2, 701 (2003) https://doi.org/10.1109/ICASSP.2003.1202463
  * Sharma et al. [2023] Sharma, M., Farquhar, S., Nalisnick, E., Rainforth, T.: Do Bayesian Neural Networks Need To Be Fully Stochastic? PMLR (2023). https://proceedings.mlr.press/v206/sharma23a.html


(4.5, 1) ellipse [y radius=0.75*1.8pt, x radius=0.75pt];
(4.887298, 1) ellipse [y radius=0.75*1.8pt, x radius=0.75pt];

Discrete regions in the definition of the ghost-penalty of the CG method.

The black circles represent the  of the discrete function space for the first time step.

In order to introduce the discrete function spaces next, we pose an assumption on the basis of the function space $\mathcal{P}^{\ktsol}([t_{n-1},t_n])$ contained in $\bgVhnks{n}{k}$:

The basis $\{ p_{0},\dots,p_{\ktsol} \}$ of the space $\mathcal{P}^{\ktsol}([t_{n-1},t_n])$ contains two functions $p_{i_\text{lo}}(t)$ and $p_{i_\text{up}}(t)$ which can be associated to the lower and upper domain boundary $t_{n-1}$ and $t_n$, s.t. $p_{i}(t_{n-1}) = \delta_{i_\text{lo}i}$ and $p_{i}(t_{n}) = \delta_{ i_\text{up}i}$
$\forall i=0, \dots, \ktsol$.

The discrete trial function space is then defined as follows in terms of the nomenclature known from the DG method: $U_h^n = U_{h,0}^n + u_{\text{init}}$ with

\begin{align*}
U_{h,0}^{n} := \{ v \in &\Whcut{\ksol}{ \EQlinn \cup (\EplOmn \times \{ t_n \})} \circ (\Theta_h^{n, \text{st}})^{-1} \mid v^{n-1}_+ = 0 \} \vspace*{-0.1cm}
\end{align*}

and $u_{\text{init}}(x,t) = \Piu (u^{n-1}_-(x)) \cdot p_{i_{\text{lo}}}(t) $
where $\Piu (u^{n-1}_-)$ serves as initial data for the solution as for the DG method, but with a strong imposition.

The relevant discrete region for each time slice is given as $\EQlinn \cup (\EplOmn \times \{ t_n \})$, where the expression $\EplOmn \times \{ t_n \}$ will allow for the discrete extension motivated above at the end of each time slice. The  from $\bgVhks{\ksol}$ related to the region $\EplOmn \times \{ t_n \}$ will stem from $p_{i_{\text{up}}}$, so the construction is well-posed by <ref>. The structure of the domain is illustrated in <ref> by black circles.

For the space of test functions, continuity beyond time slice boundaries is not required, so we define:

\begin{align*}
V_h^{n} := \Whcut{(\khsol,\ktsol-1)}{  \EQlinn \cup (\EplOmn \times \{ t_n \})} \circ (\Theta_h^{n, \text{st}})^{-1}.\vspace*{-0.1cm}
\end{align*}

Note that the polynomial order in time is one less as usual with CG methods.

Now we have all tools in place to state the variational formulation of the CG method within one time slab. It reads: Find $u = u_0 + u_{\text{init}}$ with $u_0 \in U_{h,0}^{n}$ s.t.

\begin{equation*}
B^{n}(u_0,v) + J^n(u_0,v) + j_h(\mathcal{F}_R^{n,+},t_n; (u_0)_-^n, v_-^n)  = f^n(v) - B^{n}(u_{\text{init}},v) \quad \forall v \in V_h^{n}. \vspace*{-0.1cm}
\end{equation*}

Here, both ghost penalty parts only act on the   of $u_0$ that needs to be solved for, but not on $u_{\text{init}}$. We further note that the decomposition $u = u_0 + u_{\text{init}}$ and hence the solution $u$ depends on definition of the basis function $p_{i_{\text{lo}}}$.

There exists a variant of the suggested CG method which we denote as CG$\square$, where only the ghost penalty stabilisation $J^{n\square}(u,v) = \int_{t_{n-1}}^{t_n} $ $j_h(\mathcal{F}_R^{n,+},t; u, v)$ $\mathrm{d}t$ is applied, which is defined on the smallest tensor-product domain which contains both domains used before.

This can be regarded as an implementational simplification, although more  are needed as elaborated on in <Ref>. Further details about this variant are given in the supplementary material.

§.§ Galerkin-Collocation method

The GCC method builds on the same idea as the CG method, but extends the strategy by one step: In addition to the continuity of the discrete solution function, we ask for continuity also in higher order time derivatives up to order $\ktreg$.

Correspondingly, we strengthen the assumption of <ref> w.r.t. to the basis of $\mathcal{P}^{\ktsol}(\In{n})$:

The basis $\{ p_{1},\dots,p_{\ktsol+1} \}$ of the space $\mathcal{P}^{\ktsol}([t_{n-1},t_n])$ contains $\ktreg +1$ functions $p_{i_\ell}(t)$ for $\ell=0,\dots,\ktreg$ which can be associated to the $\ell$th derivative at the boundary $t_{n-1}$: $\forall \ell=0,\dots,\ktreg, i=1,\dots,\ktsol+1, ~ \partial^{\ell} p_i (t_{n-1}) = \delta_{i_\ell i}$.

We start by introducing a corresponding trial function space: $U_h^n = U_{h,0}^n + u_{\text{init}}$ with

\begin{align*}
U_{h,0}^{n} \!\!:= \{ v \in &\Whcut{\ksol}{ \EplOmn \!\times\! I_n} \circ\! (\Theta_h^{n, \text{st}})^{-1} \!\mid
( \partial_t^{\ell}  v )^{n-1}_+ = 0, \ell = 0,\dots,\ktreg\},
\end{align*}

and $u_{\text{init}}(x,t) =\sum_{\ell=0}^{\ktreg}     \Piu (\partial_t^{\ell} u_-^{n-1}(x)) \cdot p_{i_\ell}(t)$
$\Piu (\partial_t^\ell u^{n-1}_-)$ serve as initial data for the solution as for the CG and the DG method. We note however that we have to prescribe $\ktreg+1$ functions as initial data here.
For the space of test functions, we now pick only order $\kttest\leq k_t -\ktreg -1$:

\begin{align*}
V_h^{n} := \Whcut{(\khsol,\kttest)}{\EplOmn \times I_n} \circ (\Theta_h^{n, \text{st}})^{-1}.
\end{align*}

For the ghost penalty we use the tensor-product variant CG$\square$ from <ref>[One could also construct GCC methods which generalise the original CG method and not CG$\square$ regarding the structure of function space and stabilisation. However, this would necessitate further extensions/ stabilisations in time derivatives of $u$, which is why we opt for the presented variant.].  The formulation of the GCC method then reads: Find $u = u_0 + u_{\text{init}}$ with $u_0 \in U_{h,0}^{n}$ s.t.

\begin{align*}
B^{n}(u,v) + J^{n\square}(u_0,v) &= f^n(v) \quad \forall v \in V_h^{n} \text{  and  } \\
B^{n}_{\text{col}}(t_l^n;u,v) + j_h(\mathcal{F}_R^{n,+},t_l^n;u_0,v) &= f^n_{\text{col}}(t_l^n;v) \quad \forall v \in \Vhks{\khsol}, l=1,..,\ktcol,\\%~ n=1,..,N,\\
% \end{align*}
% \begin{align*}
\text{with } B^n_{\text{col}}(t;u,v) := (\partial_t u + \vect{w}\cdot \nabla u, &v)_{\Oht{t}} + (\nabla u, \nabla v)_{\Oht{t}},~f^n_{\text{col}}(t;v) := (f,v)_{\Oht{t}},  %\\
% \text{and } J^{\text{col}}_n(t;u,v) &:= \sum_{F \in \mathcal{F}_R^{n,+}} \int_{\omega_F^h(t)} \frac{\gamma}{h^2} \jump{u}_{\omega_F^h(t)} \jump{v}_{\omega_F^h(t)} \mathrm{d}x.
\end{align*}

with collocation points $\{t^n_1,\dots,t^n_{\ktcol}\}$. Here, similar to the CG method, we made sure that the ghost penalty only acts on the  of $u_0$. Again, the solution $u$ will depend on the concrete choice of the basis functions $p_{i_\ell},~\ell=0,..,\ktreg$.

A possible choice for the triple $(\ktsol,\ktreg,\ktcol$) is $(3,1,1)$ which we also consider in the numerical experiments below. The one collocation point in this setting is chosen as the end point of each time interval, $t^n_1 = t^n_{\ktcol} = t^n$.

§ IMPLEMENTATION OF NUMERICAL INTEGRATION IN TIME

The discrete (bi)-linear forms presented in the previous section involve summands defined on the space-time domains $\Qhn$, c.f. e.g. <ref>. Hence, a numerical integration procedure is needed for those domains. More specifically, as $\Qhn$ was defined as the image of $\Qlinn$, cf. <ref>, and finite element computations are done element-wise, we apply the following calculation:

\begin{align}
(f, v)_{\Qhn} &= \int_{\Qhn} f v \, \mathrm{d}x = \int_{\Qlinn} |\det{D\Thehst}| (f \circ \Thehst) \hat v \, \mathrm{d}\hat x, \quad v = \hat v \circ \Thehst \\
&= \sum_{T \in \mathcal{T}_h} \int_{\Qlinn \cap (T \times I_n) } |\det{D\Thehst}| (f \circ \Thehst) \hat v \, \mathrm{d} \hat x
\end{align}

For the other summands, similar equations hold. Hence, we only need to implement an algorithm to approximate

\begin{equation}
\int_{\Qlinn \cap (T \times I_n) }  \tilde f(x,t) \mathrm{d}(\hat x,t) = \int_{t_{n-1}}^{t_n}  \int_{\Omlin{t} \cap T} \, \tilde f(\hat x,t) \, \mathrm{d}\hat x \,\mathrm{d}t \,\label{taskoftint}
\end{equation}

for each $T \in \Th$ for some integrand $\tilde f$. For a given time $\hat t$, an integration rule on $\Omlin{\,\hat t\,} \cap T$ is assumed to be available, e.g. on simplices by tessellation. Hence, the inner integral in <ref> is well-suited for numerical approximation, and it remains to replace the integration $\int_{t_{n-1}}^{t_n}$ by a numerical quadrature. In the following, we will define two variants of doing so, where the first approach is mentioned for comparison with existing literature [18, 38] and the second is our suggested approach.

In the topology-insensitive time integration, $\int_{t_{n-1}}^{t_n}$ is approximated by Gaussian quadrature. Let $(t_i, \omega_i)_i$ be the Gaussian quadrature rule on $[t_{n-1},t_n]$ of order $k_{t,int}$. Then, <ref> will be approximated by

\begin{equation}
\sum_i \omega_i \int_{\Omlin{t_i} \cap T} \, f(\hat x,t_i) \, \mathrm{d}\hat x.
\end{equation}

Note however, that the topology of $\Omlin{t} \cap T$ might change over $t \in I_n$; for instance the domain boundary can leave or enter an element. This will cause the integrand $\int_{\Omlin{t} \cap T} \, f(\hat x,t)\, \mathrm{d}\hat x$ to show a discontinuity in the first derivative even if  $f$ is smooth, which impedes the accuracy of the topology-insensitive time integration. To circumvent this, we define a time integration that takes the cut topology changes of $\Qlinn$ into account.

To approximate <ref>, we start by calculating the set $\mathcal{R}^\ast$ of subintervals of $I_n$ represented by a tuple $(t_a,t_b)$ in which no topology changes take place with <ref>.

Subdivision of the time interval

[1]
Set $\mathcal{R} = \{ t_{n-1}, t_n\} $ (Set of time points with topology changes).
Let $V$ be the set of vertices of $T$.
$v \in V$
Define the scalar function $\phi_v: I_n \to \mathbb{R}, t \mapsto \phi^{\text{lin}} (v, t)$.
Search for the roots $\mathcal{R}_v$ of $\phi_v$.
Set $\mathcal{R} \leftarrow \mathcal{R} \cup \mathcal{R}_v$.
Define $\mathcal{R}^\ast$ as the set of intervals with endpoints according to $\mathcal{R}$.

Denote by $I_n^s = (t_l, t_u) \in \mathcal{R}^\ast$ a subinterval in $\mathcal{R}^\ast$.
Then, <ref> is approximated by the topology-preserving time integration as

\begin{equation}
\sum_{I_n^s \in \mathcal{R}^\ast} \sum_i \omega_i^s \int_{\Omlin{t_i} \cap T} \, f(\hat x,t_l^s) \, \mathrm{d}\hat x,
\end{equation}

where $(t_i^s,\omega_i^s)$ is the Gaussian integration rule transformed to the interval $I_n^s$.

This time integration procedure exploits the fact that a (multi-)linear-in-space level set function is fully described by its function values on vertices of the finite element $T$. The functions $\phi_v$ within <ref> will in general be of higher polynomial order.

The topology-preserving time integration is our method of choice and for the numerical studies in this work we choose integration rules in time of exactness degree $2 (k_t + 1)$. In <ref>, we show the quadrature points generated by both approaches for simple examples in one and two space dimensions. Later in this paper, we investigate numerically differences to the topology-insensitive time integration in <Ref>.

[xmin=0,xmax=1,ymin=0,ymax=1, ylabel=$t$, xlabel=$x$, legend entries =$t_i$, legend style=at=(0.5,-0.2),anchor=north,legend columns=1, draw=none]
[dashed] 0.330009; [dashed] 0.0694318;
[dashed] 1 - 0.330009; [dashed] 1 - 0.0694318;
[Set1-A,ultra thick] table num_exp/visu/interf_out_172.dat;
[Set1-B, only marks, mark options=scale=1.5] table num_exp/visu/rule_out_172_naive.dat;

[xmin=0,xmax=1,ymin=0,ymax=1, yticklabels=,,, xlabel=$x$, legend entries =$\mathcal{R}^\ast$, legend style=at=(0.5,-0.2),anchor=north,legend columns=1, draw=none]
[dashdotted] 0.27; [dashdotted] 0.52;
[Set1-A,ultra thick] table num_exp/visu/interf_out_172.dat;
[Set1-B, only marks, mark options=scale=1.5] table num_exp/visu/rule_out_172_tsint.dat;

(A1) at (-0.50, 0.50, 0.00);
(B1) at ( 0.50, 0.50, 0.00);
(C1) at ( 0.00, 0.50,-0.80);

(A2) at (-0.50, 1.50, 0.00);
(B2) at ( 0.50, 1.50, 0.00);
(C2) at ( 0.00, 1.50,-0.80);

(A3) at (-0.50, 1.05, 0.00);
(B3) at ( 0.50, 1.05, 0.00);
(C3) at ( 0.00, 1.05,-0.80);

(R1) at ($(A2)!0.5!(C2)$);
(R2) at ($(B2)!0.25!(C2)$);
(R3) at ($(A1)!0.4!(A2)$);
(R4) at ($(B1)!1.0!(B3)$);

(S3) at ($(A3)!0.31!(C3)$);

[left,scale=0.7] at (A3) $t^2_\ast$;
[left,scale=0.7] at (R3) $t^1_\ast$;
[left,scale=0.7] at (A1) $t^{n\!-\!1}$;
[left,scale=0.7] at (A2) $t^{n}$;

[thin,densely dotted, opacity=0.5] (B1) – (C1);
[thin,densely dotted, opacity=0.5] (C1) – (A1);

[thin,densely dotted, opacity=0.5] (C1) – (C2);

[thin,densely dotted] (B3) – (C3);
[thin,densely dotted] (C3) – (A3);

[thin,fill=Set1-C,opacity=0.25] (S3) – (B3) – (C3) – cycle;
[thin,fill=Set1-D,opacity=0.25] (S3) – (A3) – (B3) – cycle;

[thin,fill=Set1-C,opacity=0.25] (R1) – (R2) – (C2) – cycle;
[thin,fill=Set1-D,opacity=0.25] (R1) – (R2) – (B2) – (A2) – cycle;

[thin,fill=Set1-C,opacity=0.25] (A1) – (B1) – (C1) – cycle;

[fill=Set1-A,opacity=0.1] (R1) – (R2) to[in=105,out=-85] (R4) to[in=5,out=-165] (R3) to[in=-100,out=33] (R1);
[draw=Set1-A!80!black,fill=none] (R1) – (R2) to[in=105,out=-85] (R4) to[in=5,out=-165] (R3) to[in=-100,out=33] (R1);

[thin] (A3) – (B3);

[Set1-A!60!Set1-C,thick] (B3) – (S3);

[thick] (A1) – (B1);

[thick] (A2) – (B2);
[thick] (B2) – (C2);
[thick] (C2) – (A2);

[thick] (A1) – (A2);
[thick] (B1) – (B2);

(A3B3) at ($(A3)!0.5!(B3)$);
(A2B2) at ($(A2)!0.5!(B2)$);
(A1B1) at ($(A1)!0.5!(B1)$);

(q) at ($(R1)!0.5!(R2)!0.5!(C2)$);

[right=0.5cm,scale=0.7] (ql) at (q) $\Omlin{t_i} \cap T$;
[->] (ql) to[in=50,out=150] (q);

(p) at ($(R1)!0.5!(R2)!0.1!(B2)!0.5!(A2)$);

[below=0.15cm,scale=0.7] at (A1B1) $T \times I_n$;

Comparison of the topology-insensitive time integration in 1+1D (left) with the topology-preserving time integration in 1+1D (center) and 2+1D (right). In the left and middle picture, blue dots represent integration points, and the red line represents the space-time interface. In the left picture, vertical lines indicate the Gaussian quadrature points (in time), which are independent of the cut configuration. In the middle, vertical lines represent the subdivision of the time interval into subintervals without topology change, on which the Gaussian rules (in time) are then applied. On each quadrature node in time a quadrature rule in space on the cut domain is applied resulting in the space-time integration points indicated by blue dots.

§ NUMERICAL INVESTIGATIONS

In this section, we want to demonstrate the numerical properties of the proposed methods. We consider three problem geometries to investigate specific aspects, each in one of the following subsections.
Recall that the notation $\ksol=r$ means $\khsol\!=\!\ktsol\!=\!r$ and analogously $\klset=r$ means $\khlset\!=\!\ktlset\!=\!r$.

All experiments are performed with  [24], an unfitted finite element extension of . Reproduction data are available at <https://gitlab.gwdg.de/fabian.heimann/repro-ho-unf-space-time-fem>.

The linear systems are solved with the direct solvers  and  of IntelMKL. As a basis for $\mathcal{P}^{k_t}([t_{n-1},t_n])$, we implement a Lagrange basis with the Gauß-Lobatto points between $t_{n-1}$ and $t_n$ of order $k_t$ as Lagrange points for DG and CG. These satisfy <ref>. For GCC of order 3, we implement cubic Hermite spline polynomials, which satisfy <ref>.

§.§ Kite geometry

The first geometry is a circle, which deforms with time into a kite-shaped object in two spatial dimensions. It is described in terms of the following level set function $\phi$

\begin{equation*}
\rho(t, y) = (1 - y^2) \cdot t, \quad r = \sqrt{(x- \rho)^2 + y^2}, \quad \phi = r - r_0, \quad r_0 = 1,
\end{equation*}

where $\tilde \Omega = [-3.5, 3.5] \times [-1.5,1.5]$, $\vect{w} = (\frac{\partial \rho}{\partial t}, 0)$, $T=0.5$. The geometry as well as some discrete regions for the DG method are shown in <ref>. The right-hand-side is calculated such that $u = \cos (\pi \frac{r}{r_0}) \cdot \sin (\pi t)$, and the stabilisation constant is chosen as $\gamma_J = 0.05$. We use unstructured simplicial meshes with mesh size $h=0.9 \cdot 0.5^{i_s}$, $\Delta t = 0.5 \cdot 2^{-i_t-1}$. In the CG and GCC methods, the extension factor is $\epsilon_f = 1.1$.

§.§.§ Geometry approximation quality

In a first experiment, we investigate the quality of the discrete domain approximation by measuring the distance to the exact interface, $\max_{t \in [0,T]} \mathrm{dist}(\partial \Omega^h(t), \partial \Omega(t))$, for the kite geometry with $\ksol=\klset$ and refinements in space and time, $i=i_t=i_s$. As can be seen in <ref>, the discrete regions satisfy the asymptotic error bound given in <ref>.

§.§.§ Convergence of methods of matching space and time order

First, we evaluate convergence of the numerical error for $i_s=i_t=i$, $\ksol=\klset$. The range of parameters $k \in \{1,2,\dots,6\}$ is applied to the DG and CG methods, whilst the GCC method is investigated for $k= 3$. In <ref> (left) and <ref>, we display the numerical errors for simultaneous space-time refinements $i=1,2,\dots$.

[ xlabel=$i$, ylabel=$\| u_h(T) - u (T)\|_{L^2(\Omega(T))}$, legend entries = $k=1$, $k=2$, $k=3$, $k=4$, $k=5$,$k=6$

, legend style=anchor=north,legend columns=1, draw=none, legend cell align=left, fill=none, legend pos =south west,
x label style=at=(axis description cs:0.5,0.15),anchor=center,
ymax=3, ymin=3e-11
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_DG_ks1_kt1_both_nref8_gamma0.05_sm2.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_DG_ks2_kt2_both_nref8_gamma0.05_sm2.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_DG_ks3_kt3_both_nref8_gamma0.05_sm2.dat;
+[thick,mark options=scale=1.5] table[skip coords between index=67] num_exp/out/conv_kite_DG_ks4_kt4_both_nref7_gamma0.05_sm2.dat;
+[thick,mark options=scale=1.5] table[skip coords between index=56] num_exp/out/conv_kite_DG_ks5_kt5_both_nref6_gamma0.05_sm2.dat;
+[thick,mark options=scale=1.5] table[skip coords between index=45] num_exp/out/conv_kite_DG_ks6_kt6_both_nref5_gamma0.05_sm2.dat;
[gray, dashed, domain=0:7] (1/2^(x-0.2)))^2;
[gray, dashed, domain=0:7] (1/2^(x+0.5)))^3;
[gray, dashed, domain=0:7] (1/2^(x+1)))^4;
[gray, dashed, domain=0:5] (1/2^(x+1)))^5;
[gray, dashed, domain=0:4] (1/2^(x+1)))^6;
[gray, dashed, domain=0:3] (1/2^(x+1)))^7;

[scale=0.75] at (4.5,5.25) $O(h^{k+1})= O(\Delta t^{k+1})$;
[scale=0.75, gray, dash=on 2.25pt off 2.25pt phase 0pt, line width=0.4*0.75pt] (6.25/0.75,5.25/0.75) – (6.8/0.75,5.25/0.75);
[ xlabel=$i$, ylabel=$\| u_h(T) - u (T)\|_{L^2(\Omega(T))}$, legend entries =$\gamma_J = 5 \cdot 10^4$,$\gamma_J = 5$, $\gamma_J = 5 \cdot 10^{-2}$,$\gamma_J = 5 \cdot 10^{-4}$, legend style=anchor=north,legend columns=1, draw=none, fill=none,legend pos =south west, x label style=at=(axis description cs:0.6,0.15),anchor=center]
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_DG_ks4_kt4_both_nref5_gamma50000.0_sm2.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_DG_ks4_kt4_both_nref5_gamma5.0_sm2.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_DG_ks4_kt4_both_nref5_gamma0.05_sm2.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_DG_ks4_kt4_both_nref5_gamma0.0005_sm2.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_DG_ks4_kt4_both_nref5_gamma5e-10_sm2.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_DG_ks4_kt4_both_nref3_gamma0.0_sm2.dat;
[gray, dashed, domain=0:4] (1/2^(x+1.3)))^5;

at (rel axis cs: 0.75,0.8) [r]
$\gamma_J = 5 \cdot 10^{-10}$ <ref>
$\gamma_J = 0$ <ref>;

[scale=0.75] at (4.5,5.25) $O(h^5)=O(\Delta t^5)$;
[scale=0.75, gray, dash=on 2.25pt off 2.25pt phase 0pt, line width=0.4*0.75pt] (6./0.75,5.25/0.75) – (6.55/0.75,5.25/0.75);

Left: Convergence of the DG method for the kite case. Right: Influence of the stabilisation parameter $\gamma_J$ on the convergence behaviour of the DG method for $\ksol=\klset = 4$.

[ xlabel=$i$, ylabel=$\| u_h(T) - u (T)\|_{L^2(\Omega(T))}$, legend entries = $k=1$, $k=2$, $k=3$, $k=4$, $k=5$, $k=6$, legend style=anchor=north,legend columns=1, draw=none, legend cell align=left, fill=none,legend pos =south west,
x label style=at=(axis description cs:0.5,0.15),anchor=center,
ymax =2.5
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_CG_ks1_kt1_both_nref8_gamma0.05_sm2.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_CG_ks2_kt2_both_nref8_gamma0.05_sm2.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_CG_ks3_kt3_both_nref8_gamma0.05_sm2.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_CG_ks4_kt4_both_nref7_gamma0.05_sm2.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_CG_ks5_kt5_both_nref6_gamma0.05_sm2.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_CG_ks6_kt6_both_nref5_gamma0.05_sm2.dat;
[gray, dashed, domain=0:7] (1/2^(x+0.2)))^2;
[gray, dashed, domain=0:7] (1/2^(x+0.5)))^3;
[gray, dashed, domain=0:7] (1/2^(x+0.6)))^4;
[gray, dashed, domain=0:6] (1/2^(x+1)))^5;
[gray, dashed, domain=0:5] (1/2^(x+1)))^6;
[gray, dashed, domain=0:4] (1/2^(x+1)))^7;

[scale=0.75] at (4.5,5.25) $O(h^{k+1})= O(\Delta t^{k+1})$;
[scale=0.75, gray, dash=on 2.25pt off 2.25pt phase 0pt, line width=0.4*0.75pt] (6.25/0.75,5.25/0.75) – (6.8/0.75,5.25/0.75);

[ xlabel=$i$, ylabel=$\| u_h(T) - u (T)\|_{L^2(\Omega(T))}$, legend entries = $k=3$, legend style=anchor=north,legend columns=1, draw=none, legend pos =south west, x label style=at=(axis description cs:0.5,0.15),anchor=east]
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_GCC_ks3_kt3_both_nref8_gamma0.05_sm2.dat;
[gray, dashdotted, domain=0:7] (1/2^(x+1.)))^4;

[scale=0.75] at (4.5,5.25) $O(h^4)=O(\Delta t^4)$;
[scale=0.75, gray, dash=on 2.25pt off 2.25pt phase 0pt, line width=0.4*0.75pt] (6./0.75,5.25/0.75) – (6.55/0.75,5.25/0.75);

Convergence of the higher order CG (left) and GCC (right) methods for the kite case.

We measure the numerical error in terms of the following norms:

\begin{align*}
\| u \|_{L^2(\Omega(T))}^2 := \int_{\Omega(T)} u(\mathbf{x}, T)^2 \ \mathrm{d}\mathbf{x}, \quad
\| u \|_{L^2( L^2(\Omega(t)), 0, T)}^2 := \int_0^T \int_{\Omega(t)} u(\mathbf{x}, t)^2 \ \mathrm{d}\mathbf{x} \ \mathrm{d}t.
\end{align*}

We find both errors to decrease as follows, where $k = \khsol = \ktsol = \khlset = \ktlset$, $h \sim \Delta t$,

\begin{equation}
\| u - u_h \|_{L^2(\Omega(T))} + \| u - u_h \|_{L^2( L^2(\Omega(t)), 0, T)} = \mathcal{O}(h^{k+1}) = \mathcal{O}(\Delta t^{k+1}). \label{eq:generalerrorbound} %\mathcal{O}(h^{k_s+1} + \Delta t^{k_t + 1}). \label{eq:generalerrorbound}
\end{equation}

This confirms the expected higher order convergence property of the suggested methods.

§.§.§ Choice of the stabilisation parameter

The stabilisation parameter of the ghost penalty has to be chosen large enough to fulfill sufficient conditions for stability of the method and to obtain reasonable bounds on the condition numbers of arising linear systems [10, 35, 34].

In <ref> (right) we illustrate the dependence of the numerical error on $\gamma_J$ for the problem considered in this subsection using $i=i_s=i_t$ and $\ksol = \klset = 4$ for the DG method.

Overall, $\gamma_J$ influences the absolute values of the errors, whilst the asymptotic behaviour remains intact even if the parameter is chosen relatively large. A particular limit is that of stabilisation constant $\gamma_J \to 0$. There, we observe that the method is stable in general. However, condition numbers may blow up, which impedes the practical solution of the linear system at some refinement level/ polynomial degree. (In this case, $i=2$ is the last level solved properly.) We explain this behavior as follows: For mere stability, i.e. assuming exact arithmetics, ghost penalty stabilization is not even necessary. Note that in [26] even robust, i.e. cut-position independent, a priori error bounds have been derived without ghost penalty stabilization in a comparable setting (involving the assumption of exact integration). However, taking computer arithmetics, i.e. round-off errors and their impact, into account the possible ill-conditioning of linear systems is practically hardly acceptable. A more detailed investigation of these observations (potentially also in combination with different time integration strategies) in terms of a full numerical analysis would be an interesting task for future research.

In conclusion, we may pick $\gamma_J$ small as long as a sufficiently accurate (and efficient) solution of linear systems is possible. We choose $\gamma_J=0.05$ for the remaining studies.

§.§.§ Superconvergence

The investigations in the first part of this subsection confirmed the numerical error to scale at least with $\mathcal{O}(\Delta t^{\ktsol + 1})$ in the mentioned norms. A well-known property of DG type discretisations is that of superconvergence, i.e. the convergence in time with a higher order (at specific time instances, e.g. nodes; cf. also [4] for a superconvergence result for a CG method of the heat equation.). We investigate numerically whether our method shows such a behaviour in a two-stage procedure: In a first experiment, we perform only time-refinements for a fine mesh and a high order discretisation in space. This experiment allows to estimate the order of (super)-convergence in time, although asymptotically this order will be hard to reach exactly because a saturation from the spatial error will occur. Hence, we validate the first estimate for the superconvergence in a second step by a space-time-refinement study involving the suiting spatial discretisation order. The results for the first study, involving the choices of order $(\ktsol,\khsol)=(1,3), (2,5)$, $\klset = \khsol$ and a fixed mesh in space ($i_s = \mathrm{const}$) are shown in <ref>.

[ xlabel=$i_t$, ylabel=$\| u_h(T) - u (T)\|_{L^2(\Omega(T))}$, legend entries =DG $\ktsol=1$, DG $\ktsol=2$, CG $\ktsol=1$, CG $\ktsol=2$, legend style=,anchor=north,legend columns=1, draw=none,legend pos =south west, x label style=at=(axis description cs:0.5,0.15),anchor=center]
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_DG_ks3_kt1_time_nref6_gamma0.05_sm2_ktls3_ro-1.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_DG_ks5_kt2_time_nref5_gamma0.05_sm2_ktls5_ro-1.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_CG_ks3_kt1_time_nref6_gamma0.05_sm2_ktls3_ro-1.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_CG_ks5_kt2_time_nref5_gamma0.05_sm2_ktls5_ro-1.dat;
[very thick, dashed, domain=2:4] (1/2^(x+3.2)))^3.8;
[very thick, dashdotted, domain=3:5] (1/2^(x+2.5)))^2.8;
[gray, dashed, domain=0:4.25] (1/2^(x+3.7)))^3;
[gray, dashed, domain=0:4] (1/2^(x+3.3)))^4;

[gray, dashed, domain=0:5] (1/2^(x+2.5)))^2;

[gray] at (rel axis cs: 0.85,0.67) $O(\Delta t^2)$;
[gray] at (rel axis cs: 0.9,0.2) $O(\Delta t^3)$;
[gray] at (rel axis cs: 0.65,0.07) $O(\Delta t^4)$;

[black] at (rel axis cs: 0.75,0.85) [r]  <ref> $O(\Delta t^{2.8})$
<ref> $O(\Delta t^{3.8})$;

[ xlabel=$i_t$, ylabel=$\| u_h - u\|_{L^2( L^2(\Omega(t)), 0, T) }$, legend entries =DG $\ktsol=1$, DG $\ktsol=2$, CG $\ktsol=1$, CG $\ktsol=2$, legend style=anchor=north, legend columns=1, draw=none,legend pos =south west, x label style=at=(axis description cs:0.5,0.15),anchor=center]
+[thick,mark options=scale=1.5] table [x index = 0, y index =2] num_exp/out/conv_kite_DG_ks3_kt1_time_nref6_gamma0.05_sm2_ktls3_ro-1.dat;
+[thick,mark options=scale=1.5] table [x index = 0, y index =2] num_exp/out/conv_kite_DG_ks5_kt2_time_nref5_gamma0.05_sm2_ktls5_ro-1.dat;
+[thick,mark options=scale=1.5] table [x index = 0, y index =2] num_exp/out/conv_kite_CG_ks3_kt1_time_nref6_gamma0.05_sm2_ktls3_ro-1.dat;
+[thick,mark options=scale=1.5] table [x index = 0, y index =2] num_exp/out/conv_kite_CG_ks5_kt2_time_nref5_gamma0.05_sm2_ktls5_ro-1.dat;
[gray, dashed, domain=0:5] (1/2^(x+3.0)))^2;
[gray, dashed, domain=0:4] (1/2^(x+3.2)))^3;

[gray] at (rel axis cs: 0.85,0.55) $O(\Delta t^2)$;
[gray] at (rel axis cs: 0.88,0.22) $O(\Delta t^3)$;

Superconvergence study of the DG and CG methods for $\ktsol=1,2$ where $\khsol = 2 \ktsol +1$. Refinement in time for a fine mesh. On the left-hand side, the error is measured in the $L^2(\Omega(T))$-norm, whilst on the right-hand side the same results are investigated in the $L^2( L^2(\Omega(t)), 0, T)$-norm.

Whilst in the space-time norm $L^2( L^2(\Omega(t)), 0, T)$ the estimate of <ref> is confirmed, for $L^2(\Omega(T))$ we observe for the DG method and the CG method with $\ktsol = 2$, $\| u - u_h \|_{L^2(\Omega(T))} \approx \mathcal{O}(\Delta t^{\ktsol + 1.8})$, which in light of the aforementioned saturation effect suggests

\begin{equation}
\| u - u_h \|_{L^2(\Omega(T))} = \mathcal{O}(\Delta t^{\ktsol + 2}), \textnormal{for DG, CG with } \ktsol = 2. \label{eq:superconverrorbound}
\end{equation}

For CG with $\ktsol = 1$, the estimate <ref> is also observed to be sharp, there is no superconvergence. In a second step, we confirm <ref> to hold with a space-time refinement study involving $(\ktsol,\khsol)=(1,2), (2,3)$, $\klset = \khsol$, see <ref>.

[ xlabel=$i$, ylabel=$\| u_h(T) - u (T)\|_{L^2(\Omega(T))}$, legend entries =DG $\ktsol=1$, DG $\ktsol=2$, CG $\ktsol=2$, legend style=,anchor=north,legend columns=1, draw=none,legend pos =south west, x label style=at=(axis description cs:0.5,0.15),anchor=center]
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_DG_ks2_kt1_both_nref7_gamma0.05_sm2_ktls2.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_DG_ks3_kt2_both_nref6_gamma0.05_sm2_ktls3.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_kite_CG_ks3_kt2_both_nref6_gamma0.05_sm2_ktls3.dat;

[gray, dashed, domain=0:6] (1/2^(x+1.3)))^3;
[gray, dashed, domain=0:5] (1/2^(x+1.3)))^4;

[gray] at (rel axis cs: 0.9,0.35) $O(\Delta t^3)$;
[gray] at (rel axis cs: 0.65,0.07) $O(\Delta t^4)$;

Superconvergence study of the DG and CG methods for $\ktsol=1,2$ where $\khsol = \ktsol +1$. Refinement in space and time. The error is measured in the $L^2(\Omega(T))$-norm.

We mention that theoretical results of findings similar to <ref> exist for static domains, c.f. e.g. [14, 15], and leave a further theoretical study of these phenomena in this particular setting for future research.

§.§ Moving $n$-sphere

The second example geometry is a moving circle in 2D/ sphere in 3D/ interval in 1D. These cases involve no change in deformation, but exemplify a moving domain setting. In dimension $d$, they are described by

\begin{align}
\rho(t) &= \frac{1}{\pi} \sin(2 \pi t), \quad \Omega = [-1,1] \times [-0.6,0.6]^{d-1}, \quad \mathbf{w} = (\frac{\partial \rho}{\partial t}, 0, \dots, 0), \\
r(\mathbf{x}, t) &= [ (\mathbf{x}_1 - \rho(t))^2 + {\sum}_{i=2}^d \mathbf{x}_i^2]^{1/2} \quad \phi = r - 0.5.
\end{align}

The manufactured solution $u$ remains $u = \cos (\pi \frac{r}{r_0}) \cdot \sin (\pi t)$, leading to slightly different expressions for $f$ in each dimension. Moreover, $h= 0.5^{(i_s + 1)}, \Delta t = 2^{-i_t - 2}, \gamma_J = 0.05, \epsilon_f=1.1$.

§.§.§ On 1D and 3D

In this part, we want to show that the convergence results <ref> also hold for spatially one- and three-dimensional problems. To this end, we apply the methods of consideration to the mentioned problems and study the decay of the numerical error for $\ksol=\klset$, $i_s=i_t$. In <ref>, the DG method is applied in three spatial dimensions and the CG method is considered in 1D.

[ xlabel=$i$, ylabel=$\| u_h(T) - u (T)\|_{L^2(\Omega(T))}$, legend entries = $k=1$, $k=2$, $k=3$, $k=4$, legend style=anchor=north,legend columns=1, draw=none, fill=none,legend pos =south west, x label style=at=(axis description cs:0.5,0.15),anchor=center, ymax=1]
+[thick,mark options=scale=1.5] table num_exp/out/conv_moving_sphere_DG_ks1_kt1_both_nref6_gamma0.05_sm2.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_moving_sphere_DG_ks2_kt2_both_nref5_gamma0.05_sm2.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_moving_sphere_DG_ks3_kt3_both_nref4_gamma0.05_sm2.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_moving_sphere_DG_ks4_kt4_both_nref3_gamma0.05_sm2.dat;
[gray, dashed, domain=0:5] (1/2^(x-0.2)))^2;
[gray, dashed, domain=0:4] (1/2^(x+0.5)))^3;
[gray, dashed, domain=0:3] (1/2^(x+0.75)))^4;
[gray, dashed, domain=0:2] (1/2^(x+1.5)))^5;

[scale=0.75] at (4.5,5.25) $O(h^{k+1})= O(\Delta t^{k+1})$;
[scale=0.75, gray, dash=on 2.25pt off 2.25pt phase 0pt, line width=0.4*0.75pt] (6.25/0.75,5.25/0.75) – (6.8/0.75,5.25/0.75);

[ xlabel=$i$, ylabel=$\| u_h(T) - u (T)\|_{L^2(\Omega(T))}$, legend entries = $k=1$, $k=2$, $k=3$, $k=4$, legend style=legend columns=1, draw=none, legend pos =south west, x label style=at=(axis description cs:0.8,0.15),anchor=center, ymin=5e-15]
+[thick,mark options=scale=1.5] table num_exp/out/conv_moving_int_CG_ks1_kt1_both_nref12_gamma0.05_sm1.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_moving_int_CG_ks2_kt2_both_nref11_gamma0.05_sm1.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_moving_int_CG_ks3_kt3_both_nref10_gamma0.05_sm1.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_moving_int_CG_ks4_kt4_both_nref9_gamma0.05_sm1.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_moving_int_CG_ks5_kt5_both_nref8_gamma0.05_sm1.dat;
+[thick,mark options=scale=1.5] table[skip coords between index=67] num_exp/out/conv_moving_int_CG_ks6_kt6_both_nref7_gamma0.05_sm1.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_moving_int_CG_ks7_kt7_both_nref6_gamma0.05_sm1.dat;
+[thick,mark options=scale=1.5] table num_exp/out/conv_moving_int_CG_ks8_kt8_both_nref5_gamma0.05_sm1.dat;
[gray, dashed, domain=0:11] (1/2^(x-0.2)))^2;
[gray, dashed, domain=0:10] (1/2^(x+0.5)))^3;
[gray, dashed, domain=0:9] (1/2^(x+0.6)))^4;
[gray, dashed, domain=0:8] (1/2^(x+0.7)))^5;
[gray, dashed, domain=0:7] (1/2^(x+0.25)))^6;
[gray, dashed, domain=0:5] (1/2^(x+0.25)))^7;
[gray, dashed, domain=0:5] (1/2^(x+0.25)))^8;
[gray, dashed, domain=0:4] (1/2^(x+0.25)))^9;

at (rel axis cs: 0.65,0.78) [r]
<ref> $O(h^{k+1})= O(\Delta t^{k+1})$
<ref> $k=5$

<ref> $k=6$

<ref> $k=7$
<ref> $k=8$;

Convergence of DG method for 3+1D (left) and the CG method for 1+1D (right).

The respective choice of refinement levels and tested discretisation orders reflects the different computational demands. The results in the $L^2( L^2(\Omega(t)), 0, T)$-norm are equivalent and shown in
the supplement, as well as results for the DG method in 1D and the CG method in 3D. In all cases, the convergence orders of <ref> are confirmed.

§.§.§ Non-zero entries of the matrices

For an evaluation of the computational costs, we compare the number of non-zero entries (s) of the discrete problem matrix of the respective methods. We fix the setting of the moving circle, $\ksol = \klset $, $i_s=2$. On each time-slice, after the system matrix is assembled, s are counted and the minimal and maximal value of this number over all slices is considered.

For the CG method, in comparison to the DG method, s will depend on the factor $\epsilon_f$, as we evaluate in a first study. To this end, we choose $i_t=2$ and evaluate s for the CG method with $\epsilon_f \in \{ 1, 1.1, 1.5, \dots, 100\}$, and for comparison for the DG method, for the polynomial degrees $k=1,2,3$. The results are displayed in <ref>.
[ xlabel=$\epsilon_f$, ylabel=nzes, legend entries =CG $k = 1$, CG$\square$ $k = 1$, DG $k=1$, CG $k=2$, CG$\square$ $k=2$, DG $k=2$, CG $k=3$, CG$\square$ $k=3$, DG $k=3$, CG $k=4$, CGtg $k=4$, DG $k=4$, legend style=at=(1.4,0.5),anchor=center,legend columns=1, draw=none, legend cell align=left,x label style=at=(axis description cs:0.5,0.15),anchor=east,y label style=at=(axis description cs:0.05,0.5)]
+[mark options=scale=1.5,thick] table [x index = 0, y index =3] num_exp/out/nze_moving_circle_CG_ks1_kt1_ktls1tp_gp0.dat;
table [x index = 0, y index =3] num_exp/out/nze_moving_circle_CG_ks1_kt1_ktls1tp_gp1.dat;
[no markers, thick] table [x index = 0, y index =3] num_exp/out/nze_moving_circle_DG_ks1_kt1_ktls1.dat;
+[mark options=scale=1.5,thick] table [x index = 0, y index =3] num_exp/out/nze_moving_circle_CG_ks2_kt2_ktls2tp_gp0.dat;
+[mark options=scale=1.5,thick] table [x index = 0, y index =3] num_exp/out/nze_moving_circle_CG_ks2_kt2_ktls2tp_gp1.dat;
[no markers, dashed, thick] table [x index = 0, y index =3] num_exp/out/nze_moving_circle_DG_ks2_kt2_ktls2.dat;
+[mark options=scale=1.5,thick] table [x index = 0, y index =3] num_exp/out/nze_moving_circle_CG_ks3_kt3_ktls3tp_gp0.dat;
+[mark options=scale=1.5,thick] table [x index = 0, y index =3] num_exp/out/nze_moving_circle_CG_ks3_kt3_ktls3tp_gp1.dat;
[no markers, dotted, thick] table [x index = 0, y index =3] num_exp/out/nze_moving_circle_DG_ks3_kt3_ktls3.dat;

Non-zero entries for DG and CG in comparison for different parameter choices for $\epsilon_f$. Each value represents the maximal number of non-zero entries over all time steps.

First, we observe that over all polynomial degrees the CG method leads to less non-zero entries as the DG method, whilst the relative difference is specifically strong for small orders $k$. This aligns with the elimination of one  of the time finite element space. Second, this benefit of the CG method is reduced by a high choice of the factor $\epsilon_f$, as a larger domain of extension will induce more entries in the matrix. Third, the original CG method of <Ref> performs better than the variant mentioned in <ref>. Hence, to exploit the computational benefits of the CG method fully, we suggest to use the non-tensor-product variant of the ghost penalty on a region as small as possible.

In a second study, we fix $\epsilon_f=1.1$ and compare the non-zero entries for all three methods for different time step sizes $\Delta t$ and different polynomial orders. The numbers are given in <ref>.

3c$i_t = 1$
3c$i_t = 3$
3c$i_t = 5$

2cs    1c
2cs    1c
2cs    1c

um_exp/nze2_table_k1.tex um_exp/nze2_table_k3.tex um_exp/nze2_table_k5.tex

Range for non-zero entries (s) of system matrix (depends on time step) and absolute $L^2(T)$-error () in comparison for DG, CG and GCC. In this whole table, $i_s=2$.

Summarizing the results within each cell, the number of non-zero entries decreases from DG to CG and in turn to GCC. In particular for discretisation parameter $i_t=3,5$, this holds within a regime of similar absolute numerical errors. Within each row of the table, we note that smaller time steps decrease the numbers of non-zero entries for all methods, as the interface movement per time step decreases also. This effect is stronger for CG and GCC because of the extended regions of the ghost penalty stabilisation. For increasing polynomial order, we again observe a relatively smaller benefit of the CG and GCC method as compared to DG. We conclude that those methods are specifically interesting for orders such as $k=1,2,3$ from a computational point of view.

§.§.§ Different time integration strategies

In this part we investigate numerical differences between the topology-insensitive and the topology-preserving method of numerical quadrature in time suggested in <Ref>. For this purpose, we implement both approaches and compare the respective numerical errors. By and large, we observe competitive numerical results over most of the problem settings introduced so far. This is in agreement with the good results obtained in the literature, cf. [38, 18]. However, we observed differences in the following setup with almost no ghost-penalty stabilisation: Consider the DG method applied to the moving interval in 1D with a radius $R=0.505$, constant velocity $\rho = 0.5 t$, manufactured right-hand side for $u= (x-\rho +R)^2 \cdot (x - \rho -R)^2$, and polynomial order $\ksol = \klset = 4$. For this example, we fix $i_t = 0$ and refine in space. From $i_s=2$ on, the discrete interface is captured exactly so that the method can approximate the solution exactly as a polynomial of order 4.

[ xlabel=$i_s$, ylabel=$\| u_h(T) - u (T)\|_{L^2(\Omega(T))}$,
title style=at=(0.4,0.96),draw=black, fill=white,anchor=north,
x label style=at=(0.97,0.18),anchor=west, below=5mm,
title=$\gamma_J = 0.05$, legend entries =top.-pres., top.-insen., top.-insen. double order, top.-insen. 10 substeps, legend style=legend columns=2, draw=none, legend to name=named, ymax=10, ymin=1e-17]
+[thick, mark options=scale=1.5] table num_exp/out/conv_moving_int_sp_DG_ks4_kt4_space_nref7_gamma0.05_sm1_tint_proper_ro-7.dat;
+[thick, mark options=scale=1.5] table num_exp/out/conv_moving_int_sp_DG_ks4_kt4_space_nref7_gamma0.05_sm1_tint_naive_ro-7.dat;
+[thick, mark options=scale=1.5] table num_exp/out/conv_moving_int_sp_DG_ks4_kt4_space_nref7_gamma0.05_sm1_tint_naivedo_ro-7.dat;
+[thick, mark options=scale=1.5] table num_exp/out/conv_moving_int_sp_DG_ks4_kt4_space_nref7_gamma0.05_sm1_tint_naiveS10_ro-7.dat;

[ xlabel=$i_s$,
title style=at=(0.4,0.96),draw=black, fill=white,anchor=north,
x label style=at=(0.97,0.18),anchor=west, below=5mm,
title=$\gamma_J = 5 \cdot 10^{-10}$, ymax=10, ymin=1e-17]
+[thick, mark options=scale=1.5] table num_exp/out/conv_moving_int_sp_DG_ks4_kt4_space_nref7_gamma5e-10_sm1_tint_proper_ro-7.dat;
+[thick, mark options=scale=1.5] table num_exp/out/conv_moving_int_sp_DG_ks4_kt4_space_nref7_gamma5e-10_sm1_tint_naive_ro-7.dat;
+[thick, mark options=scale=1.5] table num_exp/out/conv_moving_int_sp_DG_ks4_kt4_space_nref7_gamma5e-10_sm1_tint_naivedo_ro-7.dat;
+[thick, mark options=scale=1.5] table num_exp/out/conv_moving_int_sp_DG_ks4_kt4_space_nref7_gamma5e-10_sm1_tint_naiveS10_ro-7.dat;

[ xlabel=$i_s$,
title style=at=(0.4,0.96),draw=black, fill=white,anchor=north,
x label style=at=(0.97,0.18),anchor=west, below=5mm,
title=$\gamma_J = 5 \cdot 10^{-20}$, ymax=10, ymin=1e-17]
+[thick, mark options=scale=1.5] table num_exp/out/conv_moving_int_sp_DG_ks4_kt4_space_nref7_gamma5e-20_sm1_tint_proper_ro-7.dat;
+[thick, mark options=scale=1.5] table num_exp/out/conv_moving_int_sp_DG_ks4_kt4_space_nref7_gamma5e-20_sm1_tint_naive_ro-7.dat;
+[thick, mark options=scale=1.5] table num_exp/out/conv_moving_int_sp_DG_ks4_kt4_space_nref7_gamma5e-20_sm1_tint_naivedo_ro-7.dat;
+[thick, mark options=scale=1.5] table num_exp/out/conv_moving_int_sp_DG_ks4_kt4_space_nref7_gamma5e-20_sm1_tint_naiveS10_ro-7.dat;

Topology-insensitive and topology-preserving time integration in comparison.

The results of this study are shown in <ref>. For the method with the topology-preserving time integration, we note that there appear only errors below a reasonable round-off error tolerance such as $10^{-12}$. Meanwhile, the method of topology-insensitive time integration leads to significant errors when very small stabilisation constants and mesh sizes are used. This is caused by cut configurations where an element is included in the active mesh because of a small cut, but the corresponding matrix entries vanish because of the numerical integration in time. In the presence of a moderate ghost penalty stabilisation, the  are still controlled by this stabilisation, otherwise stability issues appear. For an interpretation of this behaviour in terms of numerical analysis, we refer the reader to the comment made about $\gamma_J \to 0$ in <Ref>. Moreover, from a numerical analysis point of view we expect that the topology-preserving time integration is important for the geometrical consistency. In the numerical experiments however, we did not observe corresponding issues with the topology-insensitive time integration, at least not in combination with sufficient ghost penalty stabilization.
We are not sure if there is an underlying mechanism that can explain this behaviour or if the test case is still too mild and leave further investigations of this for future research.

In addition we compared the computational effort for both time integration strategies for the two-dimensional configuration of <ref>. As differences only occur on a small subset of elements in the computation domain, the overall computation time is only marginally affected by the choice of the integration strategy.

We conclude that the topology-preserving time integration represents the discrete geometry properly, and hence is our method of choice, although the topology-insensitive time integration seems to perform equally well with sufficient stabilisation. Finally, we note that by doubling the order or introducing 10 evenly spaced substeps, the errors for the topology-insensitive time integration approach with degenerate stabilisation can be improved.[Here, doubling the order amounts to taking a plain Gaussian integration rule of order $4(k_t +1) =20$ instead of the plain rule of order $2 (k_t + 1) =10$ (c.f. Section 5). Moreover, introducing 10 evenly spaced substeps amounts to subdividing the interval $I_n$ into 10 subintervals of the same length and applying the Gaussian rule of order $2 (k_t +1) = 10$ to each of those.]

§.§ Colliding circles

Thirdly, we want to investigate the stability of the method in a topologically challenging example taken from [25]: Two circles collide and separate again afterwards. This setting can be described by the following level set function:

\begin{equation}
\phi(\mathbf{x},t) = \mathrm{min}(\| \mathbf{x} - s_1(t) \|, \| \mathbf{x} - s_2(t) \|) - R, \nonumber
s_1(t) = (0, t - 3/4), s_2(t) = (0, 3/4 - t).
\end{equation}

The end time is set to $T = 3/2$, such that $\phi(\dots, 0) = \phi(\dots, T)$, the diffusion constant $\alpha = 0.1$, $R=0.5$, $\Omega=[-0.6,0.6] \times [-1.35,1.35]$, and the associated convection velocity is

\begin{equation}
\mathbf{w}(x,y,t) = \begin{cases}
(0,-1)^T \quad \textnormal{if } (y > 0 \textnormal{ and } t \leq T/2) \textnormal{ or } (y < 0 \textnormal{ and } t > T/2) \\
(0,1)^T \quad \textnormal{if } (y \leq 0 \textnormal{ and } t \leq T/2) \textnormal{ or } (y > 0 \textnormal{ and } t > T/2).
\end{cases}
\end{equation}

Here, $u_0$ is defined to be $u_0(x,y) := \mathrm{sign}(y)$, i.e.  1 in the upper circle and -1 in the lower one.

First, we investigate the numerical method with $\ksol = \klset = 4$, $i_s = 2$, $\Delta t= T/16$ to generate a reference solution, which is shown on the left-hand side of <ref>. As one would expect physically, first the circles are transported with the convection field. After they merged, the diffusion successively evens out the concentration difference. By and large, the simulation does not show any numerical instability or large error. As a small detail, we notice that some numerical diffusion will appear shortly before the contact of the circles, when the numerically extended regions merge. This was also observed in [25] for the time stepping method presented there. But for fine meshes or small stabilisation constants, this effect can be kept small. We conclude that the method performs well, despite a challenging topology change, without the need of any specific adaptation.

Discrete solutions for the colliding circles test case. Left: Reference configuration. Center: Single-time-step run with $\ksol=4$. Right: Single-time-step run with $\ksol=2$.

Second, we consider a more extreme version by setting $\Delta t = T$, so that the whole time interval will be solved in one time step. We show the numerical solutions for the parameter choices $\ksol = \klset = 2, 4$ in <ref>. As the overall topology of the problem is still resolved properly, the method appears capable of handling topologically demanding configurations properly.

§ CONCLUSION

In this final section, we want to summarise the findings of this work and highlight some questions that remain for future research.

We introduced unfitted space-time methods for a
problem posed on a moving domain. These are based on discontinuous Galerkin, continuous Galerkin, and higher order continuous Galerkin-Collocation discretisations in time. We further introduced techniques to achieve a high order geometry handling of the space-time domains based on isoparametric mappings of the background mesh.

In extensive numerical experiments we discussed and compared several important properties of the proposed methods. Especially, the expectations of robustness and higher order convergence in space and time have been confirmed. Even in topologically challenging configurations, our methods perform reliably.

These results motivate the following directions for future investigations:

First, the methods have been investigated from a computational point of view in this paper. This suggests to investigate whether the numerically established convergence estimates can be proven from the point of view of numerical analysis. As far as the DG method is concerned, the master's theses [34] and [19] present such error estimates, although in the former an exact handling of the geometries is assumed, which is weakened in the latter in regard to the spatial approximation.
This leaves a numerical analysis relative to the spatially and temporally discrete geometry an open task. Moreover, for the CG and GCC method no such analysis has been performed to the author's best knowledge.

Next, the numerical experiments w.r.t. the different numerical integration approaches suggest that the topology-insensitive time integration is also a reasonable and accurate approach
and the topology preserving time integration seems to be superfluous in the presence of ghost penalty stabilisation. This behavior is not fully understood and leaves room for future research.

Finally, let us mention that although we constructed a space-time mesh deformation that can be discontinuous between time slabs, it may be possible to construct a continuous mesh deformation that preserves the accuracy of the scheme. This would allow to remove the transfer operation $\Pi_u$ and hence simplify the method and potentially its analysis.

[1]
M. Anselmann and M. Bause, CutFEM and ghost stabilization
techniques for higher order space-time discretizations of the
Navier-Stokes equations,  (2021),

[2]
M. Anselmann and M. Bause, Higher order Galerkin–collocation
time discretization with Nitsche’s method for the Navier–Stokes
equations, Math. Comput. Simul., 189 (2021), pp. 141–162,

[3]
Anselmann, Mathias, Bause, Markus, Becher, Simon, and Matthies,
Gunar, Galerkin-collocation approximation in time for the wave
equation and its post-processing, ESAIM: M2AN, 54 (2020), pp. 2099–2123,

[4]
A. K. Aziz and P. Monk, Continuous Finite Elements in Space
and Time for the Heat Equation, Math. Comput., 52 (1989),
pp. 255–274, <https://doi.org/10.2307/2008467>.

[5]
A. Barrett, A. L. Fogelson, and B. E. Griffith, A hybrid
semi-Lagrangian cut cell method for advection-diffusion problems with
Robin boundary conditions in moving domains, J. Comput. Phys., 449 (2022),

[6]
S. Becher and G. Matthies, Variational time discretizations of
higher order and higher regularity, BIT Numer. Math., 61 (2021),
pp. 721–755, <https://doi.org/10.1007/s10543-021-00851-6>.

[7]
S. Becher, G. Matthies, and D. Wenzel, Advanced Computing in
Industrial Mathematics: BGSIAM 2017, Springer International Publishing,
Cham, 2019, ch. Variational Methods for Stable Time Discretization of
First-Order Differential Equations, pp. 63–75,

[8]
A. Brenner, E. Bänsch, and M. Bause, A priori error analysis for
finite element approximations of the Stokes problem on dynamic meshes, IMA
Journal of Numerical Analysis, 34 (2014), pp. 123–146,

[9]
E. Burman, Ghost penalty, C.R. Math. Acad. Sci. Paris, 348 (2010),
pp. 1217–1220, <https://doi.org/10.1016/j.crma.2010.10.006>.

[10]
E. Burman, S. Claus, P. Hansbo, M. G. Larson, and A. Massing, Cutfem: discretizing geometry and partial differential equations, Int. J.
Numer. Methods Eng., 104 (2015), pp. 472–501,

[11]
E. Burman, S. Frei, and A. Massing, Eulerian time-stepping schemes
for the non-stationary Stokes equations on time-dependent domains, Numer.
Math.,  (2022), <https://doi.org/10.1007/s00211-021-01264-x>.

[12]
E. Burman and P. Hansbo, Fictitious domain finite element methods
using cut elements: I. a stabilized Lagrange multiplier method, Comput.
Methods Appl. Mech Eng., 199 (2010), pp. 2680–2686,

[13]
E. Burman and P. Hansbo, Fictitious domain finite element methods
using cut elements: II. a stabilized Nitsche method, Appl. Numer. Math.,
62 (2012), pp. 328–341, <https://doi.org/10.1016/j.apnum.2011.01.008>.

[14]
K. Eriksson, C. Johnson, and V. Thomée, Time discretization of
parabolic problems by the discontinuous Galerkin method, ESAIM:
Mathematical Modelling and Numerical Analysis, 19 (1985), pp. 611–643,

[15]
A. Ern and F. Schieweck, Discontinuous Galerkin method in time
combined with a stabilized finite element method in space for linear
first-order PDEs, Mathematics of Computation, 85 (2016), pp. 2099–2129,

[16]
T.-P. Fries and T. Belytschko, The extended/generalized finite
element method: An overview of the method and its applications, Int. J.
Numer. Methods Eng., 84 (2010), pp. 253–304,

[17]
S. Groß and A. Reusken, Numerical Methods for Two-phase
Incompressible Flows, Springer, Berlin, 2011,

[18]
P. Hansbo, M. G. Larson, and S. Zahedi, A cut finite element method
for coupled bulk-surface problems on time-dependent domains, Comput. Methods
Appl. Mech. Engrg., 307 (2016), pp. 96–116,

[19]
F. Heimann, On Discontinuous- and Continuous-In-Time
Unfitted Space-Time Methods for PDEs on Moving Domains,
Master's thesis, University of Göttingen,  (2020),

[20]
F. Heimann and C. Lehrenfeld, Numerical Integration on
Hyperrectangles in Isoparametric Unfitted Finite Elements, in
European Conference on Numerical Mathematics and Advanced Applications,
Springer, 2019, pp. 193–202,

[21]
C. Hirt, A. Amsden, and J. Cook, An arbitrary
Lagrangian-Eulerian computing method for all flow speeds, J. Comput.
Phys., 14 (1974), pp. 227–253,

[22]
C. Lehrenfeld, The Nitsche XFEM-DG space-time method and its
implementation in three space dimensions, SIAM J. Sci. Comp., 37 (2015),
pp. A245–A270, <https://doi.org/10.1137/130943534>.

[23]
C. Lehrenfeld, High order unfitted finite element methods on level
set domains using isoparametric mappings, Comput. Methods Appl. Mech. Eng.,
300 (2016), pp. 716 – 733, <https://doi.org/10.1016/j.cma.2015.12.005>.

[24]
C. Lehrenfeld, F. Heimann, J. Preuß, and H. von Wahl, `ngsxfem`:
Add-on to NGSolve for geometrically unfitted finite element
discretizations, J. Open Source Softw., 6 (2021), p. 3237,

[25]
C. Lehrenfeld and M. A. Olshanskii, An Eulerian finite element
method for PDEs in time-dependent domains, ESAIM: M2AN, 53 (2019),
pp. 585–614, <https://doi.org/10.1051/m2an/2018068>.

[26]
C. Lehrenfeld and A. Reusken, Analysis of a Nitsche XFEM-DG
discretization for a class of Two-Phase Mass Transport Problems,
SIAM J. Numer. Anal., 51 (2013), pp. 958–983,

[27]
C. Lehrenfeld and A. Reusken, Analysis of a high order unfitted
finite element method for an elliptic interface problem, IMA J. Numer.
Anal., 38 (2018), pp. 1351–1387,

[28]
Y. Lou and C. Lehrenfeld, Isoparametric unfitted BDF – Finite
element method for PDEs on evolving domains,

[29]
C. Lundholm, On a Space-Time Extended Finite Element
Method for the Solution of a Class of Two-Phase Mass Transport
Problems, PhD thesis, Chalmers University of Technology and University of
Gothenburg, 2021, <https://research.chalmers.se/en/publication/524200>.

[30]
C. Ma, Q. Zhang, and W. Zheng, A high-order fictitious-domain method
for the advection-diffusion equation on time-varying domain,  (2021),

[31]
R. Mittal, J. H. Seo, V. Vedula, Y. J. Choi, H. Liu, H. H. Huang, S. Jain,
L. Younes, T. Abraham, and R. T. George, Computational modeling of
cardiac hemodynamics: Current status and future outlook, J. Comput. Phys.,
305 (2016), pp. 1065–1082, <https://doi.org/10.1016/j.jcp.2015.11.022>.

[32]
M. A. Olshanskii and A. Reusken, Trace finite element methods for
PDEs on surfaces, in Geometrically unfitted finite element methods and
applications, Springer, 2017, pp. 211–258.

[33]
J. Parvizian, A. Düster, and E. Rank, Finite cell method,
Comput. Mech., 41 (2007), pp. 121–133,

[34]
J. Preuß, Higher order unfitted isoparametric space-time FEM
on moving domains, Master's thesis, University of Göttingen,  (2018),

[35]
B. Schott, Stabilized Cut Finite Element Methods for Complex
Interface Coupled Flow Problems, dissertation, Technische Universität
München, München, 2017,

[36]
V. Thomée, Galerkin finite element methods for parabolic
problems, vol. 25 of Springer Series in Computational Mathematics,
Springer-Verlag, Berlin, 1997,

[37]
H. von Wahl, T. Richter, and C. Lehrenfeld, An unfitted Eulerian
finite element method for the time-dependent Stokes problem on moving
domains, IMA J. Numer. Anal.,  (2021),

[38]
S. Zahedi, A space-time cut finite element method with quadrature in
time, in Geometrically Unfitted Finite Element Methods and Applications -
Proceedings of the UCL Workshop 2016, Lecture Notes in Computational Science
and Engineering, Springer, 2018,