Energy-consistent integration of mechanical systems
based on Livens principle

Philipp L. Kinon, Peter Betsch

Institute of Mechanics
---
Karlsruhe Institute of Technology (KIT)
Otto-Ammann-Platz 9, 76131 Karlsruhe, Germany
[philipp.kinon<EMAIL_ADDRESS>

> ## ABSTRACT
>
> In this work we make us of Livens principle (sometimes also referred to as
> Hamilton-Pontryagin principle) in order to obtain a novel structure-
> preserving integrator for mechanical systems. In contrast to the canonical
> Hamiltonian equations of motion, the Euler-Lagrange equations pertaining to
> Livens principle circumvent the need to invert the mass matrix. This is an
> essential advantage with respect to singular mass matrices, which can yield
> severe difficulties for the modelling and simulation of multibody systems.
> Moreover, Livens principle unifies both Lagrangian and Hamiltonian
> viewpoints on mechanics. Additionally, the present framework avoids the need
> to set up the system’s Hamiltonian. The novel scheme algorithmically
> conserves a general energy function and aims at the preservation of momentum
> maps corresponding to symmetries of the system. We present an extension to
> mechanical systems subject to holonomic constraints. The performance of the
> newly devised method is studied in representative examples.
>
> Keywords: Livens principle, Structure-preserving integration, Singular mass
> matrix, Energy-momentum method, Discrete gradients, Constrained mechanical
> systems.

##  1 INTRODUCTION

Over time, two approaches to the description of dynamical systems have been
developed: The well-known Lagrangian and Hamiltonian formalisms both consider
descriptive energetic scalars and deploy certain operations on them to
generate the system’s equations of motion. However, another formulation, which
unifies both frameworks by means of independent position, velocity and
momentum quantities has been proposed by Livens [1]. This Livens principle has
been recently taken up (cf. [2, 3]) under the name of Hamilton-Pontryagin
principle due to its close relation to the Pontryagin principle from the field
of optimal control [4]. Livens principle allows for an advantageous universal
description due to its mixed character.

Additionally, when simulating multibody system dynamics, the mathematical
properties of the formulation heavily depend on the particular choice of
coordinates. In certain frameworks, the mass matrix can even become singular
or does depend on the coordinates themselves (cf. [5, 6]). These cases pose
some difficulties for Hamiltonian methods, since the inversion of the mass
matrix might not be admissible.

Thus, in the present work, we propose the application of Livens principle to
multibody system dynamics which circumvents the problem of singular mass-
matrices. This extends previous works [7, 8], in which a novel variational
principle has been derived based on Livens principle, with respect to more
general mass matrices.

In recent times, the focus of scientific research has become to derive
numerical integration methods, which are capable of solving the equations of
motion approximately. Therefore, the class of structure-preserving integrators
seeks to inherit the conservation principles of dynamical systems in a
discrete sense (cf. monographs such as Hairer et al. [9]). In the field of
mechanics, structure-preserving integration schemes can be mainly divided into
two different groups: symplectic methods, of which variational integrators are
an important subclass, and energy-momentum integrators.

In the present contribution we focus on the approach of energy-momentum
methods: Based on the direct discretization of the Livens equations of motion
in conjunction with Gonzalez discrete derivatives [10], we obtain an energy-
consistent time-stepping scheme. This method is able to take into account both
singular and configuration-dependent mass matrices. It achieves discrete
energy preservation and satisfies holonomic constraint functions exactly.

The outline of the present work is as follows: In Sec. 2 the fundamentals of
Hamiltonian dynamics are briefly recapitulated. Thereafter, we introduce
Livens principle in Sec. 3 and extend it to holonomically constrained systems
in Sec. 4. In Sec. 5 we discretize the deduced equations of motion to obtain
an energy-consistent time-stepping scheme, which is used in Sec. 6 for the
simulation of representative example problems. Sec. 7 concludes the findings
and gives a short outlook for future research.

##  2 HAMILTONIAN DYNAMICS

Consider the motion of a dynamical system with $d$ degrees of freedom with
positions $\mbox{$\mathbf{q}$}\in Q$ and velocities
$\dot{\mbox{$\mathbf{q}$}}\in T_{\mathbf{q}}Q$, where $Q\subset\mathbb{R}^{d}$
denotes the configuration space and $T_{\mathbf{q}}Q\subset\mathbb{R}^{d}$ the
tangent space to $Q$ at $\mathbf{q}$.

######  Definition 2.1.

The system’s Lagrangian $L\colon TQ\to\mathbb{R}$ is given by

$\displaystyle
L(\mbox{$\mathbf{q}$},\dot{\mbox{$\mathbf{q}$}})=T(\mbox{$\mathbf{q}$},\dot{\mbox{$\mathbf{q}$}})-V(\mbox{$\mathbf{q}$})=\frac{1}{2}\dot{\mbox{$\mathbf{q}$}}\cdot\mbox{$\mathbf{M}$}(\mbox{$\mathbf{q}$})\dot{\mbox{$\mathbf{q}$}}-V(\mbox{$\mathbf{q}$}),$
(1)

where $T:TQ\rightarrow\mathbb{R}$ is the kinetic energy,
$\mbox{$\mathbf{M}$}(\mbox{$\mathbf{q}$})\in\mathbb{R}^{d\times d}$ denotes
the symmetric and positive-definite mass matrix and $V:Q\rightarrow\mathbb{R}$
is a potential function.

The total energy of the system is given by the sum of kinetic and potential
energy, i.e.

$\displaystyle
E_{\mathrm{tot}}(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$})=T(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$})+V(\mbox{$\mathbf{q}$}).$
(2)

######  Definition 2.2.

The Hamiltonian $H:T^{*}Q\rightarrow\mathbb{R}$ is obtained by employing a
Legendre transformation
$\mathbb{F}_{L}:(\mbox{$\mathbf{q}$},\dot{\mbox{$\mathbf{q}$}})\mapsto(\mbox{$\mathbf{q}$},\mbox{$\mathbf{p}$})$,
where the conjugate momenta are defined throught the fiber derivative

$\displaystyle\mbox{$\mathbf{p}$}:=D_{2}L(\mbox{$\mathbf{q}$},\dot{\mbox{$\mathbf{q}$}})$
(3)

such that

$\displaystyle
H(\mbox{$\mathbf{q}$},\mbox{$\mathbf{p}$})=\mbox{$\mathbf{p}$}\cdot\dot{\mbox{$\mathbf{q}$}}(\mbox{$\mathbf{q}$},\mbox{$\mathbf{p}$})-L(\mbox{$\mathbf{q}$},\dot{\mbox{$\mathbf{q}$}}(\mbox{$\mathbf{q}$},\mbox{$\mathbf{p}$})).$
(4)

Given a Lagrangian (1), the Hamiltonian reads

$\displaystyle
H(\mbox{$\mathbf{q}$},\mbox{$\mathbf{p}$})=T(\mbox{$\mathbf{q}$},\mbox{$\mathbf{p}$})+V(\mbox{$\mathbf{q}$})=\frac{1}{2}\mbox{$\mathbf{p}$}\cdot\mbox{$\mathbf{M}$}(\mbox{$\mathbf{q}$})^{-1}\mbox{$\mathbf{p}$}+V(\mbox{$\mathbf{q}$}).$
(5)

Eventually, the well-known Hamiltonian equations of motion appear in their
canonical form

$\displaystyle\dot{\mbox{$\mathbf{q}$}}$
$\displaystyle=\operatorname{D}_{2}H(\mbox{$\mathbf{q}$},\mbox{$\mathbf{p}$})$
(6a) $\displaystyle\dot{\mbox{$\mathbf{p}$}}$
$\displaystyle=-\operatorname{D}_{1}H(\mbox{$\mathbf{q}$},\mbox{$\mathbf{p}$})$
(6b)

which can be specified for the systems with Lagrangian (1) as

$\displaystyle\dot{\mbox{$\mathbf{q}$}}$
$\displaystyle=\mbox{$\mathbf{M}$}(\mbox{$\mathbf{q}$})^{-1}\mbox{$\mathbf{p}$},$
(7a) $\displaystyle\dot{\mbox{$\mathbf{p}$}}$
$\displaystyle=-\operatorname{D}_{1}T(\mbox{$\mathbf{q}$},\mbox{$\mathbf{p}$})-\operatorname{D}V(\mbox{$\mathbf{q}$}),$
(7b)

Note that this approach requires to set up the Hamiltonian which in turn
requires the inversion of the mass matrix. Thus, we now want to introduce a
different and more general framework.

##  3 LIVENS PRINCIPLE

From Hamilton’s principle of least action

$\displaystyle\delta\int_{0}^{T}L(\mbox{$\mathbf{q}$},\dot{\mbox{$\mathbf{q}$}})\mathrm{d}t=0$
(8)

one can proceed by allowing the velocities to be independent variables
$\mbox{$\mathbf{v}$}\in T_{\mathbf{q}}Q$. The kinematic relation
$\dot{\mbox{$\mathbf{q}$}}=\mbox{$\mathbf{v}$}$ can be enforced by means of
Lagrange multipliers $\mbox{$\mathbf{p}$}\in T^{*}_{\mathbf{q}}Q$.

######  Definition 3.1.

The augmented action integral for Livens principle reads

$\displaystyle
S(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$})=\int_{0}^{T}\left[L(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$})+\mbox{$\mathbf{p}$}\cdot(\dot{\mbox{$\mathbf{q}$}}-\mbox{$\mathbf{v}$})\right]\,\mathrm{d}t$
(9)

for the simulation time $t\in[0,T]$.

By stating the stationary condition

$\displaystyle\delta
S(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$})=0$ (10)

and executing the variations with respect to every independent variable, one
obtains the equations of motion

$\displaystyle\dot{\mbox{$\mathbf{q}$}}$ $\displaystyle=\mbox{$\mathbf{v}$},$
(11a) $\displaystyle\dot{\mbox{$\mathbf{p}$}}$
$\displaystyle=\operatorname{D}_{1}L(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$}),$
(11b) $\mathbf{p}$
$\displaystyle=\operatorname{D}_{2}L(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$}),$
(11c)

