*format=pnasformat

# Neural parameter calibration and uncertainty quantification for epidemic
forecasting

Thomas Gaskin1,2,⋆ Tim Conrad3 Grigorios A. Pavliotis2 Christof Schütte3,4

###### Abstract

The recent COVID-19 pandemic has thrown the importance of accurately
forecasting contagion dynamics and learning infection parameters into sharp
focus. At the same time, effective policy-making requires knowledge of the
uncertainty on such predictions, in order, for instance, to be able to ready
hospitals and intensive care units for a worst-case scenario without
needlessly wasting resources. In this work, we apply a novel and powerful
computational method to the problem of learning probability densities on
contagion parameters and providing uncertainty quantification for pandemic
projections. Using a neural network, we calibrate an ODE model to data of the
spread of COVID-19 in Berlin in 2020, achieving both a significantly more
accurate calibration and prediction than Markov-Chain Monte Carlo (MCMC)-based
sampling schemes. The uncertainties on our predictions provide meaningful
confidence intervals e.g. on infection figures and hospitalisation rates,
while training and running the neural scheme takes minutes where MCMC takes
hours. We show convergence of our method to the true posterior on a simplified
SIR model of epidemics, and also demonstrate our method’s learning
capabilities on a reduced dataset, where a complex model is learned from a
small number of compartments for which data is available.

Keywords: Epidemic modelling, neural networks, parameter estimation,
uncertainty quantification, COVID-19

  * 1Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge CB3 0WA, United Kingdom; 2Department of Mathematics, Imperial College London, London SW7 2AZ, United Kingdom; 3Zuse Institute Berlin, Berlin, Germany; 4Department of Mathematics and Computer Science, Freie Universität Berlin, Berlin, Germany

⋆To whom correspondence should be addressed<EMAIL_ADDRESS>
###### Contents

  1. Introduction
  2. Methodology
  3. Modelling the spread of COVID-19 in Berlin
  4. Discussion

## Introduction

The COVID-19 pandemic has underscored the need for comprehensive
epidemiological models. In a crisis, an effective government response is
predicated on such models (1) being formulated using interpretable infection
parameters, (2) being capable of accurately and quickly forecasting the
dynamics of contagion, and (3) meaningfully capturing the uncertainty inherent
in their projections [1, 2]. It is this last point in particular—the inclusion
of uncertainty quantification—that enables informed and transparent cost-
benefit analyses, since it lets policy-makers assess the probable efficacy of
different intervention strategies.

Ordinary differential equations (ODEs) play a key role in the study of
mathematical epidemiology, since they often serve as the foundation for so-
called _compartmental models_. Compartmental models capture the dynamics of
patients progressing through the various stages of a disease—from, e.g., a
susceptible and exposed state all the way to being infected, recovering,
becoming immunised, or falling critically ill. The transition rates between
these compartments model the biological and behavioral drivers of disease
transmission and recovery. Based on Bernoulli’s seminal work in the 1760s [3],
the first compartmental models were introduced by Kermack & McKendrick [4].
Since then, compartmental ODE models have played a crucial role in offering a
profound understanding of the transmission dynamics of infectious illnesses.
However, ODE models are inherently deterministic and thus fail to account for
stochastic characteristics of real-world transmission processes. Over the
years they have therefore been evolved into stochastic variants, based e.g. on
stochastic differential equations (SDEs) [5, 6], and, more recently, agent-
based micro-models (ABMs), incorporating demographic data with full spatial
resolution [7, 8]. Such models have been vital in assisting public health
authorities in predicting outbreaks and implementing efficient disease control
strategies in a variety of instances, including influenza, West Nile virus,
childhood illnesses, SARS-CoV-1, rabies, sexually transmitted infections such
as AIDS, and—most recently—COVID-19 [9, 10, 11, 12]. In this work, we present
a computational method that combines the traditional ODE approach with a
recently developed neural parameter estimation scheme [13, 14] to recover the
model parameters driving the dynamics of a disease from observable data. We
employ a neural network to optimally calibrate a given ODE model to available
data, thereby capturing both the deterministic and stochastic components of
contagion within a single modelling framework.

The existing body of research covers a variety of methods for data-based
parametrisation of ODE models quantifying associated uncertainties, including
hierarchical models, non-parametric techniques, and Bayesian approaches [15,
16, 17, 18]. The ubiquitous Bayesian paradigm again unfolds into a rich
tapestry of techniques, including Markov-Chain Monte Carlo methods (MCMC) [19]
such as Hamiltonian Monte Carlo [20], or samplers based on Langevin dynamics,
such as the Metropolis-adjusted Langevin algorithm (MALA) and its
preconditioned variants [21, 22, 23, 24], to name just a few. The commonality
of all Bayesian parameter estimation and uncertainty quantification schemes is
that they require sampling from a posterior distribution. However, the
sampling paradigm has several major disadvantages: first, since the likelihood
of a parameter is represented by its sampling frequency, high-dimensional
inference can quickly become a costly affair. This is particularly true when
the likelihood of a sample must be computed by solving the underlying ODE.
Random-walk behaviour may render sampling in high dimensions computationally
infeasible, and it is for this reason that modern MCMC samplers are often
strongly gradient-driven. Second, the complex geometry of high-dimensional
distributions may lead to samplers getting caught in local minima, and chains
not reaching stationarity. This at times can be remedied by incorporating the
topology of the distribution into a preconditioner, using e.g. the Hessian or
Fisher information matrices [22]. However, calculating the Riemannian metric
of a high-dimensional space can again become computationally expensive,
especially if such information is only defined locally and thus needs to be
recalculated at each point. Third, almost all sampling strategies require some
sample rejection and burn-in periods, during which samples are discarded to
ensure convergence of the Markov chains. Here again, whenever likelihoods must
be obtained through expensive simulations, having to discard samples means
wasting computational resources.

The strategy proposed in this paper falls within the Bayesian paradigm and yet
circumvents these problems. It thus represents a notable improvement over
existing techniques, both in terms of computational efficiency and accuracy,
and has previously been applied to a diverse set of problems, from estimating
low-dimensional ODE parameters to entire network adjacency matrices. Being
purely gradient-driven, no samples are rejected, a burn-in period is
unnecessary, and the method quickly finds the modes of the distribution. The
neural network parametrises the parameter space without the need for
calculating Riemannian metrics. The loss function contains knowledge of the
model equations, and likelihoods are estimated using simulation. Crucially, by
performing multiple (parallelisable) training iterations, our methodology not
only improves the accuracy of individual parameter estimations, but also
offers a comprehensive representation of the likelihood distribution over the
parameter space and thereby an understanding of the inherent uncertainty—a
vital aspect in the formulation of policies in the face of intrinsic
unpredictability, as is usually the case in a rapidly developing pandemic
scenario.

The main purpose of this article is to demonstrate that this new method is
straightforwardly applicable to applied epidemiological modelling of real-
world data, while going beyond merely enhancing parameter accuracy; instead,
we systematically address the frequently overlooked aspect of uncertainty in
disease forecasting. We begin by presenting the mathematical foundations of
our method on a simple model of epidemics, before revisiting a sophisticated
compartmental model and recalibrating it to observation data of the spread of
COVID-19 in Berlin.

## Methodology

In this section, we delve into the specifics of our methodology, starting with
the underlying ideas and using a simple epidemic model as an example. The
basic concepts have previously been described in [13], but will be summarised
again in the following with a view to applying it to epidemic modelling.

Consider an Itô stochastic differential model with $N$ compartments of the
dynamics of some contagious disease,

$\dfrac{\mathrm{d}\mathbf{y}}{\mathrm{d}t}=f(\mathbf{y}(t);\bm{\Lambda})\,+\sigma(\mathbf{y}(t);\bm{\Lambda})\bm{\xi}_{t}.$
(1)

Here, $\mathbf{y}(t)\in\mathbb{R}^{N}$ is the $N$-dimensional state vector
describing the model at time $t$, $f$ and $\sigma$ are the drift vector and
diffusion matrix, respectively, $\bm{\Lambda}\in\mathbb{R}^{p}$ a vector of
scalar parameters, and $\bm{\xi}_{t}$ an $N$-dimensional white noise process.
Spatial dependence of $\mathbf{y}$, leading to (stochastic) partial
differential equations, is not considered in this work but has been elsewhere
[14]. We note that our method is not dependent on the choice of the stochastic
integral used in the model equations.

Given a time series $\mathbf{T}$ comprising $L$ observations of $\mathbf{y}$,
$\mathbf{T}=(\mathbf{y}_{1},...,\mathbf{y}_{L})$, our goal is to infer the
parameters $\bm{\Lambda}$. To this end we train a neural network
$u_{\theta}:\mathbb{R}^{N\times B}\to\mathbb{R}^{p}$, where the batch size
$B\geq 1$ represents the number of time series steps that are passed as input
and $\bm{\theta}$ the neural network parameters, to produce a parameter
estimate $\bm{\hat{\Lambda}}=(\hat{\lambda}_{1},...,\hat{\lambda}_{p})$ that,
when inserted into the model equations (1), reproduces the observations
$\mathbf{T}$. The neural network is trained using a loss function (such as a
weighted least squares residual)

$J\left(\bm{\hat{\Lambda}}\;\middle|\;\mathbf{T}\right)=J\left(\mathbf{\hat{T}}(\bm{\hat{\Lambda}})\;\middle|\;\mathbf{T}\right),$
(2)

where $\mathbf{\hat{T}}(\bm{\hat{\Lambda}})$ is the time series obtained by
integrating equation (1) using the estimated parameters. . The likelihood of
any estimate is then simply proportional to

$\rho\left(\bm{\hat{\Lambda}}\;\middle|\;\mathbf{T}\right)\propto e^{-J}.$ (3)

As $\bm{\hat{\Lambda}}=\bm{\hat{\Lambda}}(\bm{\theta})$, we may calculate the
gradient $\nabla_{\bm{\theta}}J$ and use it to optimise the internal
parameters of the neural net using a backpropagation method of choice.
Calculating $\nabla_{\bm{\theta}}J$ thus requires differentiating the
predicted time series $\mathbf{\hat{T}}$, and thereby the system equations,
with respect to $\bm{\hat{\Lambda}}$. In other words: the loss function
contains knowledge of the dynamics of the model. Finally, the true data is
once again input to the neural net to produce a new parameter estimate
$\bm{\hat{\Lambda}}$, and the cycle starts afresh. Note that the gradient
descent is not applied to the parameter space directly, but to its
reparametrisation as a neural network. The optimal reparametrisation can be
obtained through optimising the hyperparameters of the neural network
architecture, that is, the depth, size, and structure of the layers, the use
of biases, and the choice of activation functions.

As the net trains, it traverses the parameter space, calculating a loss at
each point. Unlike in MCMC, the posterior density in our approach is not
constructed by considering the frequency with which each point is sampled, but
rather calculated directly via the loss function at that point (cf. [13]).
This entirely eliminates the need for rejection sampling or a burn-in time: at
each point, the true value of the likelihood is obtained, and sampling a
single point multiple times gives no additional information, leading to a
significant improvement in computational speed. Since the stochastic sampling
process is entirely driven by the gradient of $J$, the regions of high
probability are typically found much more rapidly than with a random sampler,
leading to a high sample density around the modes of the target distribution.

We thus track the neural network’s path through the parameter space and gather
the loss values it calculates along the way. Multiple training runs can be
performed in parallel, and each chain is terminated once it reaches a stable
minimum. The likelihood for each parameter is given by

$\rho(\hat{\lambda}_{i}\;|\;\mathbf{T})=\int
p\left(\bm{\hat{\Lambda}}\;\middle|\;\mathbf{T}\right)\mathrm{d}\bm{\hat{\Lambda}}_{-i},$
(4)

where the ${-i}$ subscript indicates omission of the $i$-th component in
$\bm{\hat{\Lambda}}$ in the integration. In high dimensions, calculating the
joint distribution can become computationally infeasible, and we can
approximate the likelihood function by calculating the two-dimensional joint
density of the parameter estimate and the likelihood, $p(\lambda_{i},e^{-J})$
and then integrating over the likelihood,

$\rho(\hat{\lambda}_{i}\;|\;\mathbf{T})\approx\int
p\left(\hat{\lambda}_{i},e^{-J}\right)\mathrm{d}(e^{-J}).$ (5)

By Bayes’ rule, the posterior marginal is then

$p(\hat{\lambda}_{i}\;|\;\mathbf{T})=\rho(\hat{\lambda}_{i}\;|\;\mathbf{T})\times\pi^{0}(\hat{\lambda}_{i})$

with $\pi^{0}$ the prior density [25]. The only prior information available
about the values of the parameters is that they are positive, hence in the
following we will always assume uniform priors on $\mathbb{R}_{+}$. Running
multiple chains in parallel increases the sampling density on the domain,
ensuring convergence to the posterior distribution in the limit of infinitely
many chains, independently of the choice of the prior.

(a) Marginals $p(\beta\;|\;\mathbf{T})$

(b) Marginals $p(\tau\;|\;\mathbf{T})$

(c) Marginals $p(\alpha\;|\;\mathbf{T})$

(d) True data and neural prediction

Figure 1: (a)–(c) Marginal densities on $(\beta,\tau,\alpha)$ for noisy SIR
data, obtained from the neural scheme (blue) and the MALA sampler (pink). The
ground truth (dotted line) was calculated using a simple grid search on
$(\beta,\tau)\in[0,1]\times[1,30]$ with 10.000 grid points. (d) Predicted
average time series $\langle\mathbf{\hat{T}}\rangle$ for the S (lightgreen), I
(red), and R compartments. Shown are the predictions (solid line) and standard
deviation (shaded area) generated by drawing 10.000 samples from the predicted
joint distribution of $(\beta,\tau)$ using the neural scheme and solving the
noiseless ODE model eq. (6) each. Also shown are the true data for each
compartment (dots). The neural network was trained for 100 epochs from 300
different initial conditions, and 50 MALA chains were run until stationarity
was reached, with a thinning factor of 5. Here, stationarity is defined via a
Gelman-Rubin statistic of below 1.2, see fig. S2 in the SI. Both the neural
and MCMC samplers are parallelised.