where the arbitraryness of variations and standard endpoint conditions
$\delta\mbox{$\mathbf{q}$}(t=0)=\delta\mbox{$\mathbf{q}$}(t=T)=\mbox{$\mathbf{0}$}$
have been taken into account.

With regard to (11c) the Lagrange multiplier $\mathbf{p}$ can be identified as
the conjugate momentum. Thus, Livens principle automatically accounts for the
Legendre transformation (4), whereas within the framework of Hamiltonian
dynamics momentum variables have to be defined a priori using the fiber
derivative. In analogy to (4), a generalized energy function can be introduced
as follows.

######  Definition 3.2.

Given a Lagrangian (1) the quantity

$\displaystyle
E(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$})=\mbox{$\mathbf{p}$}\cdot\mbox{$\mathbf{v}$}-L(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$})$
(12)

is called a generalized energy function.

######  Remark 3.3 (Alternative form of Livens principle).

Using the generalized energy function, (9) can be recast as

$\displaystyle
S(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$})=\int_{0}^{T}\left[\mbox{$\mathbf{p}$}\cdot\dot{\mbox{$\mathbf{q}$}}-E(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$})\right]\mathrm{d}t,$
(13)

such that the corresponding Euler-Lagrange equations (11) emerge as

$\displaystyle\dot{\mbox{$\mathbf{q}$}}$
$\displaystyle=\operatorname{D}_{3}E(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$}),$
(14a) $\displaystyle\dot{\mbox{$\mathbf{p}$}}$
$\displaystyle=-\operatorname{D}_{1}E(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$}),$
(14b) $\mathbf{0}$
$\displaystyle=\operatorname{D}_{2}E(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$}).$
(14c)

######  Proposition 3.4.

The generalized energy (12) is a conserved quantity along solutions of (11).

###### Proof.

We compute

$\displaystyle\frac{\,\mathrm{d}}{\mathrm{d}t}E(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$})$
$\displaystyle=\operatorname{D}_{1}E(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$})\cdot\dot{\mbox{$\mathbf{q}$}}+\operatorname{D}_{2}E(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$})\cdot\dot{\mbox{$\mathbf{v}$}}+\operatorname{D}_{3}E(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$})\cdot\dot{\mbox{$\mathbf{p}$}}$
(15)
$\displaystyle=-\operatorname{D}_{1}L(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$})\cdot\mbox{$\mathbf{v}$}+(\mbox{$\mathbf{p}$}-\operatorname{D}_{2}L(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$}))\cdot\dot{\mbox{$\mathbf{v}$}}+\mbox{$\mathbf{v}$}\cdot\operatorname{D}_{1}L(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$}),$
(16)

where the partial derivatives of $E$ have been used and (11) has been
inserted. Thus, the second term vanishes while the first and third term cancel
each other out, such that $\,\mathrm{d}E/\mathrm{d}t=0$. ∎

######  Remark 3.5 (Equivalence with Lagrangian dynamics).

Note that after reinserting (11c) into (11b) and making use of (11a), Livens
principle traces back to the Lagrangian equations of the second kind

$\displaystyle\frac{\,\mathrm{d}}{\mathrm{d}t}\left(\operatorname{D}_{2}L(\mbox{$\mathbf{q}$},\dot{\mbox{$\mathbf{q}$}})\right)-\operatorname{D}_{1}L(\mbox{$\mathbf{q}$},\dot{\mbox{$\mathbf{q}$}})=\mbox{$\mathbf{0}$}.$
(17)

######  Remark 3.6 (Equivalence with Hamiltonian dynamics).

For mechanical systems with Lagrangian (1), relation (11c) yields
$\mbox{$\mathbf{p}$}=M(\mbox{$\mathbf{q}$})\mbox{$\mathbf{v}$}$, so that, if
$M$ is nonsingular, (11a) and (11b) directly lead to the canonical Hamiltonian
equations of motion (7). Moreover, in this case, the generalized energy
function can be identified as the Hamiltonian, such that
$E(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$}(\mbox{$\mathbf{q}$},\mbox{$\mathbf{p}$}),\mbox{$\mathbf{p}$})=H(\mbox{$\mathbf{q}$},\mbox{$\mathbf{p}$})$.

Initially termed Livens principle (cf. Sec. 26.2 in Pars [11]) after G.H.
Livens, who proposed this functional for the first time (cf. Livens [1]),
Marsden and co-workers (e.g. [2]) coined the name Hamilton-Pontryagin
principle for this functional due to its close relation to the classical
Pontryagin principle from the field of optimal control. Due to its mixed
character with three independent fields ($q,v,p$), it resembles the Hu-Washizu
principle from the area of elasticity theory. More recently, preliminary works
[7, 8] by the authors of this constribution have taken up Livens principle in
order to obtain a novel variational principle, the GGL principle which extends
Livens principle to holonomically constrained systems. In these works, new
structure-preserving integrators for systems with constant mass matrix have
been developed.

##  4 HOLONOMICALLY CONSTRAINED SYSTEMS

Assume that the coordinates $\mathbf{q}$ are redundant due to the presence of
$m$ independent scleronomic, holonomic constraints. Arranging the constraint
functions $g_{k}\colon\mathbb{R}^{d}\rightarrow\mathbb{R}$ ($k=1,\ldots,m$) in
a column vector $\mbox{$\mathbf{g}$}\in\mathbb{R}^{m}$, the constraints are
given by

$\displaystyle\mbox{$\mathbf{g}$}(\mbox{$\mathbf{q}$})=\mbox{$\mathbf{0}$}\
\text{.}$ (18)

The constraint functions are assumed to be independent, so that the constraint
Jacobian $\operatorname{D}\mbox{$\mathbf{g}$}(\mbox{$\mathbf{q}$})$ has rank
$m$. Note that the constraints give rise to the configuration manifold

$Q=\left\\{\mbox{$\mathbf{q}$}\in\mathbb{R}^{d}\;|\;\mbox{$\mathbf{g}$}(\mbox{$\mathbf{q}$})=\mbox{$\mathbf{0}$}\right\\}\
\text{.}$ (19)

Correspondingly, Livens principle can be augmented to account for the
constraints (18). Accordingly, the functional assumes the form

$\hat{S}(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$},{\bm{\lambda}})={S}(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$})+\int\nolimits_{0}^{T}{\bm{\lambda}}\cdot\mbox{$\mathbf{g}$}(\mbox{$\mathbf{q}$})\,\mathrm{d}t\
\text{,}$ (20)

where ${\bm{\lambda}}\in\mathbb{R}^{m}$ represents Lagrange multipliers for
the enforcement of the constraints. Imposing the stationary conditions on
functional (20), the resulting Euler-Lagrange equations are obtained as

$\displaystyle\dot{\mbox{$\mathbf{q}$}}$ $\displaystyle=\mbox{$\mathbf{v}$}\
\text{,}$ (21a) $\displaystyle\dot{\mbox{$\mathbf{p}$}}$
$\displaystyle=\operatorname{D}_{1}L(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$})-\operatorname{D}\mbox{$\mathbf{g}$}(\mbox{$\mathbf{q}$})^{\mathrm{T}}{\bm{\lambda}}\
\text{,}$ (21b) $\mathbf{p}$
$\displaystyle=\operatorname{D}_{2}L(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$})\
\text{,}$ (21c) $\mathbf{0}$
$\displaystyle=\mbox{$\mathbf{g}$}(\mbox{$\mathbf{q}$})\ \text{.}$ (21d)

######  Remark 4.1.

Due to the constraints, coordinates become redundant. Hence, the mass matrix
$\mbox{$\mathbf{M}$}(\mbox{$\mathbf{q}$})$ contained in the Lagrangian can be
singular such that the mass matrix is positive semidefinite in general.

Note that (18) at every time induces the secondary constraints or hidden
velocity constraints
$\,\mathrm{d}\mbox{$\mathbf{g}$}/\mathrm{d}t=\mbox{$\mathbf{0}$}$. Taking into
account (21a), one obtains

$\displaystyle\operatorname{D}\mbox{$\mathbf{g}$}(\mbox{$\mathbf{q}$})\mbox{$\mathbf{v}$}=\mbox{$\mathbf{0}$}$
(22)

It is well-known that the above differential algebraic equations (DAEs) (21)
have differentiation index 3 (see, for example, [12]), which can lead to
numerical instabilities. An index reduction using the classical Gear-Gupta-
Leimkuhler stabilization ([13]) or an expansion of (20) to account also for
the hidden constraints (cf. GGL principle in [7, 8]) can circumvent the
arising problems.

We now want to demonstrate that the conservation principle also holds in the
constrained setting.

######  Proposition 4.2.

The generalized energy (12) is a conserved quantity along solutions of (21).

###### Proof.

As above, we compute

$\displaystyle\frac{\,\mathrm{d}}{\mathrm{d}t}E(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$})$
$\displaystyle=\operatorname{D}_{1}E(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$})\cdot\dot{\mbox{$\mathbf{q}$}}+\operatorname{D}_{2}E(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$})\cdot\dot{\mbox{$\mathbf{v}$}}+\operatorname{D}_{3}E(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$})\cdot\dot{\mbox{$\mathbf{p}$}}$
(23)
$\displaystyle=-\operatorname{D}_{1}L(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$})\cdot\mbox{$\mathbf{v}$}+(\mbox{$\mathbf{p}$}-\operatorname{D}_{2}L(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$}))\cdot\dot{\mbox{$\mathbf{v}$}}+\mbox{$\mathbf{v}$}\cdot(\operatorname{D}_{1}L(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$})-\mbox{$\mathbf{G}$}(\mbox{$\mathbf{q}$})^{\mathrm{T}}{\bm{\lambda}}),$
(24)

where the partial derivatives of $E$ have been used and (11) has been
inserted. Thus, the second term vanishes while the first and third term cancel
each other out. Eventually,

$\displaystyle\frac{\,\mathrm{d}}{\mathrm{d}t}E(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$},\mbox{$\mathbf{p}$})$
$\displaystyle={\bm{\lambda}}\cdot(\operatorname{D}\mbox{$\mathbf{g}$}(\mbox{$\mathbf{q}$})\mbox{$\mathbf{v}$})$
(25)

remains. Having in mind that (22) holds, also this last part vanished, i.e.
$\,\mathrm{d}E/\mathrm{d}t=0$. ∎

##  5 STRUCTURE-PRESERVING DISCRETIZATION

In this work, a novel integration method is proposed, which conserves first
integrals of the equations of motion at hand, e.g. the generalized energy
function $E$. This scheme results from a direct discretization of the Euler-
Lagrange equations (21) emanating from constrained Livens principle.
Particularly, discrete derivatives in the sense of Gonzalez [10] are taken
into account and are defined as follows.

######  Definition 5.1 (Gonzalez discrete gradient).

For a given function $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$ the discrete
derivative
$\overline{\mathsf{D}}:\mathbb{R}^{m}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$
according to [10] is defined by

$\displaystyle\overline{\mathsf{D}}f(\mbox{$\mathbf{x}$},\mbox{$\mathbf{y}$})=\operatorname{D}f(\mbox{$\mathbf{z}$})+\frac{f(\mbox{$\mathbf{y}$})-f(\mbox{$\mathbf{x}$})-\operatorname{D}f(\mbox{$\mathbf{z}$})\cdot\mbox{$\mathbf{v}$}}{||\mbox{$\mathbf{v}$}||^{2}}\mbox{$\mathbf{v}$}$
(26)

with
$\mbox{$\mathbf{z}$}=\frac{1}{2}(\mbox{$\mathbf{x}$}+\mbox{$\mathbf{y}$})$,
$\mbox{$\mathbf{v}$}=\mbox{$\mathbf{y}$}-\mbox{$\mathbf{x}$}$ and
$||\bullet||$ denoting the standard Euclidean norm on $\mathbb{R}^{m}$.

Those discrete derivatives represent second-order approximation to the exact
gradients such that discrete derivatives satisfy the directionality condition

$\displaystyle\overline{\mathsf{D}}f(\mbox{$\mathbf{x}$},\mbox{$\mathbf{y}$})\cdot\mbox{$\mathbf{v}$}=f(\mbox{$\mathbf{y}$})-f(\mbox{$\mathbf{x}$})$
(27)

as well as the consistency condition
$\overline{\mathsf{D}}f(\mbox{$\mathbf{x}$},\mbox{$\mathbf{y}$})=\operatorname{D}f(\mbox{$\mathbf{z}$})+O(||\mbox{$\mathbf{v}$}||)$.
Eventually, we propose the energy-consistent integrator, governed by

$\displaystyle\mbox{$\mathbf{q}$}^{n+1}-\mbox{$\mathbf{q}$}^{n}$
$\displaystyle=h\,\mbox{$\mathbf{v}$}^{n+1/2},$ (28a)
$\displaystyle\mbox{$\mathbf{p}$}^{n+1}-\mbox{$\mathbf{p}$}^{n}$
$\displaystyle=h\overline{\mathsf{D}}_{1}L(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$})-h\,\sum_{k=1}^{m}\lambda_{k}\overline{\mathsf{D}}g_{k}(\mbox{$\mathbf{q}$}),$
(28b) $\displaystyle\mbox{$\mathbf{p}$}^{n+1/2}$
$\displaystyle=\overline{\mathsf{D}}_{2}L(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$}),$
(28c) $\displaystyle\mbox{$\mathbf{g}$}(\mbox{$\mathbf{q}$}^{n+1})$
$\displaystyle=\mbox{$\mathbf{0}$},$ (28d)

where $h$ denotes the time step size,
$\mbox{$\mathbf{v}$}^{n+1/2}=\frac{1}{2}(\mbox{$\mathbf{v}$}^{n}+\mbox{$\mathbf{v}$}^{n+1})$
and
$\mbox{$\mathbf{p}$}^{n+1/2}=\frac{1}{2}(\mbox{$\mathbf{p}$}^{n}+\mbox{$\mathbf{p}$}^{n+1})$.
Moreover, $\overline{\mathsf{D}}_{i}$ denotes the partitioned discrete
derivative with respect to the $i$-th argument. Directionality condition (27)
extends to partitioned versions such that

$\displaystyle\overline{\mathsf{D}}_{1}T(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$})\cdot\Delta\mbox{$\mathbf{q}$}+\overline{\mathsf{D}}_{2}T(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$})\cdot\Delta\mbox{$\mathbf{v}$}=T(\mbox{$\mathbf{q}$}^{n+1},\mbox{$\mathbf{v}$}^{n+1})-T(\mbox{$\mathbf{q}$}^{n},\mbox{$\mathbf{v}$}^{n}),$
(29)

where
$\Delta\mbox{$\mathbf{q}$}=\mbox{$\mathbf{q}$}^{n+1}-\mbox{$\mathbf{q}$}^{n}$
and
$\Delta\mbox{$\mathbf{v}$}=\mbox{$\mathbf{v}$}^{n+1}-\mbox{$\mathbf{v}$}^{n}$.

######  Remark 5.2.

Assuming (1), the present method (28) can be specified as

$\displaystyle\mbox{$\mathbf{q}$}^{n+1}-\mbox{$\mathbf{q}$}^{n}$
$\displaystyle=h\,\mbox{$\mathbf{v}$}^{n+1/2},$ (30a)
$\displaystyle\mbox{$\mathbf{p}$}^{n+1}-\mbox{$\mathbf{p}$}^{n}$
$\displaystyle=h\,\overline{\mathsf{D}}_{1}T(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$})-h\,\overline{\mathsf{D}}V(\mbox{$\mathbf{q}$})-h\,\sum_{k=1}^{m}\lambda_{k}\overline{\mathsf{D}}g_{k}(\mbox{$\mathbf{q}$}),$
(30b) $\displaystyle\mbox{$\mathbf{p}$}^{n+1/2}$
$\displaystyle=\overline{\mathsf{D}}_{2}T(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$}),$
(30c) $\displaystyle\mbox{$\mathbf{g}$}(\mbox{$\mathbf{q}$}^{n+1})$
$\displaystyle=\mbox{$\mathbf{0}$}.$ (30d)

######  Definition 5.3.

A discrete version of the energy function (12) is given by

$\displaystyle E^{n}$
$\displaystyle=\mbox{$\mathbf{p}$}^{n}\cdot\mbox{$\mathbf{v}$}^{n}-\frac{1}{2}\mbox{$\mathbf{v}$}^{n}\cdot\mbox{$\mathbf{M}$}(\mbox{$\mathbf{q}$}^{n})\mbox{$\mathbf{v}$}^{n}+V(\mbox{$\mathbf{q}$}^{n}).$
(31)

######  Proposition 5.4.

The generalized energy (31) is conserved by time-stepping method (28) in a
discrete sense.

###### Proof.

Scalar multiplying (28b) with $\Delta\mbox{$\mathbf{q}$}$ and adding the
scalar product of (28c) with $\Delta\mbox{$\mathbf{v}$}$ yields

$\displaystyle(\mbox{$\mathbf{p}$}^{n+1}-\mbox{$\mathbf{p}$}^{n})\cdot(\mbox{$\mathbf{q}$}^{n+1}-\mbox{$\mathbf{q}$}^{n})+h\,\mbox{$\mathbf{p}$}^{n+1/2}\cdot(\mbox{$\mathbf{v}$}^{n+1}-\mbox{$\mathbf{v}$}^{n})$
$\displaystyle=$ (32) $\displaystyle
h\,\overline{\mathsf{D}}_{1}L(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$})\cdot(\mbox{$\mathbf{q}$}^{n+1}-\mbox{$\mathbf{q}$}^{n})-h\,\sum_{k=1}^{m}\lambda_{k}$
$\displaystyle\overline{\mathsf{D}}g_{k}(\mbox{$\mathbf{q}$})\cdot(\mbox{$\mathbf{q}$}^{n+1}-\mbox{$\mathbf{q}$}^{n})+\overline{\mathsf{D}}_{2}L(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$})\cdot(\mbox{$\mathbf{v}$}^{n+1}-\mbox{$\mathbf{v}$}^{n})$

Inserting (28a) and making use of the directionality property gives

$\displaystyle(\mbox{$\mathbf{p}$}^{n+1}-\mbox{$\mathbf{p}$}^{n})\cdot\mbox{$\mathbf{v}$}^{n+1/2}+\mbox{$\mathbf{p}$}^{n+1/2}\cdot(\mbox{$\mathbf{v}$}^{n+1}-\mbox{$\mathbf{v}$}^{n})$
$\displaystyle=L(\mbox{$\mathbf{q}$}^{n+1},\mbox{$\mathbf{v}$}^{n+1})-L(\mbox{$\mathbf{q}$}^{n},\mbox{$\mathbf{v}$}^{n})$
(33)
$\displaystyle\qquad\qquad-\sum_{k=1}^{m}\lambda_{k}\left(g_{k}(\mbox{$\mathbf{q}$}^{n+1})-g_{k}(\mbox{$\mathbf{q}$}^{n})\right).$

After cancelling some terms on the left-hand side of the previous relation and
taking into account that the constraints are identically satisfied in each
time step (see (28d)), one eventually arrives at

$\displaystyle
E(\mbox{$\mathbf{q}$}^{n+1},\mbox{$\mathbf{v}$}^{n+1},\mbox{$\mathbf{p}$}^{n+1})=E(\mbox{$\mathbf{q}$}^{n},\mbox{$\mathbf{v}$}^{n},\mbox{$\mathbf{p}$}^{n})$
(34)

which concludes the proof. ∎

######  Proposition 5.5.