#### Illustration of a simple epidemic model.

We first consider a synthetic example of a simple model of epidemics with
three compartments, before turning to the main goal of this paper, which is to
present an extensive analysis on a dataset of the spread of COVID-19 in
Berlin. Given observations of susceptible (S), infected (I), and recovered (R)
agents, assume the dynamics of the epidemic are given by

$\displaystyle\mathrm{d}\text{S}=-\beta\text{SI}\mathrm{d}t-\sigma\text{I}\mathrm{d}W_{\mathrm{S}}$
$\displaystyle\mathrm{d}\text{I}=(\beta\text{S}-\tau^{-1})\text{I}\mathrm{d}t+\sigma\text{I}\mathrm{d}W_{\mathrm{I}}$
(6) $\displaystyle\mathrm{d}\text{R}=\tau^{-1}\text{I}\mathrm{d}t,$

where $(\beta,\tau,\sigma)\in\mathbb{R}^{3}_{+}$ are the infection, recovery,
and noise parameters respectively, and $W_{\mathrm{S}},\,W_{\mathrm{I}}$ are
independent Wiener processes. We generate noisy observations of the time
series $\mathbf{T}(t)=(\mathrm{S}(t),\mathrm{I}(t),\mathrm{R}(t))$ using the
parameters $\beta=0.2$, $\tau=14$, and $\sigma=0.1$ (cf. fig. 1d), and try to
recover the marginal densities on $\beta$ and $\tau$ given the observations,
that is $\rho(\beta\;|\;\mathbf{T})$ and $\rho(\tau\;|\;\mathbf{T})$. The
ground truth marginals we obtain by running a grid search on $(\beta,\tau)$
and calculating the likelihood at each grid point. We train the neural net
using the loss function

$J=\|\mathbf{\hat{T}}-\mathbf{T}\|_{2}^{2}=-\log\rho\left(\bm{\Lambda}\;\middle|\;\mathbf{T}\right)$
(7)

and compare the predicted marginals to those generated by a preconditioned
Metropolis-adjusted Langevin scheme (MALA) [22, 23]. The preconditioned MALA
draws samples $\bm{\hat{\Lambda}}^{i}$ from the distribution

$\displaystyle\bm{\hat{\Lambda}}^{i+1}$
$\displaystyle\sim\dfrac{\epsilon^{i}}{2}\Bigg{[}$ $\displaystyle
G(\bm{\hat{\Lambda}}^{i})\bigg{(}\nabla\log\rho\left(\bm{\hat{\Lambda}}^{i}\;\middle|\;\mathbf{T}\right)$
$\displaystyle+\dfrac{L}{B}\sum_{k=1}^{L}\nabla\log\rho\left(\mathbf{T}_{k}\;\middle|\;\bm{\hat{\Lambda}}^{i}\right)\bigg{)}+\Gamma(\bm{\hat{\Lambda}}^{i})\Bigg{]}$
$\displaystyle+G^{\frac{1}{2}}(\bm{\hat{\Lambda}}^{i})\mathcal{N}(0,\epsilon^{i}\mathbf{I}).$
(8)

Here, $\epsilon^{i}$ is the step size at iteration $i$, $G(\bm{\Lambda})$ a
preconditioning matrix as given in [23], and
$\Gamma(\mathbf{x})=\sum_{j}\frac{\partial G_{ij}}{\partial x_{i}}$. The
series $\\{\epsilon^{i}\\}$ is chosen to be a decaying series with coefficient
$\alpha<1$. The algorithm is given in [23], which uses an approximate Fisher
information matrix for the preconditioner, thereby maintaining good
computational performance.

A hyperparameter sweep showed a two-layer neural network with 20 nodes per
layer to provide optimal results. We use hyperbolic tangent activation
functions on all except the last layer, where we use the modulus $|\cdot|$.
The net is optimised using Adam [26] and a learning rate of 0.002. Best
results were obtained using a batch size of $B=L$ (batch gradient descent).
Results are shown in figure 1a–b. Both the neural and MCMC approaches are
initialised from multiple different values with sampling chains run in
parallel. The initial values for $\beta$ are drawn from a uniform distribution
on $[0,1]$, and those for $\tau$ from uniform distribution on $[1,30]$ (see
Supplementary Information). The MCMC scheme is run with a burn-in time of 500
steps per chain, and a thinning factor of 5, meaning only every fifth sample
is retained. We see that both the neural and Langevin schemes find the
posterior marginals, though the neural estimates are more accurate.

To additionally test the neural method’s ability to perform sensitivity
analyses, we modify eq. (6) by adding a small perturbation to each term,

$\mathrm{d}\tilde{\text{S}}=\mathrm{d}\text{S}+\dfrac{\mathrm{d}t}{1000+\alpha}$

(and similarly I and R), where $\alpha\in[0,1]$, meaning that $\alpha$ is
essentially irrelevant to the dynamics, and its marginal posterior
approximately uniform on $[0,1]$. As shown in figure 1c, both schemes obtain
the expected result.

|  Neural |  MALA
---|---|---
Hellinger distance $\beta$ |  5e–4 |  2e–2
Hellinger distance $\tau$ |  5e–4 |  2e–2
Hellinger distance $\alpha$ |  9e–3 |  1e–2
Time (2D) |  3 min 57 sec |  29 min 45 sec
Time (3D) |  4 min 13 sec |  43 min 31 sec

Table 1: Hellinger distances (eq. (9)) between the estimated and true
marginals. Generating the ground truth distributions on $(\beta,\tau)$ via a
grid search took 60 minutes. Also shown are the CPU run times for the neural
and MALA schemes, both for estimating the marginals in the two-dimensional
$(\beta,\tau)$ and three-dimensional case ($\beta,\tau,\alpha$).

We more formally compare the distribution accuracy in terms of the Hellinger
distance to the ground truth,

$d_{\rm
H}(\hat{p},p)=\dfrac{1}{2}\int\left(\sqrt{\hat{p}(x)}-\sqrt{p(x)}\right)^{2}\mathrm{d}x,$
(9)

where $\hat{p}$ is the estimated density and $p$ the ground truth. We find the
neural marginals for $\beta$ and $\tau$ to be two orders of magnitude closer
to the ground truth than the MCMC estimates, see table 1. At the same time,
the neural scheme runs about an order of magnitude faster than the MCMC
sampler, a fact that was previously observed in [13, 14]. As mentioned in the
introduction, the MCMC scheme is slowed down by (1) a burn-in period, (2)
redundant samples being discarded in the Metropolis step, and (3) the long
sampling time required to reach stationarity. Each of these drawbacks the
neural scheme manages to avoid.

Using the resulting posteriors $p$, we can now generate a predicted time
series with uncertainty quantification by randomly selecting $n$ samples
$\bm{\hat{\Lambda}}^{i}$, $i=1,...,n$ gathered during the training process,
running the noiseless ODE model eq. (6) with each sample (i.e. using a noise
strength of $\sigma=0$), and calculating the mean densities

$\langle\mathbf{\hat{T}}\rangle=\sum^{N}_{i=1}p(\bm{\hat{\Lambda}}^{i})\mathbf{\hat{T}}(\bm{\hat{\Lambda}}^{i}),$
(10)

and analogously a standard deviation. The predicted densities are shown in
fig. 1d: we see the predicted parameters capture the observed dynamics well,
with all data points lying within a single standard deviation from the mean.

SEIR${\rm Q}_{\rm S}$${\rm Q}_{\rm E}$${\rm Q}_{\rm I}$SYHCDCT$\lambda_{\rm
E}(t)$$\lambda_{\rm I}$$\lambda_{\rm R}$$\lambda_{\rm CT}$$\lambda_{\rm
S}$$\lambda_{\rm CT}$$\lambda_{\rm CT}+\lambda_{\rm SY}$$\lambda_{\rm
R}$$\lambda_{\rm R}$$\lambda_{\rm R}$$\lambda_{\rm D}$$\lambda_{\rm
SY}$$\lambda_{\rm H}$$\lambda_{\rm C}$$\lambda_{\rm SY}$$\lambda_{\rm
I}$$\lambda_{\rm Q}$$\lambda_{\rm R}$

Figure 2: Schematic illustration of the SEIRD+ model, as originally presented
in [12]. Each parameter $\lambda_{i}$ indicates the transition rate between
the respective compartments. S, E, and I are the susceptible, exposed, and
infected agents. Upon contact with an infected agent, each may be contacted by
the contact tracing agency (CT) and ordered to quarantine (${\rm Q}_{\rm S}$,
${\rm Q}_{\rm E}$, ${\rm Q}_{\rm I}$ compartments). $\lambda_{\rm Q}$ models
the rate of compliance with the contact tracing agency’s instructions. SY, H,
and C are the symptomatic, sick, and critically sick agents. Agents from these
as well as the I and ${\rm Q}_{\rm I}$ compartments can recover and transition
to the R compartment, where they are assumed to stay, at least for the period
under consideration ($<$ 9 months). Finally, critically ill patients may die
of the disease (D), though this compartment is not included in the loss
function. We assume the exposure rate $\lambda_{\rm E}$ varies as public
health measures change. The parameter $\lambda_{\rm Q}$ is further assumed to
be a function of $\lambda_{\rm CT}$ and ${\rm CT}$, and is thus not learned;
see Supplementary Information. Figure adapted from [12].

## Modelling the spread of COVID-19 in Berlin

We now turn to a sophisticated model of the spread of COVID-19 in Berlin,
previously studied in [12]. The authors presented an extended version of the
compartmental SEIRD model, modelling—among others—those infected, symptomatic,
sick (i.e. requiring medical attention or hospitalisation), and critically
sick (requiring ICU or otherwise urgent treatment), as well as a contact-
tracing mechanism responsible for notifying those previously in contact with
an infected person, and consigning them to quarantine. A model overview is
presented in fig. 2, and the ODE system is similar in structure to the SIR
model previously studied; see the Supplementary Information or [12] for the
equations. This model was calibrated to data from an agent-based model of
Berlin [27], comprising over 3 million agents, the transport system, the
geography and urban structure, as well as workflow routines and travel
patterns. The compartments obtained from the ABM data are exclusive, meaning
that e.g. a critical agent is not also classified as symptomatic or
hospitalised. The ABM was calibrated to match the case numbers of COVID-19
from February 16 to October 27 2020, and provides estimates of the hidden
infection cases which were not officially recorded but nevertheless driving
hospitalisation and mortality rates. Crucially, it assumes that the official
infection figures recorded by the Robert-Koch Institute are the sum of the SY,
H, and C compartments, and do _not_ contain the I compartment, since in the
early stages of the pandemic asymptomatic cases were not usually detected (cf.
fig. 3).

Figure 3: Evolution of COVID-19 in Berlin for the period from February 16 to
October 27 2020. Shown are the ABM data [27] for the symptomatic,
hospitalised, and critical compartments (orange, red, purple), the sum of all
three (light brown, solid line), as well as the official infection figures
(light brown, dots) [28]. The red period is the calibration period, with the
shades representing varying levels of government restrictions and
correspondingly different exposure levels $\lambda_{\rm E}$: from mid-March,
businesses and factories started closing; in late March, the German government
imposed broad contact restrictions; in early May, schools and kindergartens
started reopening across the country, followed by further loosening of
restrictions in mid-June, before the start of the summer holidays. The blue
period is the projection data on which we evaluate the prediction. The ABM
data only contains a single Q compartment and no CT compartment. It also does
not produce a D compartment, for the reasons given in the text. See
Supplementary Information for details.

Public health measures naturally had an impact on the virus dynamics: on March
12, factories, theatres, and concert halls started closing, and the German
Bundesliga suspended all football games. Ten days later, the federal
government prohibited all gatherings of more than two people, exempting single
households. These measures lasted through April 2020, after which retail,
schools, and kindergartens gradually started reopening, the federal government
giving states broad autonomy to set their own policies on May 6. [29, 30].
Starting in mid-June, restrictions on social gatherings were further relaxed
across the country. These measures will primarily have affected the
population’s exposure to the virus, and we thus assume that the exposure
parameter $\lambda_{\rm E}$ is piecewise linear on the intervals [Feb 16, Mar
12, Mar 22, May 6, June 15, Oct 27]. This is not taking into account virus
mutations changing its infectivity or lethality. The vector of parameters
$\bm{\Lambda}$ we wish to estimates is thus 13-dimensional, comprising the 9
parameters shown in figure 2, one of which is time-dependent.

Figure 4: Marginals on the parameters $\bm{\Lambda}=(\lambda_{\rm
S},...,\lambda_{\rm CT})$. Shown are the neural marginals (blue, left side)
and MCMC estimates (pink, right side), which in both cases were smoothed using
a Gaussian kernel. Also shown are the means (green dots) and modes (yellow
dots) of the marginals. We employ a three-layer neural network with 20 neurons
per layer and sigmoid activation functions on all but the last layer, where we
again use the absolute value function.

(a) Neural

(b) MALA

Figure 5: Comparison of the neural calibration results (left) and the MCMC
calibration results (right) for the symptomatic, hospitalized, and critical
compartments. Red lines are the true data, green lines the prediction using
the estimated mean of the joint density, calculated by drawing 1000 samples
from the joint distribution. The green shaded areas represent one standard
deviation. The blue shaded area is the test period for which projections are
generated. Calibration results for the remaining compartments are shown in
fig. S3 in the Supplementary Information.