The discrete fiber derivative (30c) in general does not imply a pointwise
fulfillment of the continuous fiber derivative. That is,
$\mbox{$\mathbf{p}$}^{n}=\mbox{$\mathbf{M}$}(\mbox{$\mathbf{q}$}^{n})\mbox{$\mathbf{v}$}^{n}$
does not hold in general. However, the relation
$\mbox{$\mathbf{p}$}^{n}=\mbox{$\mathbf{M}$}\mbox{$\mathbf{v}$}^{n}$ is
satisfied for constant mass matrices.

###### Proof.

In the case of a constant mass matrix, the discrete derivative is given by
$\overline{\mathsf{D}}_{2}T(\mbox{$\mathbf{q}$},\mbox{$\mathbf{v}$})=\mbox{$\mathbf{M}$}\mbox{$\mathbf{v}$}^{n+1/2}$.
Accordingly, the discrete fiber derivative (30c) yields
$\mbox{$\mathbf{p}$}^{n+1/2}=\mbox{$\mathbf{M}$}\mbox{$\mathbf{v}$}^{n+1/2}$.
Consequently, provided that
$\mbox{$\mathbf{p}$}^{n}=\mbox{$\mathbf{M}$}\mbox{$\mathbf{v}$}^{n}$ is
satisfied,
$\mbox{$\mathbf{p}$}^{n+1}=\mbox{$\mathbf{M}$}\mbox{$\mathbf{v}$}^{n+1}$ is
fulfilled, too. ∎

######  Remark 5.6.

Using the proposed integrator (30), the discrete total energy

$\displaystyle
E_{\mathrm{tot}}^{n}=\frac{1}{2}\mbox{$\mathbf{v}$}^{n}\cdot\mbox{$\mathbf{M}$}(\mbox{$\mathbf{q}$}^{n})\mbox{$\mathbf{v}$}^{n}+V(\mbox{$\mathbf{q}$}^{n})$
(35)

is in general not identical to the discrete generalized energy function (31).
This is only the case, if
$\mbox{$\mathbf{p}$}^{n}=\mbox{$\mathbf{M}$}(\mbox{$\mathbf{q}$}^{n})\mbox{$\mathbf{v}$}^{n}$.

##  6 NUMERICAL EXAMPLE

In the following, the newly proposed scheme (30) is applied to some numerical
examples. Since this involves the solution of an implicit set of equations,
Newton’s method is used in every time step with a tolerance of
$\epsilon_{\mathrm{Newton}}$. The computations have been performed using the
code package available at [14], which can also be used for verification.

###  6.1 Mass-spring system with redundant coordinates

$k_{1}$
---
$m_{1}$
---
$w$
---
$k_{2}$
---
$m_{2}$
---
$x_{1}+l_{10}$
---
$q_{2}$
---
$x_{2}+l_{20}$
---
Figure 1: Spring-mass multibody system.

We take up an example with singular mass matrix from [5] (cf. Sec. 5, Example
3) as depicted in Fig. 1 with masses $m_{1}$ and $m_{2}$ and springs with
constants $k_{1}$ and $k_{2}$ and resting lengths $l_{10}$ and $l_{20}$,
respectively. By splitting up the system into two seperate subsystems it can
be considered a multibody system with two degrees of freedom. However, we use
two absolute coordinates ($x_{1}$ and $q_{2}$) and one relative coordinate
$x_{2}$ to define the systems kinematics, i.e. $d=3$, and

$\displaystyle\mbox{$\mathbf{q}$}=\begin{bmatrix}x_{1}\\\ q_{2}\\\
x_{2}\end{bmatrix}$ (36)

with corresponding velocities

$\displaystyle\mbox{$\mathbf{v}$}=\begin{bmatrix}v_{1}\\\ v_{2}\\\
v_{3}\end{bmatrix}=\begin{bmatrix}\dot{x}_{1}\\\ \dot{q}_{2}\\\
\dot{x}_{2}\end{bmatrix}.$ (37)

Consequently, the interconnection constraint $q_{2}=x_{1}+l_{10}+w$ arises,
which can be recast as

$\displaystyle
g(\mbox{$\mathbf{q}$})=\frac{1}{2}\left((q_{2}-x_{1})^{2}-(l_{10}+w)^{2}\right)=0,$
(38)

such that $m=1$. The total kinetic energy is given by

$\displaystyle T=\frac{1}{2}m_{1}v_{1}^{2}+\frac{1}{2}m_{2}(v_{2}+v_{3})^{2}$
(39)

and thus the mass matrix in (1) can be identified as

$\displaystyle\mbox{$\mathbf{M}$}=\begin{bmatrix}m_{1}&&0&&0\\\
0&&m_{2}&&m_{2}\\\ 0&&m_{2}&&m_{2}\end{bmatrix}$ (40)

and is singular. Thus, the inverse $\mbox{$\mathbf{M}$}^{-1}$ does not exist
and we cannot find a Hamiltonian for this problem. The potential is given by
the elastic potential of the two springs, which is assumed to be nonlinear,
such that

$\displaystyle
V(\mbox{$\mathbf{q}$})=\frac{1}{4}k_{1}(x_{1}^{2}+x_{1}^{4})+\frac{1}{4}k_{2}(x_{2}^{2}+x_{2}^{4}).$
(41)

Since $\mathbf{M}$ is constant and does not dependent of the configuration,
the discrete generalized energy (31) and the discrete total energy of the
system (35) are equivalent. Simulation parameters have been chosen as shown in
Table 1.

Table 1: Simulation parameters for mass-spring system with $i\in\\{1,2\\}$. $h$ | $T$ | $m_{i}$ | $k_{i}$ | $l_{i0}$ | $w$ | $\epsilon_{\mathrm{Newton}}$
---|---|---|---|---|---|---
$0.1$ | $10$ | $\\{1,1\\}$ | $\\{1,3\\}$ | $\\{1,1\\}$ | $0.1$ | $10^{-9}$

The simulation yielded the energy evolutions displayed in Fig. 2, where energy
transfer between kinetic and potential energy becomes visible. The energy-
consistent approximation provided by the proposed method becomes obvious in
Fig. 3, where the temporal increment of the generalized energy is close to
computer precision, such that it is a discretely conserved quantity.
Furthermore, Fig. 4 underlines the computationally exact treatment of the
constraint equation in each time step.

$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$0$$0.5$$1$$t$$T$$V$$E$ Figure 2: Energy
quantities $0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$-2$$-1$$0$$1$$2$$\cdot
10^{-15}$$t$$E^{n+1}-E^{n}$ Figure 3: Energy increments
$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$-2$$-1$$0$$1$$2$$\cdot
10^{-15}$$t$$g(\mbox{$\mathbf{q}$}^{n})$ Figure 4: Positional constraint (38)

###  6.2 Nonlinear spring pendulum with spherical coordinates

$\varphi$
---
$\theta$
---
$m$
---
$r$
---
$EA$
---
Figure 5: Spring pendulum.

We analyze the 3D spring pendulum with nonlinear spring potential formulated
in spherical coordinates. The distance of mass $m$ from the Cartesian origin
is given by $r\in\mathbb{R}_{\geq 0}$ and two angles $\theta,\varphi$ describe
the orientation, as depicted in Fig. 5. Consequently, the system has three
degrees of freedom and

$\displaystyle\mbox{$\mathbf{q}$}=\begin{bmatrix}r\\\ \theta\\\
\varphi\end{bmatrix}.$ (42)

The potential energy is given by the nonlinear internal potential with spring
constant $EA$ such that

$\displaystyle V(\mbox{$\mathbf{q}$})=\frac{1}{2}EA\epsilon^{2}(r),$ (43)

where the strain in the spring is computed via

$\displaystyle\epsilon(r)=\frac{r^{2}-l_{0}^{2}}{2l_{0}^{2}}$ (44)

and $l_{0}$ denotes the spring’s resting length. The mass matrix from (1) for
this example is given by

$\displaystyle\mbox{$\mathbf{M}$}(\mbox{$\mathbf{q}$})=m\,\mathrm{diag}(1,r^{2},r^{2}\sin^{2}(\theta))$
(45)

and becomes singular if $\theta=k\pi$ for all integers $k\in\mathbb{Z}$.
Inital conditions have been chosen as

$\displaystyle\mbox{$\mathbf{q}$}^{0}=\begin{bmatrix}1.05\\\ \pi/2\\\
0\end{bmatrix},\quad\mbox{$\mathbf{v}$}^{0}=\begin{bmatrix}0\\\ 1\\\
1\end{bmatrix}$ (46)

to achieve an initial elongation of the spring by $5\%$ and a tangential
initial velocity. Simulation parameters have been chosen as shown in Table 2.

Table 2: Simulation parameters for spring pendulum. $h$ | $T$ | $m$ | $EA$ | $l_{0}$ | $\epsilon_{\mathrm{Newton}}$
---|---|---|---|---|---
$0.01$ | $1$ | $1$ | $300$ | $1$ | $10^{-9}$

In this example, the energetic quantities have been simulated as can be seen
in Fig. 6. While the conservation of the discrete total energy is captured
down to an order of $10^{-4}$ (cf. Fig. 7), the generalized energy function is
exactly preserved by the newly proposed scheme (see Fig. 8).

$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$0$$0.5$$1$$1.5$$t$$T$$V$$E_{\mathrm{tot}}$
Figure 6: Energy quantities
$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$-2$$-1$$0$$1$$2$$\cdot
10^{-4}$$t$$E_{\mathrm{tot}}^{n+1}-E_{\mathrm{tot}}^{n}$ Figure 7: Discrete
total energy (35) increments
$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$-2$$-1$$0$$1$$2$$\cdot
10^{-14}$$t$$E^{n+1}-E^{n}$ Figure 8: Energy function (31) increments

##  7 CONCLUSION & OUTLOOK

The present approach to the numerical integration of dynamical systems is
based on Livens principle (cf. [1]). This variational principle is
characterized by combining Lagrangian and Hamiltonian viewpoints on mechanics.
Thus, the need for the inversion of the mass matrix is avoided, which allows
the simulation of systems with singular mass matrix. Moreover, the presented
framework takes account of configuration-dependent mass matrices. A
generalized energy function is introduced, which is identical to the total
energy and can be linked to the Hamiltonian if the inverse mass matrix exists.

The framework has been extended to mechanical systems subject to holonomic
constraints, and can thus be linked to previous works [7, 8], where a new
variational principle has been introduced. The restrictions therein to
constant and non-singular mass matrices have therefore been overcome.

Based on this novel formulation, a structure-preserving discretization has
been applied and the emanating time-stepping scheme features both velocity and
momentum quantities. As the scheme makes use of Gonzalez discrete gradients
(see, e.g. [10]), the proposed integrator discretely covers the conservation
of the generalized energy (which differs from the total energy of the system
if the mass matrix is configuration-dependent) and aims at the preservation of
momentum maps corresponding to the system’s symmetries. The integrator
satisfies the constraint functions down to computer precision. Moreover, due
to the relation with the midpoint rule, the scheme is at most second-order
accurate.

The newly proposed structure-preserving integration scheme enhances the method
presented by Gonzalez [15] with respect to the formulation in a more general
(not necessarily Hamiltonian) framework, which makes possible the energy-
consistent simulation of systems with configuration-dependent or singular mass
matrices.

The numerical properties of the present method have been demonstrated in
representative examples of multibody dynamics: the well-known spring pendulum
in spherical coordinates and a mass-spring system with singular mass matrix
taken from [5].

It should be noted that in the present work we have focussed on the design of
an energy-preserving scheme. The development of variational integrators based
on Livens principle with configuration dependent or singular mass matrices is
interesting, as well. Especially the design of higher order methods by
following the ideas presented in Wenger et al. [16] might be in the scope of
future research . In this connection, the approach recently developed in
Altmann & Herzog [17] should also be applicable.

## ACKNOWLEDGMENTS

We are thankful for insightful discussions with Felix Zähringer and Moritz
Hille. Financial support for this work by the DFG (German Research Foundation)
– project numbers 227928419 and 442997215 - is greatfully acknowledged.

## References

  * [1] Livens, G.H.:  IX. – On Hamilton’s principle and the modified function in analytical dynamics.  Proceedings of the Royal Society of Edinburgh 39 (1920) 113–119
  * [2] Bou-Rabee, N., Marsden, J.E.:  Hamilton–Pontryagin integrators on Lie groups, Part I: Introduction and structure-preserving properties.  Foundations of Computational Mathematics 9(2) (2009) 197–219
  * [3] Holm, D.D.:  Geometric mechanics: Part II: Rotating, translating and rolling. 2 edn.  Imperial College Press (2011)
  * [4] Pontryagin, L., Boltayanskii, V., Gamkrelidze, R., Mishchenko, E.:  The mathematical theory of optimal processes.  Wiley (1962)
  * [5] Udwadia, F., Phohomsiri, P.:  Explicit equations of motion for constrained mechanical systems with singular mass matrices and applications to multi-body dynamics.  Proc. Roy. Soc. London A 462 (2006) 2097–2117
  * [6] García de Jalón, J., Gutiérrez, M.:  Multibody dynamics with redundant constraints and singular mass matrix: existence, uniqueness, and determination of solutions for accelerations and constraint forces.  Multibody System Dynamics 30(3) (2013) 311–341
  * [7] Kinon, P., Betsch, P., Schneider, S.:  The GGL variational principle for constrained mechanical systems.  Multibody Syst Dyn (2023) doi: 10.1007/s11044-023-09889-6.
  * [8] Kinon, P., Betsch, P., Schneider, S.:  Structure-preserving integrators based on a new variational principle for constrained mechanical systems.  Nonlinear Dynamics (under revision) (2023)
  * [9] Hairer, E., Lubich, C., Wanner, G.:  Geometric numerical integration: structure-preserving algorithms for ordinary differential equations.  Springer (2006)
  * [10] Gonzalez, O.:  Time integration and discrete Hamiltonian systems.  Journal of Nonlinear Science 6 (1996) 449–467
  * [11] Pars, L.A.:  A Treatise on Analytical Dynamics.  The Mathematical Gazette 50(372) (1966) 226–227
  * [12] Kunkel, P., Mehrmann, V.:  Differential-algebraic equations: analysis and numerical solution.  EMS textbooks in mathematics. European Mathematical Society (2006)
  * [13] Gear, C.W., Leimkuhler, B., Gupta, G.K.:  Automatic integration of Euler-Lagrange equations with constraints.  Journal of Computational and Applied Mathematics 12-13 (1985) 77–90
  * [14] Kinon, P.L., Bauer, J.K.:  Metis: Computing constrained dynamical systems, version 1.0.7 (2023) GitHub repository, doi: 10.5281/zenodo.7769954.
  * [15] Gonzalez, O.:  Mechanical systems subject to holonomic constraints: Differential–algebraic formulations and conservative integration.  Physica D: Nonlinear Phenomena 132(1-2) (1999) 165–174
  * [16] Wenger, T., Ober-Blöbaum, S., Leyendecker, S.:  Construction and analysis of higher order variational integrators for dynamical systems with holonomic constraints.  Advances in Computational Mathematics 43(5) (2017) 1163–1195
  * [17] Altmann, R., Herzog, R.:  Continuous Galerkin schemes for semiexplicit differential-algebraic equations.  IMA Journal of Numerical Analysis 42(3) (2022) 2214–2237


Proceedings of the CTD 2023 PROC-CTD2023-27

Combined track finding with GNN & CKF

Lukas Heinrich1, Benjamin Huth2, Andreas Salzburger3, Tilo Wettig2

1Technical University of Munich, 2University of Regensburg, 3CERN

> ABSTRACT
>
>
>
>
> The application of Graph Neural Networks (GNN) in track reconstruction is a
> promising approach to cope with the challenges arising at the High-
> Luminosity upgrade of the Large Hadron Collider (HL-LHC). GNNs show good
> track-finding performance in high-multiplicity scenarios and are naturally
> parallelizable on heterogeneous compute architectures.
>
> Typical high-energy-physics detectors have high resolution in the innermost
> layers to support vertex reconstruction but lower resolution in the outer
> parts. GNNs mainly rely on 3D space-point information, which can cause
> reduced track-finding performance in the outer regions.
>
> In this contribution, we present a novel combination of GNN-based track
> finding with the classical Combinatorial Kalman Filter (CKF) algorithm to
> circumvent this issue: The GNN resolves the track candidates in the inner
> pixel region, where 3D space points can represent measurements very well.
> These candidates are then picked up by the CKF in the outer regions, where
> the CKF performs well even for 1D measurements.
>
> Using the ACTS infrastructure, we present a proof of concept based on truth
> tracking in the pixels as well as a dedicated GNN pipeline trained on
> $t\overline{t}$ events with pile-up 200 in the OpenDataDetector.

> PRESENTED AT
>
>
>
>
> Connecting the Dots Workshop (CTD 2023)
>
> October 10-13, 2023

##  1 Introduction and motivation

The upcoming upgrade of the _Large Hadron Collider_ (LHC) towards higher
luminosities (HL-LHC) is a major challenge for track-reconstruction software.
It is expected that the number of simultaneous proton-proton collisions (pile-
up) will increase from currently $\left\langle\mu\right\rangle=50$ in LHC Run
3 up to $\left\langle\mu\right\rangle=200$, which will bring the current CPU-
based reconstruction software to its limit and beyond. For that reason, much
effort is put into exploring new algorithms, especially machine-learning
algorithms based on neural networks because these are parallelizable very
efficiently on heterogeneous architectures like GPUs.

Recently, _Graph Neural Network_ (GNN)-based algorithms have shown very
promising results for track finding [8, 10]. However, in their basic version,
their performance suffers from low spatial resolution in detector regions like
strip modules.

Proposed solutions to this problem include more complex network architectures
like heterogeneous GNNs [5]. However, this work explores a different approach:
The idea is to use the GNN for track finding in the inner, high-resolution
parts of the detector, where it was shown to be able to resolve the track
candidates very efficiently [6]. For the outer parts, we rely on a standard
algorithm for track finding, the _Combinatorial Kalman Filter_ (CKF),111The
CKF is a combined track-finding and track-fitting algorithm that
combinatorially explores track candidates by considering multiple measurements
per detector module. for which the GNN track candidates serve as seeds. For
the track fitting, we use the CKF both in the inner and the outer parts. From
combining GNN and CKF, we expect several benefits:

  * •

Improved seed quality. Triplet-seeding algorithms usually produce many
duplicate seeds that result in duplicate tracks. GNN seeds can be larger than
3 hits and are unlikely to produce duplicate seeds.222In the case of no
merged-hit clusters, duplicate seeds with shared hits are not allowed by
construction. This should improve the computational performance and avoid
track cleaning by an ambiguity resolution. afterward.

  * •

Less branching. By construction, we do not allow the CKF to branch in the
pixel region, where it follows the track candidate predicted by the GNN. We
also expect that due to high-quality seeds, we can reduce the branching in the
outer regions. This will have similar effects on duplication rate and
performance as the previous benefit.

  * •

Natural handling of 1D measurements. The CKF has a natural, mathematical
consistent way of handling 1D measurements, whereas GNN architectures, at
least currently, rely on 3D space points that must be constructed from
measurements and the detector geometry. For a precise 3D-spacepoint
construction, a 2D measurement is necessary.

  * •

Smaller graph construction. Graph construction is a critical part of a GNN
pipeline and often expensive in both runtime and memory. By restricting the
graph to the pixels, these costs are reduced naturally.

Tests of this idea showed room for improvement of the track-finding
performance in some regions. Such improvement was achieved by also including
hits from the short-strip-barrel section in the GNN.

##  2 Experimental setup