We split the data into a training period of 200 days, spanning the period up
to September 3 (‘calibration period’), and a test period, spanning the
remaining eight weeks until October 27 (‘projection period’). We employ a deep
neural network with 3 layers, 20 neurons per layer, and the sigmoid as an
activation function on all but the last layer, where we again use the modulus.
The batch size $B$ is again equal to the length of the time series ($B=L$).
The number of agents in each compartment span several orders of magnitude,
from millions of susceptible agents down to hundreds of hospitalised or
critical agents. Using the simple loss function from the previous example, eq.
(7), would result in only the largest compartments being fitted accurately,
since their residuals dominate the loss. We therefore scale each compartment’s
contribution to the total loss,

$J=\sum_{i}\alpha_{i}\|\mathbf{\hat{T}}_{i}-\mathbf{T}_{i}\|_{2}^{2},$ (11)

by choosing coefficients $\alpha_{i}$ that ensure all summands are roughly of
equal magnitude:

$\alpha_{i}^{-1}=\int_{0}^{L}\mathbf{T}_{i}(t)\mathrm{d}t,$ (12)

where $L$ denotes the length of the time series. Note that, due to
inconsistencies in the official mortality statistics for the early period of
the pandemic (arising from the difficulty of discerning ‘death by COVID’ from
‘death with COVID’) we do not fit the ${\rm D}$ compartment.

Assuming uniform priors on all parameters $\lambda_{i}$, we show the marginal
posteriors alongside the MCMC estimates in figure 4. The means and modes of
the distributions are indicated in the plot. In general, the neural posteriors
are more sharply peaked and unimodal than the MCMC estimates, though the
expectation values and modes tend to roughly match. A low sensitivity to
$\lambda_{\rm S}$ is unsurprising given the large pool of susceptible agents
during the early stages of the pandemic in a city of 3.6 million. One notable
exception is $\lambda_{\rm SY}$, where the neural network predicts a much
lower rate than MCMC. The marginals on $\lambda_{\rm S}$ and $\lambda_{\rm
CT}$ are fairly broad, indicating low sensitivity to these transition rates.
Also observe that the neural marginals for the exposure parameter are
unimodal, with the means and modes obeying $\hat{\lambda}_{\rm
E,0}\gg\hat{\lambda}_{\rm\rm E,1}>\hat{\lambda}_{\rm\rm
E,4}>\hat{\lambda}_{\rm\rm E,3}>\hat{\lambda}_{\rm\rm E,2}$. This is
consistent with the level of government restrictions imposed, and it is
interesting to note that the measures taken between March 12 and March 22
already reduced the exposure rate by two-thirds. This pattern does not hold
for the MCMC estimates.

Figure 6: Results on a reduced training dataset, consisting only of the
symptomatic, hospitalized, and critical compartments.

As before, we now draw $n$ samples from the joint densities to produce a mean
time series $\langle\mathbf{\hat{T}}\rangle$ for each compartment. We compare
the quality of the fit on each compartment, both on the training period and
the projection period, using the $L^{2}$ residual

$r_{i}^{2}=\langle\mathbf{\hat{T}}_{i}(t)-\mathbf{T}_{i}(t)\rangle^{2}_{t},\
i\in\\{{\rm S},...,{\rm Q}\\},$ (13)

with the expectation value $\langle\cdot\rangle_{t}$ taken over time. In order
to circumvent having to calculate the full 13-dimensional joint distribution
$p(\bm{\hat{\Lambda}})$, we simply select $n=1000$ random samples previously
collected during training. Shown in fig. 4 is the true data (red), the mean
prediction $\langle\mathbf{\hat{T}}\rangle$ (green), and one standard
deviation (shaded area). The neural approach visibly calibrates the training
data to a higher accuracy than MCMC, consequently also achieving a better fit
on the projection period (blue shaded area). The residuals $r_{i}$ are given
in table 2. The neural scheme achieves an average calibration error of 25%,
representing a 50% improvement over MCMC, and a projection error of 18%, a
fourfold improvement over MCMC. At the same time, the neural network runs an
order of magnitude faster.

| Neural error | MALA error
---|---|---
|  Calibration |  Projection |  Calibration |  Projection
S |  0.0004 |  0.002 |  0.003 |  0.02
E |  0.36 |  0.17 |  0.32 |  0.67
I |  0.23 |  0.11 |  0.31 |  0.67
R |  0.22 |  0.04 |  0.49 |  0.67
SY |  0.25 |  0.12 |  0.35 |  0.97
H |  0.32 |  0.28 |  0.52 |  0.85
C |  0.28 |  0.24 |  0.64 |  0.52
Q |  0.34 |  0.11 |  0.43 |  1.33
Avg. |  0.25 |  0.18 |  0.38 |  0.71
Time | 8 mins 45 secs | 88 mins 1 sec
Table 2: Calibration and projection error of the neural and MALA schemes on
the different compartments. The calibration period refers to the test period
from Feburary 16 to September 3, while the projection period is the remaining
two months until October 27. The error is given in terms of the $L^{2}$
residual $r_{i}$ (eq. (13)). Also shown: CPU run time.

Until now we have assumed full knowledge of all compartments in the model, but
this is only due to a sophisticated and computationally expensive ABM running
in the background, laboriously calibrated to official data. Without this
machinery, the available data only covers the symptomatic, hospitalised, and
critical compartments [31, 28, 32]. We thus re-train the model on these
compartments only, and assess the calibration quality. Results are shown in
fig. 6. The neural network still calibrates each compartment with an average
error of 0.33, an 18% reduction compared to the full model. However, the
prediction error is 0.98, an almost five-fold decrease compared to the full
model. Simultaneously, the model’s confidence in its predictions decreases
visibly: the full model is thus required to make accurate predictions with
high confidence.

## Discussion

In this article we presented a method to optimally fit compartmental infection
models to observations of infection spreading. We applied a novel and powerful
computational method to the problem of learning probability densities on
contagion parameters and providing uncertainty quantification for epidemic
projections. This new methodology can be utilized based on simulation data
coming from finer scale models (as demonstrated herein) or directly on
observational data. In the former case, it allows finding optimal surrogate
models to fine-scale infection models in order to use them in optimal control
or multi-objective optimization approaches as in [12].

The strategy proposed in this paper represents a notable improvement over
conventional MCMC or Langevin sampling methods due to its superior
computational accuracy and efficiency in estimating the parameters of the
model and their uncertainty. The comprehensive understanding of uncertainty it
provides is vital to developing effective policy responses when faced with
intrinsic unpredictability. We mention, in particular, that the exploration of
a relatively high-dimensional parameter space using MCMC can be extremely
expensive, especially when—as is the case here—the likelihood must be obtained
via simulation. Furthermore, the marginals fig. 4 strongly indicate that the
parameter space is highly non-convex with many different local minima trapping
the MCMC sampler and significantly increasing the mixing times. In our
analysis, we also noted the slow convergence of the Gelman-Rubin statistic for
the Langevin sampler—see fig. S4 in the appendix. Overall, in our experiments
our method delivered a 10-fold decrease in compute times, while calibrating
and predicting the spread of COVID-19 significantly more accurately. In our
numerical experiments, even state-of-the-art MCMC schemes fail to fully
explore the parameter space, in particular if the model contains redundant
parameters. Our proposed method, by contrast, does not suffer from this
drawback.

Recently, new alternative sampling methods for Bayesian uncertainty
quantification and inversion have been proposed; one example is the Affine
Invariant Langevin Dynamics (ALDI) [33, 34], a modification of the Ensemble
Kalman Sampler [35] with significant theoretical advantages over
preconditioned MALA (such as affine invariance and convergence in total
variation to the posterior, at least for convex problems). Further schemes
include Hamiltonian Monte Carlo [20] and the bouncy particle sampler [36]. A
comprehensive comparison of the various sampling schemes and their relative
benefits for calibrating epidemiological models will be the subject of future
work.

Data, materials, and software availability Code data can be found under
https://github.com/ThGaskin/NeuralABM. It is easily adaptable to new models
and ideas. The code uses the utopya package111utopia-project.org,
utopya.readthedocs.io/en/latest [37, 38] to handle simulation configuration
and efficiently read, write, analyse, and evaluate data. This means that the
model can be run by modifying simple and intuitive configuration files,
without touching code. Multiple training runs and parameter sweeps are
automatically parallelised. The neural core is implemented using
pytorch222pytorch.org.  Acknowledgements TG was funded by the University of
Cambridge School of Physical Sciences VC Award via DAMTP and the Department of
Engineering, and supported by EPSRC grants EP/P020720/2 and EP/R018413/2. TC
and CS were funded by the Deutsche Forschungsgemeinschaft (DFG) under
Germany’s Excellence Strategy through grant EXC-2046 _The Berlin Mathematics
Research Center_ MATH+ (project no. 390685689) and via the grant MODUS-COVID
by the German Federal Ministry for Education and Research. GP is partially
supported by the Frontier Research Advanced Investigator Grant ERC grant
Machine-aided general framework for fluctuating dynamic density functional
theory. The authors would also like to thank Hanna Wulkow for her support in
acquiring the ABM data.  Author contributions All authors designed the
research and wrote the paper. TG designed and performed the numerical
experiments and wrote the code.

## References

  * [1] JM Carcione, JE Santos, C Bagaini, J Ba, A Simulation of a COVID-19 Epidemic Based on a Deterministic SEIR Model.  Frontiers in Public Health 8 (2020).
  * [2] R McCabe, et al., Communicating uncertainty in epidemic models.  Epidemics 37, 100520 (2021).
  * [3] D Bernoulli, Essai d’une nouvelle analyse de la mortalité causée par la petite vérole, et des avantages de l’inoculation pour la prévenir.  Académie Royale des Sciences pp. 1–45 (1760).
  * [4] WO Kermack, AG McKendrick, A contribution to the mathematical theory of epidemics.  Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character 115, 700–721 (1927).
  * [5] E Tornatore, SM Buccellato, P Vetro, Stability of a stochastic SIR system.  Phys.A 354, 111–126 (2005).
  * [6] A Gray, D Greenhalgh, L Hu, X Mao, J Pan, A Stochastic Differential Equation SIS Epidemic Model.  SIAM Journal on Applied Mathematics 71, 876–902 (2011).
  * [7] SA Müller, et al., Predicting the effects of COVID-19 related interventions in urban settings by combining activity-based modelling, agent-based simulation, and mobile phone data.  PLOS ONE 16, 1–32 (2021).
  * [8] CC Kerr, et al., Covasim: An agent-based model of COVID-19 dynamics and interventions.  PLOS Computational Biology 17, 1–32 (2021).
  * [9] F Brauer, P Van den Driessche, J Wu, LJ Allen, Mathematical epidemiology.  (Springer) Vol. 1945, (2008).
  * [10] K Hanke, et al., Reconstruction of the genetic history and the current spread of HIV-1 subtype A in Germany.  J. Virol. 93 (2019).
  * [11] MT Meehan, et al., Modelling insights into the COVID-19 pandemic.  Paediatric Respiratory Reviews 35, 64–69 (2020).
  * [12] H Wulkow, et al., Prediction of Covid-19 spreading and optimal coordination of counter-measures: From microscopic to macroscopic models to Pareto fronts.  PLOS ONE 16, e0249676 (2021).
  * [13] T Gaskin, GA Pavliotis, M Girolami, Neural parameter calibration for large-scale multi-agent models.  Proceedings of the National Academy of Sciences 120 (2023).
  * [14] T Gaskin, GA Pavliotis, M Girolami, Inferring networks from time series: a neural approach.  arXiv 2303.18059 [cs.LG] (2023).
  * [15] R Ghanem, D Higdon, H Owhadi, Handbook of uncertainty quantification.  (Springer), (2017).
  * [16] HJ Zimmermann, An application-oriented view of modeling uncertainty.  European Journal of operational research 122, 190–198 (2000).
  * [17] JC Helton, JD Johnson, WL Oberkampf, An exploration of alternative approaches to the representation of uncertainty in model predictions.  Reliability Engineering & System Safety 85, 39–71 (2004).
  * [18] Y Lin, MA Stadtherr, Validated solutions of initial value problems for parametric odes.  Applied Numerical Mathematics 57, 1145–1162 (2007).
  * [19] GO Roberts, JS Rosenthal, General state space Markov chains and MCMC algorithms.  Probability Surveys [electronic only] 1, 20–71 (2004).
  * [20] MD Hoffman, A Gelman, The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.  J. Mach. Learn. Res. 15, 1593–1623 (2014).
  * [21] GO Roberts, RL Tweedie, Exponential convergence of Langevin distributions and their discrete approximations.  Bernoulli 2, 341 – 363 (1996).
  * [22] M Girolami, B Calderhead, Riemann manifold Langevin and Hamiltonian Monte Carlo methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73, 123–214 (2011).
  * [23] C Li, C Chen, D Carlson, L Carin, Preconditioned Stochastic Gradient Langevin Dynamics for Deep Neural Networks in Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence, AAAI’16.  (AAAI Press), pp. 1788–1794 (2016).
  * [24] N Wulkow, R Telgmann, KD Hungenberg, C Schütte, M Wulkow, Deterministic and Stochastic Parameter Estimation for Polymer Reaction Kinetics I: Theory and Simple Examples.  Macromolecular Theory and Simulations 30, 2100017 (2021).
  * [25] AM Stuart, Inverse problems: A Bayesian perspective.  Acta Numerica 19, 451–559 (2010).
  * [26] DP Kingma, J Ba, Adam: A Method for Stochastic Optimization.  arXiv 1412.6980 [cs.LG] (2014).
  * [27] MODUS Covid, MATSim-Episim, Version 2020-November-12 (https://covid-sim.info/2020-11-12/secondLockdown) (2020).
  * [28] Robert Koch-Institut and Bundesamt für Kartographie und Geodäsie, Robert Koch-Institut: Fallzahlen in Deutschland (https://npgeo-corona-npgeo-de.hub.arcgis.com/datasets/6d78eb3b86ad4466a8e264aa2e32a2e4_0/about) (2021).
  * [29] ARD Tagesschau Online, Chronologie: Drei Jahre Pandemie (https://www.tagesschau.de/inland/gesellschaft/corona-pandemie-rueckblick-101.html) (2023).
  * [30] German Federal Health Ministry (Bundesministerium für Gesundheit), Coronavirus-Pandemie: Was geschah wann? (https://www.bundesgesundheitsministerium.de/coronavirus/chronik-coronavirus) (2023).
  * [31] Statistik Berlin-Brandenbug, Schwerpunkt Corona (https://www.statistik-berlin-brandenburg.de/corona) (2022).
  * [32] Robert-Koch Institut, DIVI Intensivregister (Register of ICU bed occupancy) (https://github.com/robert-koch-institut/Intensivkapazitaeten_und_COVID-19-Intensivbettenbelegung_in_Deutschland) (2023).
  * [33] A Garbuno-Inigo, N Nüsken, S Reich, Affine Invariant Interacting Langevin Dynamics for Bayesian Inference.  SIAM Journal on Applied Dynamical Systems 19, 1633–1658 (2020).
  * [34] A Garbuno-Inigo, F Hoffmann, W Li, AM Stuart, Interacting Langevin diffusions: gradient structure and ensemble Kalman sampler.  SIAM J. Appl. Dyn. Syst. 19, 412–441 (2020).
  * [35] S Reich, S Weissmann, Fokker–Planck Particle Systems for Bayesian Inference: Computational Approaches.  SIAM/ASA Journal on Uncertainty Quantification 9, 446–482 (2021).
  * [36] A Bouchard-Côté, SJ Vollmer, A Doucet, The Bouncy Particle Sampler: A Nonreversible Rejection-Free Markov Chain Monte Carlo Method.  Journal of the American Statistical Association 113, 855–867 (2018).
  * [37] L Riedel, B Herdeanu, H Mack, Y Sevinchan, J Weninger, Utopia: A comprehensive and collaborative modeling framework for complex and evolving systems.  Journal of Open Source Software 5, 2165 (2020).
  * [38] Y Sevinchan, B Herdeanu, J Traub, dantro: a python package for handling, transforming, and visualizing hierarchically structured data.  Journal of Open Source Software 5, 2316 (2020).

Supporting Information

## Methodology

### Initialising the neural network

We can initialise the neural network to a particular distribution on its
output by drawing random variables $\bm{\Lambda}_{0}$ from the target
distribution and training the neural network to output the target value before
the run. This can be done using the simple loss function

$J=\|\bm{\hat{\Lambda}}-\bm{\Lambda}_{0}\|_{2},$ (14)

which typically only takes a few seconds. When running the neural network
multiple times from different initialisations the resulting distribution of
initial values will approach the target distribution as the number of chains
goes to infinity. In fig. S1 we show the initial values for the SIR parameters
$(\beta,\tau,\alpha)$ for 100 different initialisations; as is visible, the
distributions are approximately uniform on the intervals as specified.

Figure S1: Initial value distribution on $\beta$, $\tau$, and $\alpha$ for the
neural network. Each parameter is initialised to a uniform distribution.

### Running the MCMC

We initialise the MCMC with the same initial distribution as the neural
network. For the SIR experiments, we run 50 randomly initialised chains for
10.000 steps each. Figure S2 shows the Gelman-Rubin statistic over time for
the 3-dimensional SIR experiment. Note that the value for $\alpha$ does not
converge to below 1.2, which seems to indicate the MCMC sampler getting caught
in local minima, preventing complete mixing.

(a) Infection parameter $\beta$

(b) Recovery parameter $\tau$

(c) Redundant parameter $\alpha$

Figure S2: Gelman-Rubin statistics for the MCMC chains for the parameters
$(\beta,\tau,\alpha)$. The red dotted line indicates a value of 1.2.

## Modelling the spread of COVID-19 in Berlin

(a) Susceptible

(b) Exposed

(c) Infected

(d) Recovered

(e) Quarantined (all)

Figure S3: Calibration and projection results pertaining to table 2 and fig.
5; shown are the results for the S, E, I, R, and Q compartments. The Q
compartment is given by the sum of the ${\rm Q}_{\rm S}$, ${\rm Q}_{\rm E}$,
and ${\rm Q}_{\rm I}$ compartments.

The ODE model used to model the spread of COVID-19 in Berlin was previously
presented in [12], and is visualised in figure 2 in the main manuscript:

$\displaystyle\dfrac{\mathrm{d}{\rm S}}{\mathrm{d}t}$
$\displaystyle=-\lambda_{\rm E}{\rm SI}+\lambda_{\rm S}{\rm Q}_{\rm
S}-\lambda_{\rm Q}{\rm S}$ $\displaystyle\dfrac{\mathrm{d}{\rm
D}}{\mathrm{d}t}=\lambda_{\rm D}{\rm C}$ $\displaystyle\dfrac{\mathrm{d}{\rm
E}}{\mathrm{d}t}$ $\displaystyle=\lambda_{\rm E}{\rm SI}-\lambda_{\rm E}{\rm
I}-\lambda_{\rm Q}{\rm E}$ $\displaystyle\dfrac{\mathrm{d}{\rm Q}_{\rm
E}}{\mathrm{d}t}=-\lambda_{\rm I}{\rm Q}_{\rm E}+\lambda_{\rm Q}{\rm E}$
$\displaystyle\dfrac{\mathrm{d}{\rm I}}{\mathrm{d}t}$
$\displaystyle=\lambda_{\rm I}{\rm E}-\lambda_{\rm R}{\rm I}-\lambda_{\rm
SY}{\rm I}-\lambda_{\rm Q}{\rm I}$ $\displaystyle\dfrac{\mathrm{d}{\rm Q}_{\rm
I}}{\mathrm{d}t}=\lambda_{\rm I}{\rm Q}_{\rm E}+\lambda_{\rm Q}{\rm
I}-\lambda_{\rm SY}{\rm Q}_{\rm I}-\lambda_{\rm R}{\rm Q}_{\rm I}$
$\displaystyle\dfrac{\mathrm{d}{\rm R}}{\mathrm{d}t}$
$\displaystyle=\lambda_{\rm R}\left({\rm I}+{\rm SY}+{\rm H}+{\rm C}+{\rm
Q}_{\rm I}\right)$ $\displaystyle\dfrac{\mathrm{d}{\rm
CT}}{\mathrm{d}t}=\lambda_{\rm SY}{\rm I}+(\lambda_{\rm CT}-\lambda_{\rm
Q})({\rm S}+{\rm E}+{\rm I})$ $\displaystyle\dfrac{\mathrm{d}{\rm
SY}}{\mathrm{d}t}$ $\displaystyle=\lambda_{\rm SY}\left({\rm Q}_{\rm I}+{\rm
I}\right)-\lambda_{\rm R}{\rm SY}-\lambda_{\rm H}{\rm SY}$
$\displaystyle\dfrac{\mathrm{d}{\rm Q}_{\rm S}}{\mathrm{d}t}=-\lambda_{\rm
S}{\rm Q}_{\rm S}+\lambda_{\rm Q}{\rm S}$ $\displaystyle\dfrac{\mathrm{d}{\rm
H}}{\mathrm{d}t}$ $\displaystyle=\lambda_{\rm H}{\rm SY}-\lambda_{\rm R}{\rm
H}-\lambda_{\rm C}{\rm H}$ $\displaystyle\dfrac{\mathrm{d}\lambda_{\rm
Q}}{\mathrm{d}t}=\lambda_{\rm q}\lambda_{\rm CT}{\rm CT}$
$\displaystyle\dfrac{\mathrm{d}{\rm C}}{\mathrm{d}t}$
$\displaystyle\lambda_{\rm C}{\rm H}-\lambda_{\rm R}{\rm C}-\lambda_{\rm
D}{\rm C}$ (15)

Here, each compartment $i={\rm S},...,{\rm CT}$ represents the density of
agents in that compartment, i.e. ${\rm S}={\rm s}/N$, where ${\rm s}$ is the
number of agents in that compartment, and $N$ is the total population of
Berlin (ca. $3.6\times 10^{6}$). $\lambda_{\rm q}$ is the rate of people going
into quarantine after being called by the contact tracing agency, and is
assumed to be independent of counter-measures. Following [12], we set
$\lambda_{\rm q}=10.25$. The ABM data does not contain a deceased compartment,
as the Robert-Koch Institute changed its classification of ‘deceased’ in
mid-2020, in order to better reflect that not all patients diagnosed with
COVID-19 at the time of their death had actually died _of_ the disease. The
mortality rates for this period are unreliable, and are excluded from the
model. A good approximation (leading to similar results) would be to set ${\rm
D}={\rm C}/3$. Also not given are the number of contact-traced agents (CT),
and furthermore the ABM data does not distinguish between the three different
quarantined compartments. The loss function used to train the neural network
thus becomes

$\displaystyle J=\sum_{t=1}^{L}\big{[}$ $\displaystyle\alpha_{\rm S}(\hat{{\rm
S}}-{\rm S})^{2}+\alpha_{\rm E}(\hat{{\rm E}}-{\rm E})^{2}+\alpha_{\rm
I}(\hat{{\rm I}}-{\rm I})^{2}+\alpha_{\rm R}(\hat{{\rm R}}-{\rm
R})^{2}+\alpha_{\rm SY}(\hat{{\rm SY}}-{\rm SY})^{2}$
$\displaystyle+\alpha_{\rm H}(\hat{{\rm H}}-{\rm H})^{2}+\alpha_{\rm
C}(\hat{{\rm C}}-{\rm C})^{2}+\alpha_{\rm Q}((\hat{\rm{Q}}_{\rm
I}+\hat{\rm{Q}}_{\rm S}+\hat{\rm{Q}}_{\rm E})-{\rm Q})^{2}\big{]},$

with the coefficients $\alpha_{i}$ given in eq. (12). In fig. S3 we show the
calibrations and projections for the S, E, I, R, and Q compartments, which
were not shown in the main manuscript.

In figure S4 we show the Gelman-Rubin statistic for the Langevin sampler. For
most of the parameters, the chains remain far from the desired value of below
1.2, which would broadly indicate convergence of the chains. Instead, most of
the values remain stubbornly high. This again is likely due to the highly non-
convex structure of the parameter space, with many different local minima
trapping the sampler and preventing mixing (cf. fig. 4).

Figure S4: Gelman-Rubin statistics for the MCMC scheme for each of the
parameters used to calibrate the ODE model eq. (15). The red dotted value
indicates a value of 1.2.


# NNLO QCD corrections to $\operatorname{W}\kern-1.0ptb\kern 1.0pt\bar{b}$
production at the LHC

Heribertus Bayu Hartanto<EMAIL_ADDRESS>Cavendish Laboratory,
University of Cambridge, Cambridge CB3 0HE, United Kingdom  Rene Poncelet
<EMAIL_ADDRESS>Cavendish Laboratory, University of Cambridge,
Cambridge CB3 0HE, United Kingdom  Andrei Popescu<EMAIL_ADDRESS>Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, United
Kingdom  Simone Zoia<EMAIL_ADDRESS>Dipartimento di Fisica and Arnold-
Regge Center, Università di Torino, and INFN, Sezione di Torino, Via P. Giuria
1, I-10125 Torino, Italy

###### Abstract

We compute theoretical predictions for the production of a W-boson in
association with a bottom-quark pair at hadron colliders at next-to-next-to-
leading order (NNLO) in QCD, including the leptonic decay of the W-boson,
while treating the bottom quark as massless. This calculation constitutes the
very first $2\to 3$ process with a massive external particle to be studied at
such a perturbative order. We derive an analytic expression for the required
two-loop five-particle amplitudes in the leading colour approximation
employing finite-field methods. Numerical results for the cross section and
differential distributions are presented for the Large Hadron Collider at
$\sqrt{s}=8$ TeV. We observe an improvement of the perturbative convergence
for the inclusive case and for the prediction with a jet veto upon the
inclusion of the NNLO QCD corrections.

††preprint: CAVENDISH-HEP-22/05

##  I Introduction

Studying vector boson production in association with multi-jet final states at
the Large Hadron Collider (LHC) offers a wide variety of interesting
phenomenological explorations. In particular, the production of a W-boson in
association with bottom quark ($b$) jets is very interesting both from
experimental and theoretical perspectives. It is crucial to scrutinise
experimental signatures for both $\operatorname{W}\kern-1.0pt+1b$ jet and
$\operatorname{W}\kern-1.0pt+2b$ jets, in order to test our knowledge of the
strong interaction at high energies and improve our modelling of bottom-quark
jets at the LHC. The cross sections for both the signatures have been measured
at the Tevatron [1, 2] and LHC [3, 4, 5, 6]. While the
$\operatorname{W}\kern-1.0pt+1b$ jet signature provides a fundamental probe of
the bottom-quark parton distribution functions (PDFs), the
$\operatorname{W}\kern-1.0pt+2b$ jets final state constitutes an irreducible
background to many important reactions studied at the LHC, such as the Higgs-
strahlung process
($\operatorname{p}\kern-1.0pt\operatorname{p}\kern-1.0pt\to\operatorname{W}\kern-1.0pt\operatorname{H}\kern-1.0pt\
(\operatorname{H}\kern-1.0pt\to b\kern 1.0pt\bar{b})$) and single top
production
($\operatorname{p}\kern-1.0pt\operatorname{p}\kern-1.0pt\to\bar{b}\,t\
(t\to\operatorname{W}\kern-1.0ptb)$), as well as many beyond the Standard
Model (BSM) searches. Moreover, the $\operatorname{W}\kern-1.0pt+b$ jets
processes are interesting from a theoretical point of view as they are a
perfect testing ground to study the different ways of treating the $b$-quark.
In particular, the choice of whether to take its mass and presence in the PDF
into account, leads to two disparate computational schemes: the four- (4FS)
and five-flavour number schemes (5FS).

In this letter we compute the NNLO QCD corrections to W-boson production in
association with a bottom-quark pair, which we henceforth call
$\operatorname{W}\kern-1.0ptb\kern 1.0pt\bar{b}$, including the leptonic decay
of the W-boson ($\operatorname{W}\kern-1.0pt\to\ell\nu)$. We work in 5FS, thus
treating the bottom quark, and additionally the charged leptons, as massless
particles. They contribute to the $\operatorname{W}\kern-1.0pt+2b$ jets
signature, as well as $\operatorname{W}\kern-1.0pt+1b$ jet production when at
least one $b$-jet is tagged. Extensive studies of
$\operatorname{W}\kern-1.0ptb\kern 1.0pt\bar{b}$ production at NLO QCD
accuracy [7, 8, 9, 10, 11, 12] indicate a poor perturbative behaviour at such
order (i.e. the corrections are large and the scale uncertainties do not
improve with respect to the leading order predictions for inclusive final
state) due to the opening of the $qg$-initiated channel. Several efforts to
assess corrections beyond NLO were done by including additional jet radiations
[13, 14]. It is clear that a fully-fledged NNLO QCD prediction is mandatory to
improve the perturbative convergence of $\operatorname{W}\kern-1.0ptb\kern
1.0pt\bar{b}$ production.

Improvement of theoretical precision is also a critical component of the
progress in Particle Physics, as we enter the precision LHC era with the
upcoming Run 3 and high-luminosity phases. We have seen spectacular
breakthroughs in perturbative QCD calculations in the recent years, with a
number of $2\to 3$ processes computed at NNLO QCD accuracy for fully massless
final states [15, 16, 17, 18, 19, 20]. This success stems from both the
advancements in the scattering amplitude computations and the developments of
NNLO subtraction schemes.

The analytic computation of the required two-loop five-particle amplitudes is
one of the main bottlenecks towards achieving NNLO QCD accuracy for $2\to 3$
processes. However, the progress for five-particle processes with a single
massive external particle has been spectacular recently. All planar two-loop
five-particle Feynman integrals are now available analytically [21, 22, 23,
24] in terms of bases of special functions which substantially simplify
computation of the finite remainders, and allow for an extremely efficient
numerical evaluation [25, 26]. Partial results for the non-planar integral
families have also become available recently [27, 28, 29]. This progress
resulted in a number of two-loop amplitudes computed at leading colour [25,
30, 31, 32]. In this work we derive the analytic form of the leading colour
two-loop amplitude contributing to
$\operatorname{W}\kern-1.0pt\,(\to\ell\nu)\,b\kern 1.0pt\bar{b}$ production,
and employ it to compute a number of observables for this process at NNLO in
QCD. Our computation marks a significant precision-calculation milestone,
since it represents the very first prediction to be derived for a $2\to 3$
process involving a massive final state.

This letter is structured as follows. We begin by discussing the derivation of
the required two-loop matrix elements and the tools used to perform the cross
section computation, and then present phenomenological results for the cross
section and a number of interesting differential distributions.

##  II Calculation

We consider
$\operatorname{p}\kern-1.0pt\operatorname{p}\kern-1.0pt\to\ell^{+}\nu\,b\kern
1.0pt\bar{b}$ production (with $\ell=e$ or $\mu$) up to
$\mathcal{O}(\alpha^{2}\alpha_{s}^{4})$. The calculation has been performed
within the Stripper framework, a C++ implementation of the four-dimensional
formulation of the sector-improved residue subtraction scheme [33, 34, 35].
The tree-level matrix elements are supplied by the AvH library [36], while the
one-loop matrix elements are provided by the OpenLoops package [37, 38]. We
compute the double virtual contribution $\mathcal{V}^{(2)}$ in the leading
colour approximation for

$u(p_{1})+\bar{d}(p_{2})\to
b(p_{3})+\bar{b}(p_{4})+\ell^{+}(p_{5})+\nu(p_{6})\,.$ (1)

It consists of the (colour and helicity summed) two-loop and one-loop squared
matrix elements,

$\mathcal{V}^{(2)}=\sum_{\text{col.}}\sum_{\text{hel.}}\left\\{2\,\mathrm{Re}\big{[}M^{(0)*}\mathcal{F}^{(2)}\big{]}+\big{|}\mathcal{F}^{(1)}\big{|}^{2}\right\\}\,,$
(2)

where $M^{(0)}$ is the tree-level amplitude, and $\mathcal{F}^{(L)}$ is the
$L$-loop finite remainder. We decompose $\mathcal{V}^{(2)}$ at leading colour
into

$\displaystyle\mathcal{V}^{(2)}_{\mathrm{LC}}$
$\displaystyle=\mathcal{V}^{(2),1}+\frac{n_{f}}{N_{c}}\mathcal{V}^{(2),n_{f}}+\frac{n_{f}^{2}}{N_{c}^{2}}\mathcal{V}^{(2),n_{f}^{2}}\,,$
(3)

where $n_{f}$ is the number of massless closed fermion loops. We note that the
leading-colour approximation is only enforced in the scale-independent part of
the double virtual contribution,

$\displaystyle\mathcal{V}^{(2)}(\mu_{R}^{2})$
$\displaystyle=\mathcal{V}^{(2)}_{\mathrm{LC}}(s_{12})+\sum_{i=1}^{4}c_{i}\ln^{i}\bigg{(}\frac{\mu_{R}^{2}}{s_{12}}\bigg{)}\,,$
(4)

where $s_{ij}=(p_{i}+p_{j})^{2}$, and the kinematic-dependent coefficients
$c_{i}$ are expressed in terms of full colour lower-order matrix elements.

The analytic computation of the two-loop amplitude follows closely Ref. [25],
with modifications implemented to incorporate the decay of the W-boson. Since
the QCD corrections do not apply to the
$\operatorname{W}\kern-1.0pt\to\ell\nu$ decay, we can separate the 6-point
squared amplitude $\mathcal{M}^{(2)}_{6}$ into the product of the 5-point
W-production squared amplitude $\mathcal{M}^{(2)}_{5\mu\nu}$ and the leptonic
tensor $\mathcal{D}^{\mu\nu}$,

$\displaystyle\mathcal{M}_{6}^{(2)}=\mathcal{M}_{5\mu\nu}^{(2)}\;\mathcal{D}^{\mu\nu}\;|P(s_{56})|^{2}\,,$
(5)

where
$P(s)=1/(s-M_{\operatorname{W}\kern-1.0pt}^{2}+iM_{\operatorname{W}\kern-1.0pt}\Gamma_{\operatorname{W}\kern-1.0pt})$
is the W-boson propagator factor. We perform tensor decomposition on the
5-point W-production squared amplitude,

$\displaystyle\mathcal{M}_{5}^{(2)\mu\nu}=\sum_{i=1}^{16}a^{(2)}_{i}v_{i}^{\mu\nu}\,,$
(6)

using $\\{p_{1},p_{2},p_{3},p_{5}+p_{6}\\}$ as the spanning basis to build the
$v_{i}^{\mu\nu}$ basis tensors [39, 40]. The coefficients $a^{(2)}_{i}$ can be
determined by contracting Eq. (6) with $v_{i\mu\nu}$ and inverting the
resulting linear system of equations. The analytic form of the contracted
squared amplitudes $v_{i\mu\nu}\mathcal{M}_{5}^{(2)\mu\nu}$ was derived using
finite-field reconstruction methods within the FiniteFlow framework [41, 42].
We expressed them in terms of the special functions of Ref. [26] and rational
coefficients, which we simplified using MultivariateApart [43] and Singular
[44]. We further implemented these amplitudes in C++ for a fast numerical
evaluation. We evaluate the special functions using the PentagonFunction++
library [26]. Our analytic result is validated numerically against the
$\operatorname{W}\kern-1.0pt+4$ quarks helicity amplitudes derived in Ref.
[31] at the level of the helicity-summed squared finite remainder. The
complete analytic expression is included in ancillary files.

##  III Phenomenology

We present numerical results for the LHC center-of-mass energy $\sqrt{s}=8$
TeV, focusing on the $\operatorname{W}\kern-1.0pt^{+}(\to\ell^{+}\nu)\,b\kern
1.0pt\bar{b}$ final state. The Standard Model input parameters are

$\displaystyle M_{\operatorname{W}\kern-1.0pt}$
$\displaystyle=80.351972\;\mathrm{GeV},$
$\displaystyle\Gamma_{\operatorname{W}\kern-1.0pt}=2.0842989\;\mathrm{GeV},$
$\displaystyle M_{\operatorname{Z}\kern-1.0pt}$
$\displaystyle=91.153481\;\mathrm{GeV},$
$\displaystyle\Gamma_{\operatorname{Z}\kern-1.0pt}=2.4942665\;\mathrm{GeV},$
$\displaystyle G_{F}$ $\displaystyle=1.16638\cdot 10^{-5}\;\mathrm{GeV}^{-2},$

from which the electromagnetic coupling $\alpha$ can be derived within the
$G_{\mu}$ scheme. We assume a diagonal CKM matrix and employ the
`NNPDF31_as_0118` PDF sets [45] with its perturbative order matching that of
the corresponding calculations. Since we treat the bottom quark as massless,
we need a flavour-sensitive jet algorithm to define the flavoured jets in an
infrared-safe way. The partons are clustered into a jet using the
flavour-$k_{\text{T}}$ jet algorithm [46] with $R=0.5$. The jets (including
$b$-jets) and charged leptons are required to fulfill the following event-
selection criteria [6]:

$p_{\text{T},\ell}>30\;\mathrm{GeV},\;\;|\eta_{\ell}|<2.1,\;\;p_{\text{T},j}>25\;\mathrm{GeV},\;\;|\eta_{j}|<2.4.$
(7)

The renormalisation and factorisation scales are set to a common value
$\mu_{R}=\mu_{F}=H_{\text{T}}$, with

$H_{\text{T}}=E_{\text{T}}(\ell\nu)+p_{\text{T}}(b_{1})+p_{\text{T}}(b_{2})\,,$
(8)

where $b_{1}$, $b_{2}$ are correspondingly the hardest and second hardest
$b$-flavoured (either $b$ or $\bar{b}$) jets. Unless otherwise specified, the
scale uncertainties are obtained using the 7-point scale variation, where
$\mu_{R}$ and $\mu_{F}$ are varied by a factor of 2 around $H_{\text{T}}$,
while satisfying the $1/2\leq\mu_{R}/\mu_{F}\leq 2$ constraint.

Based on the number of jets required in the final states, we can define the
following configurations for the NLO and NNLO predictions:

  * •

inclusive (inc): at least 2 $b$-jets;

  * •

exclusive (exc): exactly 2 $b$-jets and no other jets.

Naïve scale variation of the exclusive prediction may lead to an
underestimation of the scale uncertainties [47]. Hence, for the exclusive
configuration, we use also the uncorrelated prescription of Ref. [47], in
addition to the 7-point scale variation.

|  inclusive [fb] |  $\mathcal{K}_{\mathrm{inc}}$ |  exclusive [fb] |  $\mathcal{K}_{\mathrm{exc}}$
---|---|---|---|---
$\sigma_{\mathrm{LO}}$ |  $213.2(1)^{+21.4\%}_{-16.1\%}$ |  - |  $213.2(1)^{+21.4\%}_{-16.1\%}$ |  -
$\sigma_{\mathrm{NLO}}$ |  $362.0(6)^{+13.7\%}_{-11.4\%}$ |  1.7 |  $249.8(4)^{+3.9(+27)\%}_{-6.0(-19)\%}$ |  1.17
$\sigma_{\mathrm{NNLO}}$ |  $445(5)^{+6.7\%}_{-7.0\%}$ |  1.23 |  $267(3)^{+1.8(+11)\%}_{-2.5(-11)\%}$ |  1.067
Table 1:  Fiducial cross sections for
$\operatorname{p}\kern-1.0pt\operatorname{p}\kern-1.0pt\to\ell^{+}\nu
b\bar{b}$ production at the LHC with $\sqrt{s}=8$ TeV at LO, NLO and NNLO for
both inclusive (inc) and exclusive (exc) final states. The corresponding
$\mathcal{K}$ factor is defined as
$\mathcal{K}=\sigma_{\mathrm{N^{(n)}LO}}/\sigma_{\mathrm{N^{(n-1)}LO}}$. The
statistical errors are shown for the central predictions. Scale uncertainties
for the exclusive predictions are provided using both the standard 7-point
scale variation and uncorrelated prescription of Ref. [47]. The latter is
quoted inside parentheses in the error estimates.

In Table 1, we present numerical results for the fiducial cross section for
the inclusive and exclusive configurations at different perturbative orders.
As observed in the previous studies [9, 14], the NLO QCD corrections are large
in the case of the inclusive phase space. In our calculation this amounts to
about $70\%$ corrections. The jet veto in the exclusive selection reduces the
NLO QCD corrections to a moderate $17\%$. A similar observation holds at NNLO
QCD, where we find a positive correction of $23\%$ in the inclusive and
$6.7\%$ in the exclusive case. The NNLO QCD corrections are smaller than the
NLO QCD corrections in both cases indicating perturbative convergence. In that
respect, by using the scale dependence as the canonical way to estimate the
uncertainties from missing higher orders, we conclude that theoretical
uncertainty reduces with inclusion of higher order terms. However, for the
inclusive phase space, the NLO corrections are significantly larger than the
LO scale dependence. The situation at NNLO QCD slightly improves, but the
corrections are still only barely covered by the NLO scale band. For the
exclusive case, the NLO corrections are within the LO band, however the
estimated uncertainty from the 7-point scale variation is comparatively small,
only $5\%$. The NNLO corrections here are also smaller, but are well outside
the NLO scale uncertainty, indicating that the NLO scale dependence is
underestimated. This motivates the alternative prescription of Ref. [47] to
estimate theory uncertainties, taking into account the jet veto effect. The
uncertainties resulting from this prescription are shown in the parentheses
and are significantly larger. The higher order corrections fall well within
the uncertainty bands, implying that this method is more reliable, but also
quite conservative.

The double virtual corrections, which have been included only in the leading
colour approximation, deserve an additional comment. For the inclusive setup,
we find that the contribution of Eq. (4) to the cross section is about $5\%$.
In the exclusive case, the Born configurations are unaffected by the jet veto,
but a fraction of the hard radiative corrections are suppressed. This leads to
an enhancement of the sensitivity to the double virtual matrix element, which
contributes $\sim\\!10\%$ of the fiducial cross section in this case. The
naïve expectation for the subleading colour effects is that they are about
$10\%$ of the double virtual matrix element, implying that potential
corrections to the fiducial cross section would be about $1\%$ ($0.5\%$) for
the exclusive (inclusive) case.

Figure 1: The charged lepton’s transverse momentum distribution. The upper
panel shows the absolute predictions for the inclusive and exclusive selection
at different perturbative orders. The middle panel shows the inclusive cross
sections as a ratio with respect to the central NLO prediction, with the
coloured bands indicating the 7-point scale variation. The lower panel shows
the same ratio for the exclusive configuration. Here, the coloured bands
correspond to the decorrelated scale variation, and the hashed bands to the
standard 7-point variation. The vertical bars indicate the statistical
uncertainty.

Turning to the differential distributions, we present the transverse momentum
of the charged lepton, $p_{\text{T},\ell}$, in Fig. 1, for the inclusive, as
well as exclusive, phase space selection. Focusing on the perturbative
convergence of the spectrum, we can draw similar conclusions as for the
fiducial cross section. In the inclusive case, we find sizeable NNLO QCD
corrections of $\sim\\!20\%$, which are barely contained in the NLO
uncertainty. The corrections have a tendency to increase at higher energies,
being the largest around $p_{\text{T},\ell}\approx 100$ GeV, similarly to the
NLO corrections. For the exclusive phase space, we find positive corrections
of about $7\%$ for low $p_{\text{T},\ell}$, and negative corrections of order
$\sim\\!10\%$ for $p_{\text{T},\ell}>100$ GeV. Again, we observe that the
decorrelated prescription to estimate the uncertainty is more reliable.

The next two distributions characterise the $b\kern 1.0pt\bar{b}$ system. In
Fig. 2, we show the transverse momentum of the $b\kern 1.0pt\bar{b}$ system,
$p_{\text{T},b\bar{b}}$. In terms of perturbative corrections we find a
similar trend as for the charged lepton transverse momentum. Additionally, the
absolute distributions highlight that the inclusive spectrum is, in general,
harder than the exclusive case, confirming the intuition that the jet veto
suppresses additional large transverse momentum jets. In the case of exclusive
phase space, this differential distribution can be understood as a proxy for
the W transverse momentum.

Figure 2: Transverse momentum distribution of the $b\bar{b}$ system. Same
layout as in Fig. 1.

The distribution of the invariant mass of the $b\kern 1.0pt\bar{b}$ system,
$M_{b\bar{b}}$, is shown in Fig. 3. This observable is interesting when
considering the QCD process $\operatorname{W}\kern-1.0ptb\kern 1.0pt\bar{b}$
as background to the Higgs-strahlung process
$\operatorname{W}\kern-1.0pt\operatorname{H}\kern-1.0pt\,(\to b\bar{b})$.
Around the Higgs mass we can see that the NNLO QCD corrections are about
$20\%$ in the inclusive selection and only $\sim\\!5\%$ in the exclusive case.
By comparing the two prescriptions for estimating the uncertainty, we see that
around the Higgs mass the 7-point prescription implies a 2-3 times smaller
uncertainty than the decorrelated method.

Figure 3: Invariant mass distribution of the $b\bar{b}$ system. Same layout as
in Fig. 1.

The reader is invited to find our results for other observables in the
auxilliary files to this publication.

##  IV Conclusions

We presented fiducial and differential cross sections for the
$\operatorname{W}\kern-1.0ptb\kern 1.0pt\bar{b}$ process at the LHC with 8 TeV
center-of-mass energy. This includes the computation of the double virtual
amplitudes in the leading colour approximation with incorporated decay of the
W-boson.

We addressed the observation of large NLO QCD corrections in this process, and
found that the NNLO QCD corrections are significantly smaller. We observe a
significant reduction of the scale dependence, which indicates perturbative
convergence. We discussed the behaviour of the jet-vetoed cross section, which
exhibits much smaller corrections but suffers from accidental cancellations in
the scale dependence, rendering the theory uncertainty estimates from
canonical scale variation unreliable. At NNLO accuracy we validated the
alternative prescription of Ref. [47] for estimating the theory uncertainties.

This work constitutes the first NNLO QCD calculation for a $2\to 3$ process
including a massive final state particle. Studying this class of processes at
such accuracy is of the utmost importance for the physics programme of the
LHC. However, the steep requirements for amplitudes involving many loops and
high multiplicities used to put them beyond the reach of computational
capabilities. Our results demonstrate that the door of precision phenomenology
is finally open for these processes as well.

###### Acknowledgements.

The authors would like to thank Michał Czakon for making the Stripper library
available to us, Vasily Sotnikov for help in the comparison with the results
of Ref. [31], and Simon Badger and Alexander Mitov for many inspiring
discussions and useful comments on the draft. This project received funding
from the European Union’s Horizon 2020 research and innovation programmes New
level of theoretical precision for LHC Run 2 and beyond (grant agreement No
683211), and High precision multi-jet dynamics at the LHC (grant agreement No
772009). HBH was partially supported by STFC consolidated HEP theory grant
ST/T000694/1. SZ gratefully acknowledges the computing resources provided by
the Max Planck Institute for Physics and by the Max Planck Computing & Data
Facility. AP is also supported by the Cambridge Trust and Trinity College
Cambridge. RP acknowledges the support from the Leverhulme Trust and the Isaac
Newton Trust, as well as the use of the DiRAC Cumulus HPC facility under Grant
No. PPSP226.

## References

  * Abazov  _et al._ [2005] V. M. Abazov _et al._ (D0), A Search for $Wb\bar{b}$ and $WH$ Production in $p\bar{p}$ Collisions at $\sqrt{s}=1.96$ TeV, Phys. Rev. Lett. 94, 091802 (2005), arXiv:hep-ex/0410062 .
  * Abazov  _et al._ [2013] V. M. Abazov _et al._ (D0), Measurement of the $p\bar{p}\to W+b+X$ production cross section at $\sqrt{s}=1.96$ TeV, Phys. Lett. B 718, 1314 (2013), arXiv:1210.0627 [hep-ex] .
  * Aad  _et al._ [2012] G. Aad _et al._ (ATLAS), Measurement of the cross section for the production of a $W$ boson in association with $b^{-}$ jets in $pp$ collisions at $\sqrt{s}=7$ TeV with the ATLAS detector, Phys. Lett. B 707, 418 (2012), arXiv:1109.1470 [hep-ex] .
  * Aad  _et al._ [2013] G. Aad _et al._ (ATLAS), Measurement of the cross-section for W boson production in association with b-jets in pp collisions at $\sqrt{s}$ = 7 TeV with the ATLAS detector, JHEP 06, 084, arXiv:1302.2929 [hep-ex] .
  * Chatrchyan  _et al._ [2014] S. Chatrchyan _et al._ (CMS), Measurement of the Production Cross Section for a W Boson and Two b Jets in pp Collisions at $\sqrt{s}$=7 TeV, Phys. Lett. B 735, 204 (2014), arXiv:1312.6608 [hep-ex] .
  * Khachatryan  _et al._ [2017] V. Khachatryan _et al._ (CMS), Measurement of the production cross section of a W boson in association with two b jets in pp collisions at $\sqrt{s}=8{\,\mathrm{{TeV}}}$, Eur. Phys. J. C 77, 92 (2017), arXiv:1608.07561 [hep-ex] .
  * Ellis and Veseli [1999] R. K. Ellis and S. Veseli, Strong radiative corrections to W b anti-b production in p anti-p collisions, Phys. Rev. D 60, 011501 (1999), arXiv:hep-ph/9810489 .
  * Febres Cordero  _et al._ [2006] F. Febres Cordero, L. Reina, and D. Wackeroth, NLO QCD corrections to W boson production with a massive b-quark jet pair at the Tevatron p anti-p collider, Phys. Rev. D 74, 034007 (2006), arXiv:hep-ph/0606102 .
  * Febres Cordero  _et al._ [2009] F. Febres Cordero, L. Reina, and D. Wackeroth, W- and Z-boson production with a massive bottom-quark pair at the Large Hadron Collider, Phys. Rev. D 80, 034015 (2009), arXiv:0906.1923 [hep-ph] .
  * Badger  _et al._ [2011] S. Badger, J. M. Campbell, and R. K. Ellis, QCD Corrections to the Hadronic Production of a Heavy Quark Pair and a W-Boson Including Decay Correlations, JHEP 03, 027, arXiv:1011.6647 [hep-ph] .
  * Oleari and Reina [2011] C. Oleari and L. Reina, W +- b $\bar{b}$ production in POWHEG, JHEP 08, 061, [Erratum: JHEP 11, 040 (2011)], arXiv:1105.4488 [hep-ph] .
  * Frederix  _et al._ [2011] R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, R. Pittau, and P. Torrielli, W and $Z/\gamma*$ boson production in association with a bottom-antibottom pair, JHEP 09, 061, arXiv:1106.6019 [hep-ph] .
  * Luisoni  _et al._ [2015] G. Luisoni, C. Oleari, and F. Tramontano, $Wb\overline{b}j$ production at NLO with POWHEG+MiNLO, JHEP 04, 161, arXiv:1502.01213 [hep-ph] .
  * Anger  _et al._ [2018] F. R. Anger, F. Febres Cordero, H. Ita, and V. Sotnikov, NLO QCD predictions for $Wb\bar{b}$ production in association with up to three light jets at the LHC, Phys. Rev. D 97, 036018 (2018), arXiv:1712.05721 [hep-ph] .
  * Chawdhry  _et al._ [2020] H. A. Chawdhry, M. L. Czakon, A. Mitov, and R. Poncelet, NNLO QCD corrections to three-photon production at the LHC, JHEP 02, 057, arXiv:1911.00479 [hep-ph] .
  * Kallweit  _et al._ [2021] S. Kallweit, V. Sotnikov, and M. Wiesemann, Triphoton production at hadron colliders in NNLO QCD, Phys. Lett. B 812, 136013 (2021), arXiv:2010.04681 [hep-ph] .
  * Chawdhry  _et al._ [2021] H. A. Chawdhry, M. Czakon, A. Mitov, and R. Poncelet, NNLO QCD corrections to diphoton production with an additional jet at the LHC, JHEP 09, 093, arXiv:2105.06940 [hep-ph] .
  * Czakon  _et al._ [2021] M. Czakon, A. Mitov, and R. Poncelet, Next-to-Next-to-Leading Order Study of Three-Jet Production at the LHC, Phys. Rev. Lett. 127, 152001 (2021), arXiv:2106.05331 [hep-ph] .
  * Badger  _et al._ [2022a] S. Badger, T. Gehrmann, M. Marcoli, and R. Moodie, Next-to-leading order QCD corrections to diphoton-plus-jet production through gluon fusion at the LHC, Phys. Lett. B 824, 136802 (2022a), arXiv:2109.12003 [hep-ph] .
  * Chen  _et al._ [2022] X. Chen, T. Gehrmann, N. Glover, A. Huss, and M. Marcoli, Automation of antenna subtraction in colour space: gluonic processes, (2022), arXiv:2203.13531 [hep-ph] .
  * Papadopoulos  _et al._ [2016] C. G. Papadopoulos, D. Tommasini, and C. Wever, The Pentabox Master Integrals with the Simplified Differential Equations approach, JHEP 04, 078, arXiv:1511.09404 [hep-ph] .
  * Abreu  _et al._ [2020] S. Abreu, H. Ita, F. Moriello, B. Page, W. Tschernow, and M. Zeng, Two-Loop Integrals for Planar Five-Point One-Mass Processes, JHEP 11, 117, arXiv:2005.04195 [hep-ph] .
  * Canko  _et al._ [2021] D. D. Canko, C. G. Papadopoulos, and N. Syrrakos, Analytic representation of all planar two-loop five-point Master Integrals with one off-shell leg, JHEP 01, 199, arXiv:2009.13917 [hep-ph] .
  * Syrrakos [2021] N. Syrrakos, Pentagon integrals to arbitrary order in the dimensional regulator, JHEP 06, 037, arXiv:2012.10635 [hep-ph] .
  * Badger  _et al._ [2021a] S. Badger, H. B. Hartanto, and S. Zoia, Two-Loop QCD Corrections to $Wb\bar{b}$ Production at Hadron Colliders, Phys. Rev. Lett. 127, 012001 (2021a), arXiv:2102.02516 [hep-ph] .
  * Chicherin  _et al._ [2021] D. Chicherin, V. Sotnikov, and S. Zoia, Pentagon Functions for One-Mass Planar Scattering Amplitudes, (2021), arXiv:2110.10111 [hep-ph] .
  * Papadopoulos and Wever [2020] C. G. Papadopoulos and C. Wever, Internal Reduction method for computing Feynman Integrals, JHEP 02, 112, arXiv:1910.06275 [hep-ph] .
  * Abreu  _et al._ [2022a] S. Abreu, H. Ita, B. Page, and W. Tschernow, Two-loop hexa-box integrals for non-planar five-point one-mass processes, JHEP 03, 182, arXiv:2107.14180 [hep-ph] .
  * Kardos  _et al._ [2022] A. Kardos, C. G. Papadopoulos, A. V. Smirnov, N. Syrrakos, and C. Wever, Two-loop non-planar hexa-box integrals with one massive leg, (2022), arXiv:2201.07509 [hep-ph] .
  * Badger  _et al._ [2021b] S. Badger, H. B. Hartanto, J. Kryś, and S. Zoia, Two-loop leading-colour QCD helicity amplitudes for Higgs boson production in association with a bottom-quark pair at the LHC, JHEP 11, 012, arXiv:2107.14733 [hep-ph] .
  * Abreu  _et al._ [2022b] S. Abreu, F. Febres Cordero, H. Ita, M. Klinkert, B. Page, and V. Sotnikov, Leading-color two-loop amplitudes for four partons and a W boson in QCD, JHEP 04, 042, arXiv:2110.07541 [hep-ph] .
  * Badger  _et al._ [2022b] S. Badger, H. B. Hartanto, J. Kryś, and S. Zoia, Two-loop leading colour helicity amplitudes for $W^{\pm}\gamma+j$ production at the LHC, (2022b), arXiv:2201.04075 [hep-ph] .
  * Czakon [2010] M. Czakon, A novel subtraction scheme for double-real radiation at NNLO, Phys. Lett. B 693, 259 (2010), arXiv:1005.0274 [hep-ph] .
  * Czakon and Heymes [2014] M. Czakon and D. Heymes, Four-dimensional formulation of the sector-improved residue subtraction scheme, Nucl. Phys. B 890, 152 (2014), arXiv:1408.2500 [hep-ph] .
  * Czakon  _et al._ [2019] M. Czakon, A. van Hameren, A. Mitov, and R. Poncelet, Single-jet inclusive rates with exact color at $\mathcal{O}$ (${\alpha}_{s}^{4}$), JHEP 10, 262, arXiv:1907.12911 [hep-ph] .
  * Bury and van Hameren [2015] M. Bury and A. van Hameren, Numerical evaluation of multi-gluon amplitudes for High Energy Factorization, Comput. Phys. Commun. 196, 592 (2015), arXiv:1503.08612 [hep-ph] .
  * Buccioni  _et al._ [2018] F. Buccioni, S. Pozzorini, and M. Zoller, On-the-fly reduction of open loops, Eur. Phys. J. C 78, 70 (2018), arXiv:1710.11452 [hep-ph] .
  * Buccioni  _et al._ [2019] F. Buccioni, J.-N. Lang, J. M. Lindert, P. Maierhöfer, S. Pozzorini, H. Zhang, and M. F. Zoller, OpenLoops 2, Eur. Phys. J. C 79, 866 (2019), arXiv:1907.13071 [hep-ph] .
  * Chen [2021] L. Chen, A prescription for projectors to compute helicity amplitudes in D dimensions, Eur. Phys. J. C 81, 417 (2021), arXiv:1904.00705 [hep-ph] .
  * Peraro and Tancredi [2021] T. Peraro and L. Tancredi, Tensor decomposition for bosonic and fermionic scattering amplitudes, Phys. Rev. D 103, 054042 (2021), arXiv:2012.00820 [hep-ph] .
  * Peraro [2016] T. Peraro, Scattering amplitudes over finite fields and multivariate functional reconstruction, JHEP 12, 030, arXiv:1608.01902 [hep-ph] .
  * Peraro [2019] T. Peraro, FiniteFlow: multivariate functional reconstruction using finite fields and dataflow graphs, JHEP 07, 031, arXiv:1905.08019 [hep-ph] .
  * Heller and von Manteuffel [2022] M. Heller and A. von Manteuffel, MultivariateApart: Generalized partial fractions, Comput. Phys. Commun. 271, 108174 (2022), arXiv:2101.08283 [cs.SC] .
  * Decker  _et al._ [2021] W. Decker, G.-M. Greuel, G. Pfister, and H. Schönemann, Singular 4-2-1 — A computer algebra system for polynomial computations, http://www.singular.uni-kl.de (2021).
  * Ball  _et al._ [2017] R. D. Ball _et al._ (NNPDF), Parton distributions from high-precision collider data, Eur. Phys. J. C 77, 663 (2017), arXiv:1706.00428 [hep-ph] .
  * Banfi  _et al._ [2006] A. Banfi, G. P. Salam, and G. Zanderighi, Infrared safe definition of jet flavor, Eur. Phys. J. C 47, 113 (2006), arXiv:hep-ph/0601139 .
  * Stewart and Tackmann [2012] I. W. Stewart and F. J. Tackmann, Theory Uncertainties for Higgs and Other Searches Using Jet Bins, Phys. Rev. D 85, 034011 (2012), arXiv:1107.2117 [hep-ph] .


#  $p$-adic non-commutative analytic subgroup theorem

Duc Hiep Pham

###### Abstract.

In this paper, we formulate and prove the so-called $p$-adic non-commutative
analytic subgroup theorem. This result is seen as the $p$-adic analogue of a
recent theorem given by Yafaev in [11].

2020 Mathematical Subject Classification: 14L10 (22E35, 11F85, 11J81)
Keywords: algebraic groups, analytic subgroup theorem, $p$-adic,
transcendence.

##  1\. Introduction

It is known that the analytic subgroup theorem is considered as one of the
most powerful theorems in complex transcendental number theory. The theorem
was established by Wüstholz in the 1980’s based on a very deep auxiliary
result on multiplicity estimates on group varieties (see [8] and [9]). To
present the theorem, we start with $G$ a commutative algebraic group defined
over $\overline{\mathbb{Q}}$, that is a smooth quasi-projective variety
defined over $\overline{\mathbb{Q}}$ with a commutative group law for which
the composition map $G\times G\rightarrow G,\;(g,h)\mapsto g\circ h$ and the
inverse map $G\rightarrow G,\;g\mapsto g^{-1}$ are regular morphisms between
algebraic varieties defined over $\overline{\mathbb{Q}}$. Then the set
$G({\mathbb{C}})$ of complex points of $G$ is a complex Lie group, and one has
the exponential map
$\exp_{G({\mathbb{C}})}:\mathrm{Lie}(G({\mathbb{C}}))\rightarrow
G({\mathbb{C}})$ of $G$. We say that an element
$u\in\mathrm{Lie}(G({\mathbb{C}}))$ is an algebraic point of the exponential
map of $G$ if $\exp_{G({\mathbb{C}})}(u)\in G(\overline{\mathbb{Q}})$. The
following theorem is a direct consequence of [7, Theorem 1], and it is still
called the analytic subgroup theorem (see also [6, Theorem 1.2]).

######  Theorem 1.1 (Wüstholz).

Let $G$ be a connected commutative algebraic group defined over
$\overline{\mathbb{Q}}$. Let $u\in\mathrm{Lie}(G({\mathbb{C}}))$ be an
algebraic point of the exponential map of $G$ and $V_{u}$ the smallest
$\overline{\mathbb{Q}}$-vector subspace of $\mathrm{Lie}(G)$ such that $u$
lies in the complex vector space
$V_{u}\otimes_{\overline{\mathbb{Q}}}{\mathbb{C}}$. Then there exists a
connected algebraic subgroup $H_{u}$ of $G$ defined over
$\overline{\mathbb{Q}}$ satisfying $\mathrm{Lie}(H_{u})=V_{u}$.

Many results in transcendence theory can be deduced from Wüstholz analytic
subgroup theorem as consequences by choosing suitable commutative algebraic
groups. For instance, the theorem generalizes Baker’s famous result on linear
forms in logarithms and its elliptic analogue in full generality (see [1,
Chapter 6]). It also implies some important results concerning linear
independence on abelian varieties between periods (see [6, Chapter 1]).
Furthermore, Wüstholz himself has recently obtained a nice application on
elliptic and abelian periods spaces (see[10]).

It is able to extend Theorem 1.1 to the non-commutative case. Note that in
this case, there is still the exponential map between the Lie algebra and the
Lie group as the commutative case. The following theorem has been recently
proved by Yafaev (see [11]).

######  Theorem 1.2 (Yafaev).

Let $G$ be a connected algebraic group defined over $\overline{\mathbb{Q}}$.
Let $u\in\mathrm{Lie}(G({\mathbb{C}}))$ be an algebraic point of the
exponential map of $G$ and $V_{u}$ the smallest $\overline{\mathbb{Q}}$-vector
subspace of $\mathrm{Lie}(G)$ such that $u$ lies in the complex vector space
$V_{u}\otimes_{\overline{\mathbb{Q}}}{\mathbb{C}}$. Then there exists a
connected algebraic subgroup $H_{u}$ of $G$ defined over
$\overline{\mathbb{Q}}$ which is commutative and satisfies
$\mathrm{Lie}(H_{u})=V_{u}$.

In 2015, Fuch and Pham obtained the $p$-adic analytic subgroup theorem (see
[3]), and the aim of this paper is to give a such $p$-adic analogue for
Theorem 1.2 above. In the $p$-adic setting, similar to [3, Theorem 2.2], we
also use the $p$-adic logarithm map which is still valid in the non-
commutative case. Fix a prime number $p$, let $\overline{\mathbb{Q}}_{p}$ be
the field of $p$-adic numbers. Denote by ${\mathbb{C}}_{p}$ the completion of
the algebraic closure of $\overline{\mathbb{Q}}_{p}$ (with respect to the
normalized $p$-adic absolute value). Let $G$ now be an algebraic group defined
over $\overline{\mathbb{Q}}$. Consider the set $G({\mathbb{C}}_{p})$ of
${\mathbb{C}}_{p}$-points of $G$. This is a Lie group over ${\mathbb{C}}_{p}$.
Denote by $G({\mathbb{C}}_{p})_{f}$ the set of $x$ in $G({\mathbb{C}}_{p})$
satisfying there exists a strictly increasing sequence $(n_{i})$ of positive
integers for which $x^{n_{i}}$ tends to the identity element of
$G({\mathbb{C}}_{p})$ as $i$ tends to infinity. It follows from [2, Chapter
III, 7.6] that there is a ${\mathbb{C}}_{p}$-analytic map
$\log_{G({\mathbb{C}}_{p})}:G({\mathbb{C}}_{p})_{f}\rightarrow\mathrm{Lie}(G({\mathbb{C}}_{p}))$,
and it is called the $p$-adic logarithm map of $G$. Then the $p$-adic non-
commutative analytic subgroup theorem is formulated as follows.

######  Theorem 1.3.

Let $G$ be a connected algebraic group defined over $\overline{\mathbb{Q}}$ of
positive dimension, and
$\log_{G({\mathbb{C}}_{p})}:G({\mathbb{C}}_{p})_{f}\rightarrow\mathrm{Lie}(G({\mathbb{C}}_{p}))$
the $p$-adic logarithm map of $G$. Let $\gamma\in G({\mathbb{C}}_{p})_{f}$ be
an algebraic point of $G(\overline{\mathbb{Q}})$ with
$\log_{G(\mathbb{C}_{p})}(\gamma)\neq 0$ and $V_{\gamma}$ the smallest
$\overline{\mathbb{Q}}$-vector subspace of $\mathrm{Lie}(G)$ such that the
$p$-adic vector space
$V_{\gamma}\otimes_{\overline{\mathbb{Q}}}{\mathbb{C}}_{p}$ containing the
point $\log_{G(\mathbb{C}_{p})}(\gamma)$. Then there exists a connected
commutative algebraic subgroup $H_{\gamma}\subseteq G$ defined over
$\overline{\mathbb{Q}}$ of positive dimension satisfying $\gamma\in
H_{\gamma}(\overline{\mathbb{Q}})$ and $\mathrm{Lie}(H_{\gamma})=V$.

The proof of our theorem follows closely that of Yafaev’s theorem. We also
reduce the general case to the commutative case and then apply the $p$-adic
analytic subgroup theorem. The main difference is that, in the $p$-adic
setting, the $p$-adic exponential map is only defined locally on an open
subgroup of the $p$-adic Lie algebra of the given algebraic group. However, by
using the $p$-adic logarithm map and by taking a certain power of the
algebraic point, we are able to still keep the important similar arguments as
in the complex case, and based on this to get the derised result.

##  2\. The p-adic exponential map and related commutative algebraic subgroups

In this section, we first briefly recall some background on exponential and
logarithm maps over $p$-adic fields. The main reference we follow here is [2].
For an algebraic group $G$ defined over $\overline{\mathbb{Q}}$, by [2,
Chapter III, 7.2, 7.6] there exist an open subgroup $U_{p}$ of
$\mathrm{Lie}(G({\mathbb{C}}_{p}))$ and an analytic map
$\exp_{G({\mathbb{C}}_{p})}:U_{p}\rightarrow G({\mathbb{C}}_{p})$ such that
$\exp_{G({\mathbb{C}}_{p})}(U_{p})$ is an open subgroup of
$G({\mathbb{C}}_{p})$, and the map $\exp_{G({\mathbb{C}}_{p})}$ induces an
isomprphism between $U_{p}$ and $\exp_{G({\mathbb{C}}_{p})}(U_{p})\subset
G({\mathbb{C}}_{p})_{f}$, whose inverse is the restriction of
$\log_{G({\mathbb{C}}_{p})}$ to $\exp_{G({\mathbb{C}}_{p})}(G_{p})$. The map
$\exp_{G({\mathbb{C}}_{p})}$ is called the $p$-adic exponential map of $G$.
For an element $u$ in $U_{p}$ and for a subset $S$ in ${\mathbb{C}}_{p}$,
define $S\cdot u$ by the set $\\{su;s\in S\\}$. The following lemma plays a
central role in the proof of the main theorem which allows us to construct a
relevant commutative algebraic subgroup of $G$ defined over
$\overline{\mathbb{Q}}$.

######  Lemma 2.1.

Let $G$ be an algebraic group defined over $\overline{\mathbb{Q}}$. Let $u$ be
an element in $U_{p}$ with $\exp_{G({\mathbb{C}}_{p})}(u)\in
G(\overline{\mathbb{Q}})$. Denote by $G_{u}$ the Zariski closure of the set
$\exp_{G({\mathbb{C}}_{p})}(({\mathbb{C}}_{p}\cdot u)\cap U_{p})$ in
$G({\mathbb{C}}_{p})$. Then $G_{u}$ is a commutative algebraic subgroup of $G$
defined over $\overline{\mathbb{Q}}$ and $\exp_{G({\mathbb{C}}_{p})}(u)\in
G_{u}(\overline{\mathbb{Q}})$.

###### Proof.

Put $\Phi_{u}:=\exp_{G({\mathbb{C}}_{p})}(({\mathbb{C}}_{p}\cdot u)\cap
U_{p})$. We first see that any two elements $x,y\in{\mathbb{C}}_{p}\cdot u$
commute, and therefore

$\exp_{G({\mathbb{C}}_{p})}(x)\exp_{G({\mathbb{C}}_{p})}(y)=\exp_{G({\mathbb{C}}_{p})}(x+y)=\exp_{G({\mathbb{C}}_{p})}(y+x)=\exp_{G({\mathbb{C}}_{p})}(y)\exp_{G({\mathbb{C}}_{p})}(x),$

by Campbell-Hausdorff formula (see [5, Section 16]). In particular, this shows
that $\Phi_{u}$ is a commutative subgroup of $G({\mathbb{C}}_{p})$. It follows
from [4, Lemma 1.40] that the Zariski closure $G_{u}$ of $\Phi_{u}$ is an
algebraic subgroup of $G({\mathbb{C}}_{p})$. We now consider the commutator
morphism $f:G_{u}\times G_{u}\rightarrow G_{u}$ given by $f(x,y)=[x,y]$. Then
$f$ is an algebraic morphism and trivial on $\Phi_{u}\times\Phi_{u}$. This
gives $\Phi_{u}\times\Phi_{u}$ is contained in the algebraic set $f^{-1}(e)$,
where $e$ denotes the identity element of $G$. But $\Phi_{u}\times\Phi_{u}$ is
Zariski-dense in $G_{u}\times G_{u}$, this implies that the morphism $f$ must
be trivial on $G_{u}\times G_{u}$. In other words, $G_{u}$ is commutative. It
remains to show $G_{u}$ is defined over $\overline{\mathbb{Q}}$. If $u=0$ then
$G_{u}$ is the trivial group, and hence clearly defined over
$\overline{\mathbb{Q}}$. Assume $u\neq 0$, let $\Gamma_{u}$ denote the image
of the set $(\mathbb{Q}\cdot u)\cap U_{p}$ under the map
$\exp_{G({\mathbb{C}}_{p})}$. As above, we also have $\Gamma_{u}$ is a
subgroup of $G({\mathbb{C}}_{p})$, and therefore by [4, Lemma 1.40] again the
Zariski closure $G^{\prime}_{u}$ of $\Gamma_{u}$ in $G({\mathbb{C}}_{p})$ is
an algebraic group. On the other hand, since $\exp_{G({\mathbb{C}}_{p})}(u)\in
G(\overline{\mathbb{Q}})$, it follows that $\exp_{G({\mathbb{C}}_{p})}(u)\in
G_{u}(\overline{\mathbb{Q}})$. Note that
$\exp_{G({\mathbb{C}}_{p})}(mu)=\exp_{G({\mathbb{C}}_{p})}(u)^{m}$ for all
$m\in\mathbb{Z}$, and this leads to $\Gamma_{u}$ is contained in
$G_{u}(\overline{\mathbb{Q}})$. Hence, the algebraic group $G^{\prime}_{u}$
must be defined over $\overline{\mathbb{Q}}$. Moreover, the group $\Gamma_{u}$
is infinite (since $u\neq 0$), thus the algebraic group $G_{u}^{\prime}$ has
positive dimension.

Next, one has $\exp_{G({\mathbb{C}}_{p})}(u)\in
G^{\prime}_{u}(\overline{\mathbb{Q}})\cap G({\mathbb{C}}_{p})_{f}\subseteq
G^{\prime}_{u}({\mathbb{C}}_{p})_{f}$. This gives

$u=\log_{G({\mathbb{C}}_{p})}\big{(}\exp_{G({\mathbb{C}}_{p})}(u)\big{)}=\log_{G^{\prime}_{u}({\mathbb{C}}_{p})}\big{(}\exp_{G({\mathbb{C}}_{p})}(u)\big{)}\in\mathrm{Lie}(G^{\prime}_{u}({\mathbb{C}}_{p})).$

Since $\mathrm{Lie}(G^{\prime}_{u}({\mathbb{C}}_{p}))$ is a vector space over
${\mathbb{C}}_{p}$, it follows that the line ${\mathbb{C}}_{p}\cdot u$ is also
contained in $\mathrm{Lie}(G^{\prime}_{u}({\mathbb{C}}_{p}))$. From this, we
get

$\Phi_{u}\subseteq\exp_{G({\mathbb{C}}_{p})}\big{(}\mathrm{Lie}(G^{\prime}_{u}({\mathbb{C}}_{p}))\cap
U_{p}\big{)}=\exp_{G^{\prime}_{u}({\mathbb{C}}_{p})}\big{(}\mathrm{Lie}(G^{\prime}_{u}({\mathbb{C}}_{p}))\cap
U_{p}\big{)}\subseteq G^{\prime}_{u}({\mathbb{C}}_{p}).$

This allows us to deduce that $G_{u}({\mathbb{C}}_{p})\subset
G^{\prime}_{u}({\mathbb{C}}_{p})$. But, by definition,
$G^{\prime}_{u}({\mathbb{C}}_{p})$ is obviously contained in
$G_{u}({\mathbb{C}}_{p})$. Therefore, we conclude that $G_{u}=G^{\prime}_{u}$,
and this completes the proof of the lemma.

##  3\. Proof of the main theorem

We are now ready to prove the main theorem. Let
$\exp_{G({\mathbb{C}}_{p})}:U_{p}\rightarrow G({\mathbb{C}}_{p})$ be the
$p$-adic exponential map of $G$ defined as in Section 2. By definition of the
set $G({\mathbb{C}}_{p})_{f}$, there is a positive integer $k$ such that
$\gamma^{k}\in U_{p}$. Let $u$ be the point $\log_{G(C_{p})}(\gamma^{k})$. We
have

$\exp_{G({\mathbb{C}}_{p})}(u)=\exp_{G({\mathbb{C}}_{p})}\big{(}\log_{G(C_{p})}(\gamma^{k})\big{)}=\gamma^{k}\in
G(\overline{\mathbb{Q}}).$

Using Lemma 2.1, there exists a commutative algebraic subgroup $G_{u}$ of $G$
defined over $\overline{\mathbb{Q}}$ of positive dimension with
$\exp_{G({\mathbb{C}}_{p})}(u)\in G_{u}(\overline{\mathbb{Q}})$. This also
gives $\exp_{G({\mathbb{C}}_{p})}(u)\in G_{u}({\mathbb{C}}_{p})_{f}$, and
hence

$\log_{G({\mathbb{C}}_{p})}(\gamma^{k})=u\in\mathrm{Lie}(G_{u}({\mathbb{C}}_{p}))\subseteq\mathrm{Lie}(G({\mathbb{C}}_{p})).$

On the other hand, it follows from [2, Chapter III, 7.6] that

$\log_{G({\mathbb{C}}_{p})}(\gamma^{k})=k\log_{G({\mathbb{C}}_{p})}(\gamma)\in
V_{\gamma}\otimes_{\overline{\mathbb{Q}}}{\mathbb{C}}_{p}.$

In particular, this implies that $V_{\gamma}\subseteq\mathrm{Lie}(G_{u})$.
This enables us to apply the $p$-adic analytic subgroup theorem ([3, Theorem
2.2]) to the commutative algebraic group $G_{u}$, the algebraic point
$\gamma^{k}$ and the $\overline{\mathbb{Q}}$-vector space $V_{\gamma}$ to
obtain an algebraic (commutative) subgroup $H$ of $G$ defined over
$\overline{\mathbb{Q}}$ of positive dimension such that
$\mathrm{Lie}(H)\subseteq V_{\gamma}$ and $\gamma^{k}\in
H(\overline{\mathbb{Q}})$. Let $H^{0}$ be the connected component of $H$, then
$\mathrm{Lie}(H^{0})=\mathrm{Lie}(H)$ and there is a positive integer $m$ for
which $(\gamma^{k})^{m}\in H^{0}(\overline{\mathbb{Q}})$. Put $n=km$ and take
the algebraic subgroup $H_{\gamma}$ of $G_{u}$ with
$H_{\gamma}(\overline{\mathbb{Q}})=\\{\alpha\in
G_{u}(\overline{\mathbb{Q}}):\alpha^{n}\in H^{0}(\overline{\mathbb{Q}})\\}$.
Then $\gamma\in H_{\gamma}$ and $\mathrm{Lie}(H_{\gamma})\subseteq
V_{\gamma}$. This leads to

$\log_{G({\mathbb{C}}_{p})}(\gamma)=\log_{H_{\gamma}({\mathbb{C}}_{p})}(\gamma)\in\mathrm{Lie}(H_{\gamma}({\mathbb{C}}_{p}))=\mathrm{Lie}(H_{\gamma})\otimes_{\overline{\mathbb{Q}}}{\mathbb{C}}_{p}.$

By the minimality of the $\overline{\mathbb{Q}}$-vector space $V_{\gamma}$, we
conclude that $\mathrm{Lie}(H_{\gamma})=V_{\gamma}$. The theorem is therefore
proved. ∎ Remark. By the same method given in the proof of [11, Theorem 1.5]
(which is deduced as a consequence of the non-commutative analytic subgroup
theorem), we just use the $p$-adic exponential map instead of the complex one
to also obtain the following result as a consequence of the $p$-adic non-
commutative analytic subgroup theorem.

######  Theorem 3.1.

Let $G$ be an algebraic group defined over $\overline{\mathbb{Q}}$ and
$\exp_{G({\mathbb{C}}_{p})}:U_{p}\rightarrow G({\mathbb{C}}_{p})$ be the
$p$-adic exponential map of $G$. Let $u$ be a non-zero element of
$U_{p}\cap\mathrm{Lie}(G)$ such that $\exp_{G({\mathbb{C}}_{p})}(u)$ lies in
$G(\overline{\mathbb{Q}})$.
(1) The element $\exp_{G({\mathbb{C}}_{p})}(u)$ is contained in a subgroup
isomorphic to the additive group $\mathbb{G}_{a}$.
(2) Assume $G$ is an affine algebraic group. The element $u$ is a nilpotent
element of $\mathrm{Lie}(G({\mathbb{C}}_{p}))$, and hence
$\exp_{G({\mathbb{C}}_{p})}(u)$ is a unipotent element of
$G({\mathbb{C}}_{p})$.

## References

  * [1] A. Baker and G. Wüstholz, _Logarithmic Forms and Diophantine Geometry_ , New Mathematical Monographs, Cambridge University Press, 2007.
  * [2] N. Bourbaki, _Elements of Mathematics. Lie groups and Lie algebras. Part I: Chapters 1-3. English translantion.,_ Actualities scientifiques et industrielles, Herman. Adiwes International Series in Mathematics. Paris: Hermann, Publishers in Arts and Science; Reading, Mass.: Addison-Wesley Publishing Company. XVII, 1975.
  * [3] C. Fuchs and D. H. Pham, _The $p$-adic analytic subgroup theorem revisited_, $p$-Adic Number, Ultrametric Analysis and Applications 7 (2015), 143–156.
  * [4] J. Milne, _Algebraic Groups_ , Cambride Studies in Advanced Mathematics, Cambride University Press, 2017.
  * [5] P. Schneider, _$p$ -adic Lie groups_, Grundlehren der mathematischen Wissenschaften, Springer Verlag, 2011.
  * [6] P. Tretkoff, _Periods and Special Functions in Transcendence_ , Advanced Textbooks in Mathematics, World Scientific, 2017.
  * [7] G. Wüstholz, _Some remarks on a conjecture of Waldschmidt_ , Diophantine approximations and transcendental numbers, Birkhäuser Boston, Boston, MA 1983, 329–336.
  * [8] G. Wüstholz, _Multiplicity estimates on group varieties_ , Ann. of Math. 129 (1989), 471–500.
  * [9] G. Wüstholz, _Algebraische Punkte auf Analytischen Untergruppen algebraischer Gruppen_ , Ann. of Math. 129 (1989), 501–517.
  * [10] G. Wüstholz, _Elliptic and abelian period spaces_ , Acta Arith. 198 (2021), 329–357.
  * [11] A. Yafaev, _Non-commutative analytic subgroup theorem_ , J. Number Theory (2021), doi: https://doi.org/10.1016/j.jnt.2021.01.025.

Duc Hiep Pham
University of Education
Vietnam National University, Hanoi
144 Xuan Thuy, Cau Giay, Hanoi
Vietnam
Email<EMAIL_ADDRESS>