For the present study, a track reconstruction chain that includes a GNN track-
finding module for the pixels and a subsequent CKF that is seeded with the GNN
track candidates has been implemented using the ACTS toolkit [11]. In the
following, we refer to this chain as the GNN+CKF chain. For the detector
geometry, the OpenDataDetector (ODD) is used [2], a virtual, hermetic detector
with 3 parts: inner, high-resolution pixel modules, a short-strip section, and
double-sided long-strip modules in the outermost part.

In addition to this novel tracking chain, several other chains were studied
for comparison:

  * •

The _proof-of-concept_ chain is very similar to the GNN+CKF chain, but the GNN
is replaced by a truth-track-finder. This is helpful for figuring out whether
inefficiencies come from the GNN or the CKF part of the main chain.

  * •

A _standard CKF_ chain with classical triplet seeding allows for comparison to
state-of-the-art tracking tools. However, we did not tune this chain in this
study, and thus we do not expect cutting-edge performance.

  * •

As an upper performance bound, a _Kalman truth-tracking_ chain has been
implemented as well.

For this study, $t\overline{t}$-events with pile-up 200 were generated using
the Pythia8 [4] event generator, followed by a Geant4 [1] full simulation in
the ODD. To emulate the detector readout, a hybrid approach was chosen:

  * •

For the region covered by the GNN, a _geometric digitization_ was applied.
This digitization implements a basic simulation of the individual cell
activation based on the path of the particle through the module, and thus can
provide cluster information in addition to the hit coordinates. The only major
simplification in this setup is that cluster merging was disabled. This means
that each cluster contains only one hit and no shared hits are possible. This
simplifies the mapping of a collision event to a graph.

  * •

For the remaining regions, a simple digitization based on Gaussian smearing of
the truth hits was used, both because the CKF does not directly use cluster
information and to limit the effort for validating the digitization
configuration.

We implemented a GNN-based track-finding pipeline as follows: The first
element is a _graph-construction_ stage based on metric learning. This is
followed by three _edge-classification_ stages (one simple filter based on a
_multi-layer perceptron_ (MLP) and two GNN stages333It turned out that the
second GNN stage improved the performance. We suspect that the improvement is
due to a more balanced graph in the second GNN stage. It would be interesting
to explore this issue in future work.). The final element is a _track-
building_ stage, in which track candidates are constructed from the GNN output
by a simple connected-components algorithm.

The two GNN stages use a modified version of the _Interaction GNN_ that is
provided in the GNN4ITk Common Framework [3]. The major differences are the
implementation of the undirected graph and the recurrent implementation of the
message-passing step.

##  3 Training

The training is done with a custom branch of the GNN4ITk Common Framework,
which provides the infrastructure to train and evaluate multi-stage pipelines
based on PyTorch [9].

For the training of the GNN module, 2000 events with the configuration
described above were simulated, of which 250 events each were dedicated for
testing and validation purposes, respectively. The 3D space points in
cylindrical coordinates $(r,\varphi,z)$ are used as input features to the
neural networks. We also used additional features based on the cluster shape
(_cell count_ , _sum of activations_ and the _cluster size_ in the surface-
local coordinates) in the graph-construction stage because they lead to
performance improvements in this stage. All input features are scaled such
that they are in the range between 0 and 1. Particle trajectories that loop
inside the detector volume were split at the turning points before generating
the ground truth so that the network does not have to learn signatures of
looping particles.

In order to focus the training on particles of interest, edges of particles
with $p_{T}>$1\text{\,}\mathrm{GeV}$$ and at least 3 hits in the pixels were
assigned higher weight. We call this _target selection_ in the following.

The classifier stages each return a score between 0 and 1 for every edge. In
the different stages, different score cuts were used to classify the edges
into positive and negative edges. For the MLP filter and the first GNN, we
have chosen a very low score cut to preserve a very high efficiency and only
remove edges that are classified as negative with high confidence. In the last
stage, a balanced score cut of $0.5$ was used.

To quantify the performance of different pipeline stages, we define metrics
based on classified edges,

$\text{ edge efficiency}=\frac{\text{true positive edges}}{\text{true
edges}}\,,\qquad\text{edge purity}=\frac{\text{true positive
edges}}{\text{positive edges}}\,.$ (1)

The resulting graph is perfect if both metrics are equal to $1$. We can
evaluate these metrics for both the whole set of edges and the target edges.
Because we do not expect the pipeline to be a perfect classifier between
target and non-target particles, both _target efficiency_ and _overall purity_
should be close to 1.

Figure 1: Target-edge efficiency of the different stages as a function of
$\eta$, $p_{T}$ and the number of pixel hits of the corresponding particle for
the pipeline trained with the pixel hits only

We reach an overall purity of $>0.99$ at the second GNN, which is close to
optimal. Therefore, we focus on the target-edge efficiency in the following
(see Fig. 1).

We observe a very high and mostly flat performance for the first two
classifier stages. However, at the second GNN, the performance drops
significantly for $|\eta|\approx 1$ and $p_{T}>$10\text{\,}\mathrm{GeV}$$. A
possible explanation is the number of hits available to the GNN: The ODD has
only 4 pixel layers. For 3 hits in roughly the same $\eta$-slice, each
combination of hits is physically valid. Correspondingly, we only observe an
efficiency of $\approx 0.75$ in this case. Even a fourth hit might not provide
enough discrimination power at high $p_{T}$.

Figure 2: Same as Fig. 1 but trained with the pixel hits and the hits in the
short-strip barrel.

To provide further evidence for this explanation, we trained the pipeline
again, this time also including the short-strip-barrel section (providing 4
extra layers of hits, albeit with greatly reduced resolution). Indeed, we see
much better performance (see Fig. 2). Both the efficiency in the central-
barrel region and in the high-$p_{T}$ region have greatly improved. We will
mainly focus on this GNN pipeline in the following.

##  4 Inference and performance analysis

The implementation of the pipeline in ACTS is based on the C++ interface to
PyTorch. A modified CKF algorithm was implemented, which accepts the GNN track
candidates with their corresponding start parameters as input. The start
parameters are computed from a 3-point estimate from the first, the last and
the middle hit of the track candidate. This allows the CKF to only use hits
from these track candidates in the pixel region and the short-strip-barrel
region. All code to reproduce our results can be found in [7].

To measure the track-finding performance, we use metrics based on the
following selection and matching criteria. We select particles if they have
$p_{T}>$1\text{\,}\mathrm{GeV}$$, at least 7 hits in the full detector, and at
least 3 hits in the pixels. We select tracks if they have a fitted
$p_{T}>$1\text{\,}\mathrm{GeV}$$ and at least 7 hits in the full detector. The
particles and tracks are matched using _double matching_ : We match a track to
a particle if the track contains more than $50\text{\,}\mathrm{\char
37\relax}$ of the hits of the particle (particle efficiency) and if more than
$50\text{\,}\mathrm{\char 37\relax}$ of the hits in the track belong to the
particle (track purity). Otherwise, it is a _fake track_.444The fake track is
assigned to its majority particle (based on track purity) for further
analysis, see Figs. 3 and 5. Our performance metrics are

matching efficiency
$\displaystyle=\frac{N_{\text{particles}}(\text{selected},\text{matched})}{N_{\text{particles}}(\text{selected})}\,,$
(2) duplication rate
$\displaystyle=\frac{N_{\text{particles}}(\text{selected},\text{matched
tracks}>1)}{N_{\text{particles}}(\text{selected})}\,,$ (3) fake rate
$\displaystyle=\frac{N_{\text{particles}}(\text{selected},\text{unmatched
tracks}>0)}{N_{\text{particles}}(\text{selected})}\,.$ (4)

These three metrics measure different aspects of the track-finding
performance. In general, we seek high matching efficiency and low fake and
duplication rate.

Figure 3: Efficiency, duplication rate, and fake rate vs particle $\eta$ and
$p_{T}$, evaluated with 50 events.

Fig. 3 shows the track-finding performance for 50 $t\overline{t}$-events
(generated independently of the training set) as histograms of the particle
$\eta$. The matching efficiency is shown in Fig. 3 (left). In the central
barrel region, we observe a mostly flat performance, with the performance of
the standard CKF and the proof-of-concept chain slightly above the GNN+CKF
chain. This suggests small inefficiencies in the GNN+CKF chain. For the
standard CKF chain, we observe efficiency drops at $|\eta|\approx 2$ and
$|\eta|>3$. This is also reflected in the fake rate and suggests an
inefficiency in the classical triplet seeding.

(a) GNN+CKF chain, pixel-only pipeline.

(b) GNN+CKF chain, short-strip-extended pipeline.

Figure 4: Detailed matching efficiency for two versions of the GNN+CKF chain,
evaluated with 50 events.

We can resolve the inefficiencies in the GNN+CKF chain in more detail (see
Fig. 4). We first disregard particles that are not matched in the proof-of-
concept chain555Here the failure to match occurs in the CKF part of the chain
because the truth-track finder is perfect by construction. (shown in gray).
Then, we identify particles for which no matching track candidate was found by
the GNN (shown in red). This red area is rather large for the pixel-only GNN
pipeline (Fig. 4(a)) but almost disappears for the pipeline that includes the
hits in the short-strip barrel (Fig. 4(b)). Finally, some particles with a
valid GNN track candidate are not matched after the CKF stage (shown in
yellow), mainly because the track candidates are not pure enough to produce
high-quality seeds.666For example, there could be an additional hit included
in the track seed that leads to a wrong estimation of the start parameters of
the CKF fit. This could be improved by applying more advanced track-building
algorithms.

We return to Fig. 3 and discuss the duplication rate (middle). As expected,
the standard CKF has a high duplication rate, which is close to 1 over the
whole $\eta$ range. In the proof-of-concept and GNN+CKF chains, duplication is
caused mainly by the CKF part. The observation that both chains have similar
duplication rates implies that the GNN track candidates, compared to the
perfect proof-of-concept candidates, do not lead to a higher duplication rate.

To check if we can reduce this duplication rate further without efficiency
costs, we configured the CKF to only accept a single measurement per surface,
effectively reducing the CKF to a simple Kalman Filter. The chosen measurement
is the one with the lowest $\chi^{2}$ contribution. The resulting performance
is shown in Fig. 5. We observe that we can remove the combinatorial aspect of
the fitter with negligible matching-efficiency cost. Furthermore, the fake
rate is also lower in this setup.

Figure 5: Metrics for the GNN+CKF chain with and without the combinatorial
aspect in the CKF stage.

##  5 Computational performance

(a) Timing comparison of different chains.

(b) Zoom into the timing of the GNN+CKF chain.

Figure 6: Average per-event timing for 10 events. MCKF is short for “Modified
CKF”.

To measure the timing of the different chains, we executed the workflows in a
single-threaded environment with a single available GPU. The system comprised
an AMD EPYC 7662 64-core processor and an Nvidia A100 with 40GB memory. The
results are shown in Fig. 6.

The GNN+CKF chain can run the reconstruction of a $t\overline{t}$ event with
pile-up 200 in the ODD in $\approx$2\text{\,}\mathrm{s}$$. We observe a
slightly increased runtime when additional hits from the short-strip barrel
are included in the GNN. The standard CKF shows an unreasonably large wall
time. This drastic difference in timing is caused by the high duplication rate
and reminds us that it is important to keep the combinatorics under control
during reconstruction.777Because the GNN+CKF chain uses the same
infrastructure, we can exclude code issues as the problem here.

##  6 Conclusion

We have presented a complete track-reconstruction chain that uses a novel
hybrid approach by combining a GNN for track finding in the inner detector
regions and a subsequent CKF for track finding in the outer parts and fitting.

We observed that using the GNN in the pixels limits its performance in the
central detector region. This issue was mitigated by including additional hits
from the short-strip-barrel section. This finding suggests that there is a
minimum number of hits necessary for the GNN to reach its full discriminative
power. It would be interesting to investigate this behavior further in
different setups.

We also showed that for good-quality GNN track candidates, the combinatorial
aspect in the subsequent CKF can be omitted without significant efficiency
loss. As long as Kalman-based fitting is needed, this offers an opportunity to
optimize runtime in scenarios with high throughput demands such as triggers.

ACKNOWLEDGEMENTS

LH is supported by the Excellence Cluster ORIGINS, funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy - EXC-2094-390783311. GPU resources were made available by
the ORIGINS Data Science Lab (ODSL). BH and TW are supported in part by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), project
number 460248186.

## References

  * [1] S. Agostinelli et al.  “Geant4—a simulation toolkit”  In _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 506.3, 2003, pp. 250–303  DOI: https://doi.org/10.1016/S0168-9002(03)01368-8
  * [2] Corentin Allaire et al.  “OpenDataDetector”  In _Zenodo_ , 2022  DOI: 10.5281/zenodo.6445359
  * [3] Markus Julian Atkinson et al.  “gnn4itk”  URL: https://github.com/GNN4ITkTeam/CommonFramework
  * [4] Christian Bierlich et al.  “A comprehensive guide to the physics and usage of PYTHIA 8.3”  arXiv, 2022  DOI: 10.48550/ARXIV.2203.11601
  * [5] Sylvain Caillou  “Novel fully-heterogeneous GNN designs for track reconstruction at the HL-LHC”, 2023
  * [6] Paolo Calafiura et al.  “Applying and optimizing the Exa.TrkX Pipeline on the OpenDataDetector with ACTS”  In _PoS_ ICHEP2022, 2022, pp. 227  DOI: 10.22323/1.414.0227
  * [7] Benjamin Huth  “https://github.com/benjaminhuth/acts-ckf-gnn-combination”, 2023
  * [8] Xiangyang Ju et al.  “Performance of a geometric deep learning pipeline for HL-LHC particle tracking”  In _The European Physical Journal C_ 81.10, 2021, pp. 876  DOI: 10.1140/epjc/s10052-021-09675-8
  * [9] Adam Paszke et al.  “PyTorch: An Imperative Style, High-Performance Deep Learning Library”  arXiv, 2019  DOI: 10.48550/ARXIV.1912.01703
  * [10] C. Rougier et al.  “CTD2022: ATLAS ITk Track Reconstruction with a GNN-based Pipeline”  Zenodo, 2023  DOI: 10.5281/zenodo.8119762
  * [11] Andreas Salzburger et al.  “acts-project/acts: v22.0.0”  Zenodo, 2022  DOI: 10.5281/zenodo.7467880


SNSN-323-63

Measurements of top-quark pair and single top total cross sections, and in-
situ systematic constraints

Olga Bessidskaia Bylund
on behalf of the ATLAS and CMS Collaborations111Copyright [2018] CERN for the
benefit of the ATLAS and CMS collaborations, CC-BY-4.0 license.

Bergische Universität Wuppertal, GERMANY

> Measuring the total cross sections of top-quark pair production and single
> top processes at high precision tests the predictions of the Standard Model
> and can bring better understanding of properties such as the top-quark mass,
> electroweak couplings, lepton universality and proton parton distribution
> functions. Some of the recent measurements by the ATLAS and CMS experiments
> at the Large Hadron Collider are outlined in this document. Experimental
> methods are given particular attention, particularly the determination of
> lepton isolation efficiencies in the ATLAS dilepton measurement of top-quark
> pair production.

> PRESENTED AT
>
>
>
>
> $13^{\mathrm{th}}$ International Workshop on Top Quark Physics
> Durham, UK (videoconference), 14–18 September, 2020

## Introduction

Measuring the total production cross sections of top-quark pairs
($t\overline{t}$) and in single top-quark production tests the Standard Model
and can help us to better understand various properties of the top-quark,
electroweak couplings and the structure of the proton. Some recent results
from ATLAS [1] and CMS[2] at the LHC are outlined below, with a focus on in-
situ techniques in the $t\overline{t}$ dilepton analysis by ATLAS.

## Top-quark pair production in the dilepton channel (ATLAS)

Measuring $t\overline{t}$ in the dilepton channel while using innovative
experimental techniques, the ATLAS experiment has pushed the precision on the
measured cross section down to 2.4%. A partial Run 2 dataset of 36.1
$\mathrm{fb}^{-1}$collected at 13 TeV collision energy in 2015-2016 is used.
In the selection, an electron-muon pair of opposite-sign charge is required
along with jets, of which either one or two are $b$-tagged. The equations for
the observed number of events in the 1 $b$-tag and 2 $b$-tag regions ($N_{1}$,
$N_{2}$) are constructed:

$\displaystyle N_{1}$ $\displaystyle=$ $\displaystyle
L\sigma_{t\overline{t}}\epsilon_{e\mu}2\epsilon_{b}(1-C_{b}\epsilon_{b})+N_{1}^{bkg}$
(1) $\displaystyle N_{2}$ $\displaystyle=$ $\displaystyle
L\sigma_{t\overline{t}}\epsilon_{e\mu}2C_{b}\epsilon_{b}^{2}+N_{2}^{bkg}$ (2)

and solved for $\sigma_{t\overline{t}}$ and the $b$-tagging efficiency
$\epsilon_{b}$. The correlation coefficient $C_{b}$ for reconstructing and
tagging the two $b$-jets is determined from simulation and $N_{1,2}^{bkg}$
denotes the background contribution. The efficiency to pass the dilepton
selection $\epsilon_{e\mu}$ is determined from simulation and then corrected
using data; it is lower in 2016 than in 2015. The contamination from QCD
multijet events where leptons fail the isolation requirements is estimated
using a selection where the cut on the isolation variable
$|d_{0}|/\sigma_{d_{0}}$ has been inverted. The isolation efficiencies are
corrected, accounting for this contamination, by up to 1% for electrons and
0.4% for muons, see Fig. 1.

Figure 1: The corrections for the lepton isolation efficiencies, binned in
lepton transverse momentum $p_{\mathrm{T}}$ and pseudorapidity $|\eta|$, for
the measurement of $t\overline{t}$ in the dilepton channel by the ATLAS
experiment [3]. The baseline from simulation is shown in turquoise and the
corrections from data in black.

Using the BLUE method [4, 5], the datasets from 2015 and 2016 are combined
with weights of 0.49 and 0.51 respectively (note that the 2016 dataset is much
larger), which improves the total sensitivity by 9% compared to treating them
as one dataset. Each source of systematic uncertainty is evaluated by changing
all relevant input quantities coherently and re-solving Eqs. (1), (2). The
cross section is measured to be:

$\displaystyle\sigma_{t\overline{t}}=826.4\pm 3.6\mathrm{(stat)}\pm
11.5\mathrm{(syst)}\pm 15.7\mathrm{(lumi)}\pm
1.9\mathrm{(beam)}\,\textrm{pb}\,,$

which is the most precise inclusive measurement of $t\overline{t}$ to date and
in agreement with

$\displaystyle\sigma_{t\overline{t}}^{theo}=832\pm
35\mathrm{(pdf+\alpha_{S})}^{+20}_{-29}\mathrm{(scale)}\,\textrm{pb}\,.$

The sensitivity of this measurement is limited by the precision on the
luminosity, amounting to an uncertainty of 1.9% on the cross section. The
total systematic uncertainty of 1.39% has the main contributions from the
normalisation of the $tW$ background and the modelling uncertainties for
$t\overline{t}$. The dilepton analysis is a priori sensitive to modelling
uncertainties, but their impact has been reduced by measuring the lepton
isolation efficiencies in-situ: when comparing the predictions from different
generators, the different predictions for these efficiencies need not be
included.

Additionally, the ratio of $t\overline{t}$ to $Z$ production is evaluated to
nearly cancel the large dependency on the luminosity and constrain proton
p.d.fs. Ratios of $t\overline{t}$ production at different collision energies
and double ratios of $t\overline{t}$ to $Z$ production at different energies
are also computed, leading to further cancellations.

## Top-quark pair production in the dilepton channel (CMS)

In the $t\overline{t}$ cross section measurement in the dilepton channel by
CMS [6], in addition to $e^{\pm}\mu^{\mp}$, the opposite-sign $ee$ and
$\mu\mu$ channels are considered. 35.9 $\mathrm{fb}^{-1}$ of data collected at
13 TeV collision energy is used for the measurement. A total of 28 event
categories are defined by splitting by jet and $b$-jet multiplicity in these
channels: with 0, 1 or 2 $b$-tagged jets and 0, 1, 2 or at least 3 additional
un-tagged jets (the no $b$-jet categories are only used for $e\mu$ events). A
simultaneous template fit using all event categories is performed, fitting the
transverse momentum $p_{\mathrm{T}}$ of the additional jet with highest
$p_{\mathrm{T}}$, if defined, otherwise the event yield is fitted. A priori,
the identification efficiency is known to better precision for muons than for
electrons. By using the different lepton channels, the electron identification
efficiency is constrained to the precision for muons. Using categories with
different $b$-jet multiplicities allows constraining the efficiency for
identifying and selecting a $b$-jet. The total cross section is measured to:

$\displaystyle\sigma_{t\overline{t}}=803\pm 2\mathrm{(stat)}\pm
25\mathrm{(syst)}\pm 20\mathrm{(lumi)}\,\textrm{pb}\,,$

in agreement with theory. The main uncertainties relate to the integrated
luminosity and lepton identification and isolation, with an impact of 2.5% and
2.0% respectively.

## Top-quark pair production, lepton+jets (ATLAS)

The ATLAS experiment has performed a measurement of the $t\overline{t}$ cross
section in the lepton+jets channel using the full Run 2 dataset (139
$\mathrm{fb}^{-1}$) at 13 TeV [7]. Exactly one charged lepton is required as
well as at least four jets. Three event categories are defined, split by jet
and $b$-jet multiplicity with: 4 jets with 1 $b$-tagged in the SR1 category, 4
jets of which 2 are $b$-tagged in SR2 and at least 5 jets of which 2 are
$b$-tagged in SR3. The aplanarity is used as the fitting variable in SR1, the
minimum invariant mass over all lepton-jet pairs in SR2 and a distance
parameter between the jets in SR3. The result is:

$\displaystyle\sigma_{t\overline{t}}=830\pm 0.4\mathrm{(stat)}\pm
36\mathrm{(syst)}\pm 14\mathrm{(lumi)}\,\textrm{pb},$

with a precision of 4.6% and in agreement with the theoretical prediction.
This is the most precise inclusive $t\overline{t}$ measurement in this
channel. The measurement is systematically limited, with the main
uncertainties originating from parton shower systematics, followed by
luminosity and final state radiation.

## Top-quark pair production with a hadronic $\tau$ (CMS)

CMS has measured $t\overline{t}$ production in a selection targeting a
hadronically decaying $\tau$ lepton and one electron or muon, using 35.9
$\mathrm{fb}^{-1}$ of data collected at 13 TeV collision energy [8]. This
channel corresponds to around 5% of all $t\overline{t}$ final states and
serves as a test of lepton universality. One electron or muon, a hadronic
$\tau$ candidate of opposite-sign charge and two jets are required, of which
at least one should be $b$-tagged. The selection is further split into
categories by cutting on a distance parameter between the jets.

The $\tau$ leptons are reconstructed with the hadron-plus-strips (HPS)
algorithm. In each jet, a charged hadron is combined with other nearby charged
hadrons or photons and the decay modes of the $\tau$ are identified. Electrons
and photons are clustered in strips along the bending direction of the
trajectory, which enhances the identification of $\pi^{0}$ mesons by taking
into account the early showering of photons. A boosted decision tree (BDT) is
trained to separate HPS $\tau$ candidates from QCD multijet events. The
largest background in the measurement comes from $t\overline{t}$ events with a
jet misidentified as originating from a $\tau$. The uncertainty in the
efficiency of $\tau$ identification is 5% and has the highest impact on the
precision. To extract the cross section, the transverse mass
$M_{T}(\ell,E_{\mathrm{T}^{\mathrm{miss}}})$ is fitted in a profile likelihood
fit, giving:

$\displaystyle\sigma_{t\overline{t}}=781\pm 7\mathrm{(stat)}\pm
62\mathrm{(syst)}\pm 20\mathrm{(lumi)}\,\textrm{pb}\,,$

which is consistent with the Standard Model prediction. This is the first
measurement of $t\overline{t}$ production with a $\tau$ lepton in the final
state performed at 13 TeV. The ratio to $t\overline{t}$ production in the
dilepton channel is computed and found to be in agreement with lepton
universality.

## Single top production in the $t$-channel (CMS)

The CMS experiment has measured the single top production cross section in the
$t$-channel [9], using 35.9 $\mathrm{fb}^{-1}$ of data collected at 13 TeV.
This process is charge asymmetric due to the valence quark content of the
colliding protons. An electron or muon is required in the selection together
with jets. In the signal-enriched category exactly one $b$-tagged jet and one
additional jet are required. The cross sections are measured separately for
top and top antiquark production by splitting the events into categories by
the charge of the lepton and extracting them in a simultaneous fit. Their
ratio $R_{t}=\sigma_{t}/\sigma_{\overline{t}}$, which is sensitive to proton
p.d.fs, is also determined, see Fig. 2.

Figure 2: Feynman diagrams illustrating $t$-channel production and the
measured value for $R_{t}$ by the CMS experiment [9], compared to different
p.d.fs.

BDT algorithms are employed to enhance the separation of the signal from the
large backgrounds consisting mostly of $t\overline{t}$ and $W$+jets events.
Experimental uncertainties and uncertainties on the background rates and
modelling enter as nuisance parameter in the profile likelihood fit, which
allows them to be constrained. The uncertainties relating to signal modelling
are not profiled to avoid propagating any constraints on them to the full
phase space, where they may not be valid. The cross sections and their ratio
are measured to be:

$\displaystyle\sigma_{t}=$ $\displaystyle 130\pm 1\mathrm{(stat)}\pm
19\mathrm{(syst)}\,\textrm{pb}$ , $\displaystyle\sigma_{\overline{t}}=$
$\displaystyle 77\pm 1\mathrm{(stat)}\pm 12\mathrm{(syst)}\,\textrm{pb}$
$\displaystyle,\quad R_{t}=1.68\pm 0.02\mathrm{(stat)}\pm
0.05\mathrm{(syst)}\,.$

This is consistent with the predictions from the Standard Model:

$\displaystyle\sigma_{t}^{theo}$
$\displaystyle=136.0^{+4.1}_{-2.9}\mathrm{(scale)}\pm
3.5\mathrm{(pdf+\alpha_{S})}\,\textrm{pb}$ ,
$\displaystyle\sigma_{\overline{t}}^{theo}$
$\displaystyle=81.0^{+2.5}_{-1.7}\mathrm{(scale)}\pm
3.2\mathrm{(pdf+\alpha_{S})}\,\textrm{pb}$ $\displaystyle,\quad R_{t}=1.68\,.$

## Single top combination and measurements of $\mathbf{tW}$

ATLAS and CMS have performed a combination of single top measurements at 7 and
8 TeV [10] using the BLUE method. The results are summarised in Fig. 3.

Figure 3: The combination of single top measurements by ATLAS and CMS at 7 and
8 TeV [10]. The theoretical prediction is shown in gray, the ATLAS result in
blue, the CMS result in red and the combination in black.

Since then, ATLAS also published a measurement of $tW$ production at 8 TeV in
the challenging lepton+jets channel [11], using 20.2 $\mathrm{fb}^{-1}$of
data, obtaining:

$\displaystyle\sigma_{tW}=26\pm 7\,\textrm{pb}\,,$

with a statistical uncertainty of 4 pb. The result agrees with the theoretical
prediction:

$\displaystyle\sigma_{tW}^{theo}=22.4\pm 0.6\mathrm{(scale)}\pm
1.4\mathrm{(pdf)}\,\textrm{pb}\,.$

CMS reported the most precise measurement of $tW$ at 13 TeV [12] with 35.9
$\mathrm{fb}^{-1}$:

$\displaystyle\sigma_{tW}=63.1\pm 1.8\mathrm{(stat)}\pm 6.4\mathrm{(syst)}\pm
2.1\mathrm{(lumi)}\,\textrm{pb}\,,$

consistent with

$\displaystyle\sigma_{tW}^{theo}=71.7\pm 1.8\mathrm{(scale)}\pm
3.4\mathrm{(pdf)}\,.$

## Summary

A brief overview has been given of some recent total cross section
measurements of $t\overline{t}$ and single top production by ATLAS and CMS at
the LHC, with unprecedented precision reached in several channels. Particular
attention is paid to the in-situ techniques in the ATLAS dilepton measurement
that allows for the reduction of systematic uncertainties. A highlight from
CMS is the first measurement of $t\overline{t}$ production at 13 TeV with
selection targeting a $\tau$ lepton. Comparisons and combinations of the
results from ATLAS and CMS help verify the results and improve the precision.

## References

  * [1] ATLAS Collaboration, 2008 JINST 3 S08003.
  * [2] CMS Collaboration, 2008 JINST 3 S08004.
  * [3] ATLAS Collaboration, Eur.Phys.J.C 80 (2020) 6, 528.
  * [4] L. Lyons, D. Gibaut, P. Clifford, Nucl. Instrum. Methods A 270, 110 (1988)
  * [5] A. Valassi, Nucl. Instrum. Methods A 500, 391 (2003).
  * [6] CMS Collaboration, Eur.Phys.J.C 79 (2019) 5, 368.
  * [7] ATLAS Collaboration, Phys. Lett. B 810 (2020) 135797.
  * [8] CMS Collaboration, JHEP 02 (2020) 191.
  * [9] CMS Collaboration, Phys.Lett.B 800 (2020) 135042.
  * [10] ATLAS Collaboration and CMS Collaboration, JHEP 05 (2019) 088.
  * [11] ATLAS Collaboration, 2007.01554, submitted to EPJC.
  * [12] CMS Collaboration, JHEP 10 (2018) 117.