# M-LIO: Multi-lidar, multi-IMU odometry with sensor dropout tolerance

Sandipan Das1,2, Navid Mahabadi3, Maurice Fallon4, Saikat Chatterjee1 1 KTH
EECS, Sweden<EMAIL_ADDRESS>
2 Scania, Sweden<EMAIL_ADDRESS>
3 Stockholm, Sweden<EMAIL_ADDRESS>
4 Oxford Robotics Institute, UK<EMAIL_ADDRESS>
We present a robust system for state estimation that fuses measurements from
multiple lidars and inertial sensors with GNSS data. To initiate the method,
we use the prior GNSS pose information. We then perform incremental motion in
real-time, which produces robust motion estimates in a global frame by fusing
lidar and IMU signals with GNSS translation components using a factor graph
framework. We also propose methods to account for signal loss with a novel
synchronization and fusion mechanism. To validate our approach extensive tests
were carried out on data collected using Scania test vehicles (5 sequences for
a total of $\approx$ 7 Km). From our evaluations, we show an average
improvement of 61% in relative translation and 42% rotational error compared
to a state-of-the-art estimator fusing a single lidar/inertial sensor pair.

##  I Introduction

State estimation, which is a sub-problem of Simultaneous Localization and
Mapping (SLAM), is a fundamental building block of autonomous navigation. To
develop robust SLAM systems, proprioceptive (IMU – inertial measurement unit,
wheel encoders) and exteroceptive (camera, lidar, GNSS – Global Navigation
Satellite System) sensing are fused. Existing systems have achieved robust and
accurate results [1, 2, 3]; however, state estimation in dynamic conditions
and when there is measurement loss or noisy is still challenging.

Compared to visual SLAM, lidar-based SLAM has higher accuracy as lidar range
measurements of up to 100m directly enable precise motion tracking. Moreover,
owing to the falling cost of lidars, mobile platforms are increasingly
equipped with multiple lidars [4, 5], to give complementary and complete
$360^{\circ}$ sensor coverage (see Fig. 1, for our data collection platform).
This also improves the density of measurements which may be helpful for state
estimation in degenerate scenarios such as tunnels or straight highways. We
can also estimate the reliability of a state estimator by computing the
covariance of measurement errors produced from multiple lidar measurements.

Meanwhile, IMUs are low cost sensors which can be used to estimate a motion
prior for lidar odometry by integrating the rotation rate and accelerometry
measurements. However, IMUs suffer from bias instability and are susceptible
to noise. An array of multiple IMUs (MIMUs) could provide enhanced signal
accuracy with bias and noise compensation as well as increasing operational
robustness to sensor dropouts. Multi-lidar odometry [6, 7] and MIMU odometry
[8, 9] have been studied separately; however to the best of our knowledge the
fusion of the combined set has not yet been explored.

Figure 1: Illustration of the four lidars with their embedded IMUs positioned
around the data collection vehicle. The vehicle base frame $\mathtt{{B}}$, is
located at the center of the rear axle. The sensor frames of the lidars are:
$\mathtt{L}_{{F_{R}}}$, $\mathtt{L}_{{F_{L}}}$, $\mathtt{L}_{{R_{R}}}$ and
$\mathtt{L}_{{R_{L}}}$, whereas, the sensor frames of the IMUs are:
$\mathtt{{I}}_{{F_{R}}}$, $\mathtt{{I}}_{{F_{L}}}$, $\mathtt{{I}}_{{R_{R}}}$
and $\mathtt{{I}}_{{R_{L}}}$. ${F_{R}}$: front-right, ${F_{L}}$: front-left,
${R_{R}}$: rear-right and ${R_{L}}$: rear-left.

###  I-A Motivation

As our vehicles are equipped with a GNSS system, we perform the state
estimation in global frame and take advantage of it to limit the drift rate.
Since GNSS information is unreliable in urban environments (‘urban canyon’) or
in underground scenarios (such as mining) we also need to fuse information
from onboard sensors (which might be susceptible to signal loss) to create
robust state estimates.

To achieve this, we have identified two broad problems: Problem 1: State
estimation over long time horizons with onboard sensing inherently suffers
from drift. Problem 2: Sensor signals are susceptible to loss, due to
networking or operational issues, which might affect reliability of the state
estimates using onboard sensing.

In our work we have addressed Problem 1) by fusing the onboard state estimator
with GNSS-based estimates. For Problem 2) we provide our formulation for the
fusion of multiple lidars and MIMU to create robust signals under noisy
conditions, which can give robustness both to failure of an individual lidar
or IMU sensor as well as expanding the lidar field-of-view.

###  I-B Contribution

Our work is motivated by the broad literature in lidar and IMU based state
estimation. Our proposed contributions are:

  * •

Multi-lidar odometry using a single fused local submap and handling lidar
dropout scenarios.

  * •

MIMU fusion which compensates for the Coriolis effect and accounts for
potential signal loss.

  * •

A factor graph framework to jointly fuse multiple lidars, MIMUs and GNSS
signals for robust state estimation in a global frame.

  * •

Experimental results and verification using data collected from Scania
vehicles with the sensor setup shown in Fig. 2, with FoV schematics similar to
Fig. 1.

##  II Related Work

There have been multiple studies of lidar and IMU-based SLAM after the seminal
LOAM paper by Zhang et al. [10] which is itself motivated by generalized-ICP
work by Segal et al. [11]. In our discussion we briefly review the relevant
literature.

###  II-A Direct tightly coupled multi-lidar odometry

To develop real-time SLAM systems simple edge and plane features are often
extracted from the point clouds and tracked between frames for computational
efficiency [3, 12, 13]. Using IMU propagation, motion priors can then be used
to enable matching of point cloud features between key-frames. However, this
principle cannot be applied to featureless environments. Hence, instead of
feature engineering, the whole point cloud is often processed which has an
analogy to processing the whole image in visual odometry methods such as LSD-
SLAM, [14], which is known as direct estimation. To support direct methods,
recently Xu et al. [15] proposed the ikd-tree in their Fast-LIO2 work which
efficiently inserts, updates, searches and filters to maintain a local submap.
The ikd-tree achieves its efficiency through “lazy delete” and “parallel tree
re-balancing” mechanisms.

Furthermore, instead of point-wise operations the authors of Faster-LIO [16]
proposed voxel-wise operations for point cloud association across frames and
reported improved efficiency. In our work we also maintain an ikd-tree of the
fused lidar measurements and tightly couple the relative lidar poses, IMU
preintegration and GNSS prior in our propose estimator. Since, we jointly
estimate the state based on the residual cost function built upon the multiple
modalities we consider this a tightly coupled system.

While there have been many studies of state estimation using single lidar and
IMU, there is limited literature available for fusing multi-lidar and MIMU
systems. Our idea closely resembles M-LOAM [7], where state estimation with
multiple lidars and calibration was performed. However, the working principles
are different as M-LOAM is not a direct method and MIMU system is not
considered in their work.

Finally, most authors do not address how to achieve reliability in situations
of signal loss — an issue which is important for practical operational
scenarios.

##  III Problem Statement

###  III-A Sensor platform and reference frames

Figure 2: Reference frames conventions for our vehicle platform. The world
frame $\mathtt{W}$ is a fixed frame, while the base frame $\mathtt{{B}}$, as
shown in Fig. 1, is located at the rear axle center of the vehicle. Each
sensor unit contains the two optical frames $\mathtt{C}$, an IMU frame,
$\mathtt{I}$, and lidar frame $\mathtt{L}$. The cameras are shown for
illustration only and they are not used for this work.

The sensor platform with its corresponding reference frames is shown in Fig. 2
along with the illustrative sensor fields-of-view in Fig. 1. Each of the
sensor housings contain a lidar with a corresponding embedded IMU and two
cameras. Although we do not use the cameras in this work they are illustrated
here to show the full sensor setup. We used logs from a bus and a truck with
similar sensor housings for our experiments. The two lower mounted modules
from the rear, present in both the vehicles, are not shown here in the
picture. The embedded IMUs within the lidar sensors are used to form the MIMU
setup.

Now we describe the necessary notation and reference frames used in our system
according to the convention of Furgale [17]. The vehicle base frame,
$\mathtt{{B}}$ is located on the center of the rear-axle of the vehicle.
Sensor readings from GNSS, lidars, cameras and IMUs are represented in their
respective sensor frames as $\mathtt{G}$, $\mathtt{L}^{(k)}$,
$\mathtt{C}^{(k)}$ and $\mathtt{I}^{(k)}$ respectively. Here,
$k\in[F_{L},F_{R},R_{L},R_{R}]$ denotes the location of the sensor in the
vehicle corresponding to front-left, front-right, rear-left and rear-right
respectively. The GNSS measurements are reported in world fixed frame,
$\mathtt{W}$ and transformed to $\mathtt{{B}}$ frame by performing a
calibration routine outside the scope of this work. In our discussions the
transformation matrix is denoted as,
$\mathbf{T}=\left[\begin{array}[]{ll}\mathbf{R}_{3\times
3}&\mathbf{t}_{3\times 1}\\\
\mathbf{0}^{\top}&1\end{array}\right]\in\mathrm{SE}(3)$ and
$\mathbf{R}\mathbf{R}^{T}=\mathbf{I}_{3\times 3}$, since the rotation matrix
is orthogonal.

###  III-B Problem formulation

Our primary goal is to estimate the position
$\tensor[_{\mathtt{{W}}}]{\mathbf{t}}{{}_{\mathtt{{W}}\mathtt{{B}}}}$,
orientation
$\tensor[_{\mathtt{{W}}}]{\mathbf{R}}{{}_{\mathtt{{W}}\mathtt{{B}}}}$, linear
velocity $\tensor[_{\mathtt{W}}]{\mathbf{v}}{{}_{\mathtt{{W}}\mathtt{{B}}}}$,
and angular velocity
$\tensor[_{\mathtt{W}}]{\boldsymbol{\omega}}{{}_{\mathtt{{W}}\mathtt{{B}}}}$,
of the base frame $\mathtt{{B}}$, relative to the fixed world frame
$\mathtt{W}$. Additionally, we also estimate the MIMU biases
$\tensor[_{\mathtt{{B}}}]{\mathbf{b}}{{}^{g}},\;\tensor[_{\mathtt{{B}}}]{\mathbf{b}}{{}^{a}}$
expressed in $\mathtt{{B}}$ frame, as that is where it can be sensed. Hence,
our estimate of vehicle’s state $\boldsymbol{x}_{i}$ at time $t_{i}$, is
denoted as:

$\displaystyle\boldsymbol{x}_{i}=\left[{\mathbf{R}}_{i},\;{\mathbf{t}}_{i},\;{\mathbf{v}}_{i},\;{\boldsymbol{\omega}}_{i},\;\mathbf{b}^{a}_{i},\;\mathbf{b}^{g}_{i}\right]\in\mathrm{SE}(3)\times\mathbb{R}^{15},$
(1)

where, the corresponding measurements are in the frames mentioned above.

##  IV Methodology

###  IV-A Initialization

To provide an initial pose we use the GNSS measurements,
$\tensor[_{\mathtt{W}}]{\mathbf{T}}{{}_{\mathtt{{\mathtt{W}\mathtt{G}}}}}$ and
determine an initial estimate of the starting yaw and translation component,
$\tensor[_{\mathtt{W}}]{\mathbf{t}}{{}_{\mathtt{{\mathtt{W}\mathtt{{B}}}}}}$,
in UTM co-ordinates. The GNSS unit also provides IMU measurements in
$\mathtt{G}$ frame using its internally embedded IMU, which is transformed to
$\mathtt{{B}}$ frame. For gravity alignment we use the equations [18, Eq. 25,
26] to estimate roll and pitch after collecting IMU data when the vehicle is
static for $t$ seconds according to the GNSS estimate. Subsequently, only the
translation component
$\tensor[_{\mathtt{W}}]{\mathbf{t}}{{}_{\mathtt{W}\mathtt{G}}}$, from the GNSS
measurements are used.

Meanwhile, the noise processes and starting bias estimates for the embedded
IMU sensors were characterized in advance by estimating the Allan Variance
[19] parameters using logs collected while stationary.

###  IV-B Handling lossy signal scenario

To synchronizate the different signals and to handling signal drop we maintain
individual buffer queues for each of the lidar and IMU signals. We read the
data from the buffers in a FIFO approach and compare the header timestamps
(within a threshold of 10 $\mathrm{m}\mathrm{s}$ for lidar and 1
$\mathrm{m}\mathrm{s}$ for IMU) to associate corresponding signals across
different sensors as illustrated in Fig. 3. We add an age penalty for unused
messages in the buffer queues to avoid unbounded queue growth. Unlike the
default approximate time sync policy (in ROS) we ensure that data from a
single sensor is returned if the other corresponding signals are available. As
an illustration, in Fig. 3, let us assume we receive lossy lidar and IMU
signals from our sensor setup. The green cells represent signal availability
in the buffer queue, whereas the white cells represent signal loss. The
signals are read in a FIFO approach and fused after corresponding associations
are formed based on timestamps, which is shown with links between the messages
in the middle part of Fig. 3. The right hand side illustrates the components
of the fused signal. For example at $t_{2}$, the fused signals comprises of
$F_{L}$ and $R_{R}$ signal only.

Figure 3: Handling and fusing lossy signals. Green cells denote available data
while white cells denote a lost signal for the corresponding timestamp in the
buffer queue.

###  IV-C Lidar fusion

After using the prior extrinsic calibration and signal correspondences we fuse
the lidar data from the individual sensors in the $\mathtt{{B}}$ frame using a
ikd-tree [15] and downsample them to a 5 $\mathrm{c}\mathrm{m}$ voxel
resolution. To compensate for the motion distortion during a lidar scan, all
points are projected to a pose corresponding to the beginning of the
corresponding scan using IMU pre-integration [20] between the start and end of
the lidar scan. We use the dual quaternion (DQ) operator to handle
interpolation of translation and rotation components together using screw
linear interpolation (ScLERP) [21]. Let, $\mathbf{p}_{t}$ be the the timestamp
of a received lidar point at time $t\in[t_{i},t_{i+1}]$. The point can be
transformed to the beginning of the scan as
$\mathbf{T}_{\mathtt{L}_{i}\mathtt{L}_{t}}^{-1}\mathbf{p}_{t}$, where,
$\mathbf{T}_{\mathtt{L}_{i}\mathtt{L}_{t}}$ can be recovered in DQ form as
($\mathbf{Q}$ is the DQ form of $\mathbf{T}$),

$\displaystyle\mathbf{Q}_{\mathtt{L}_{i}\mathtt{L}_{t}}=\mathbf{Q}_{\mathtt{L}_{i}}\otimes(\mathbf{Q}_{\mathtt{L}_{i}}^{-1}\otimes\mathbf{Q}_{\mathtt{L}_{i+1}})^{\eta},\quad\eta=\frac{t-t_{i}}{t_{i+1}-t_{i}}.$
(2)

###  IV-D IMU fusion

To create a robust IMU signal we fuse the available embedded IMUs as a set,
which we call a MIMU. A naive approach for MIMU signal fusion would be
averaging. Instead we use a MLE approach for the fusion [8].

Figure 4: Illustration of IMU transformation with Coriolis force compensation.

####  IV-D1 IMU transformation with Coriolis compensation

Let,
$\left[\mathbf{R}_{\mathtt{{B}}\mathtt{I}^{(k)}},\mathbf{t}_{\mathtt{{B}}\mathtt{I}^{(k)}}\right]$
be the calibration parameters between the $\mathtt{{B}}$ and the $k^{th}$ IMU
frame, $\mathtt{I}^{(k)}$.
$\mathbf{f}_{\mathtt{I}^{(k)}},\boldsymbol{\omega}_{\mathtt{I}^{(k)}}$ are the
linear acceleration and angular velocity measurements for the $k^{th}$ IMU in
its corresponding sensor frame, $\mathtt{I}^{(k)}$. Before the fusion of the
MIMU signals, all the IMU measurements must be converted to the $\mathtt{{B}}$
frame by compensating for the Coriolis force illustrated in Fig. 4. Please
note we have omitted $k$ superscript for brevity in rest of the discussion.

$\mathbf{R}_{BI}^{-1}\mathbf{f}_{I}=\mathbf{f}_{B}+\underbrace{\boldsymbol{\omega}\times\left(\boldsymbol{\omega}\times\mathbf{t}_{BI}\right)}_{\text{Centrifugal
force
}}+\underbrace{\dot{\boldsymbol{\omega}}\times\mathbf{t}_{BI}}_{\text{Euler
force }}$ (3) $\displaystyle\mathbf{R}_{BI}^{-1}\mathbf{f}_{I}$
$\displaystyle=\mathbf{f}_{B}+[\boldsymbol{\omega}]_{\times}[\boldsymbol{\omega}]_{\times}\mathbf{t}_{BI}+[\dot{\boldsymbol{\omega}}]_{\times}\mathbf{t}_{BI}$
(4)
$\displaystyle=\mathbf{f}_{B}+[\boldsymbol{\omega}]^{2}_{\times}\mathbf{t}_{BI}-[\mathbf{t}_{BI}]_{\times}\dot{\boldsymbol{\omega}},$

where, $[\boldsymbol{.}]_{\times}$ is a skew-symmetric matrix and
$[\mathbf{a}]_{\times}\mathbf{b}=-[\mathbf{b}]_{\times}\mathbf{a}$. In the
first step of the IMU fusion, we orient all the IMU measurements in their
respective sensor frames to $\mathtt{{B}}$ frame by applying eq. 4 with
suitable calibration and noise parameters, estimated in advance while
characterizing the IMUs.

####  IV-D2 Accelerometer signal model

From eq. 4 we can also derive the following signal model for an array of $K$
accelerometers, where $\mathbf{y}_{f}$ is the acceleration measurement vector
and $\boldsymbol{\eta}^{f}$ is the measurement noise of the accelerometers:

$\mathbf{y}_{f}=\mathbf{h}_{f}(\boldsymbol{\omega})+\mathbf{H}\Phi+\boldsymbol{\eta}^{f},$
(5)

where,

$\mathbf{h}_{f}(\boldsymbol{\omega})=\begin{bmatrix}{[\boldsymbol{\omega}]^{2}_{\times}}\mathbf{t}_{BI^{(1)}}\\\
\vdots\\\
[\boldsymbol{\omega}^{2}_{\times}]\mathbf{t}_{BI^{(K)}}\end{bmatrix},\mathbf{H}=\begin{bmatrix}-[{\mathbf{t}}_{BI^{(1)}}]_{\times}&\mathbf{I}_{3}\\\
\vdots&\vdots\\\
-[{\mathbf{t}}_{BI^{(K)}}]_{\times}&\mathbf{I}_{3}\end{bmatrix},\Phi=\begin{bmatrix}\dot{\boldsymbol{\omega}}\\\
\mathbf{f}_{B}\end{bmatrix}.$

(6)

####  IV-D3 Gyroscope signal model

The angular velocity is sensed by gyroscope on a rigid body is the same no
matter where the sensor is located, the same applies for angular acceleration.
Hence we can derive a simple model for the measurement vector
$\mathbf{y}_{\omega}$, given the gyro measurement noise
$\boldsymbol{\eta}^{\omega}$:

$\displaystyle\mathbf{y}_{\omega}$
$\displaystyle=\mathbf{h}_{\omega}(\boldsymbol{\omega})+\boldsymbol{\eta}^{\omega},$
(7)

where,
$\mathbf{h}_{\omega}(\boldsymbol{\omega})=\mathbf{1}_{K}\otimes\boldsymbol{\omega}$.
Note that in eq. 7 we do not model the bias, as it has been done in IMU pre-
integration. Here $\mathbf{1_{K}}$ is a column vector of ones and $\otimes$ is
the kronecker product.

####  IV-D4 MIMU signal model

Concatenation of eq. 5 and eq. 7 yields:

$\mathbf{y}=\mathbf{h}(\boldsymbol{\omega})+\mathbf{H\Phi}+\boldsymbol{\eta},$
(8)

where,

$\mathbf{y}=\begin{bmatrix}\mathbf{y}_{f}\\\
\mathbf{y}_{\omega}\end{bmatrix},\mathbf{h}(\boldsymbol{\omega})=\begin{bmatrix}\mathbf{h}_{f}(\boldsymbol{\omega})\\\
\mathbf{h}_{\omega}(\boldsymbol{\omega})\end{bmatrix},\mathbf{H}=\begin{bmatrix}\mathbf{H}_{f}\\\
\mathbf{0}_{3K,6}\end{bmatrix},\boldsymbol{\eta}=\begin{bmatrix}\boldsymbol{\eta}^{f}\\\
\boldsymbol{\eta}^{\omega}\end{bmatrix},$

(9)

and $\mathbf{0}_{m,n}$ denotes a zero matrix of size $m$ by $n$.

####  IV-D5 Maximum likelihood estimator

The objective is to use eq. 8 to form a log-likelihood function with the
assumption of zero-mean Gaussian distribution of measurement noises:

$L(\boldsymbol{\omega},\Phi)=\frac{1}{2}\left\|\mathbf{y}-\mathbf{h}(\boldsymbol{\omega})-\mathbf{H\Phi}\right\|_{\mathbf{Q}^{-1}}^{2}.$
(10)

Here, $\mathbf{Q}$ is the covariance matrix of the system and the norm is
given by
$\left\|\mathbf{X}\right\|_{\Sigma}^{2}=\mathbf{X}^{T}\Sigma\mathbf{X}$. The
maximum likelihood estimator is then defined as:

$\displaystyle{\begin{Bmatrix}\boldsymbol{\omega}^{\star},\Phi^{\star}\end{Bmatrix}}$
$\displaystyle=\underset{\boldsymbol{\omega},\Phi}{\arg\max}\left[L\left(\boldsymbol{\omega},\boldsymbol{\Phi}\right)\right].$
(11)

To solve eq. 11, a trivial way is to fix the parameter for
$\boldsymbol{\omega}$, and to maximize the likelihood function with respect to
$\Phi$, and to then replace the solution in the function. As the angular
velocity over the sensor array is the same and independent of the location, we
can easily calculate the fused $\boldsymbol{\omega}^{\star}$ as weighted least
square of gyro measurements as:

$\displaystyle\leavevmode\resizebox{381.5877pt}{}{$\boldsymbol{\omega}^{\star}=\left(\left(\mathbf{1}_{K}^{T}\otimes\mathbf{I}_{3}\right)\mathbf{Q}^{-1}\left(\mathbf{1}_{K}\otimes\mathbf{I}_{3}\right)\right)^{-1}\left(\mathbf{1}_{K}^{\top}\otimes\mathbf{I}_{3}\right)\mathbf{Q}^{-1}\mathbf{y}_{\boldsymbol{\omega}}$}.$
(12)

Using the estimated $\boldsymbol{\omega}^{\star}$, we can now solve eq. 11 as
follows:

$\displaystyle\widehat{\Phi}\left(\boldsymbol{\omega}^{\star}\right)$
$\displaystyle=\underset{\Phi}{\arg\max}\left[L\left(\boldsymbol{\omega}^{\star},\Phi\right)\right]$
$\displaystyle\stackrel{{\scriptstyle
eq.\ref{eq:weighted_LS}}}{{=}}\left(\mathbf{H}^{\top}\mathbf{Q}^{-1}\mathbf{H}\right)^{-1}\mathbf{H}^{\top}\mathbf{Q}^{-1}\left(\mathbf{y}-\mathbf{h}(\boldsymbol{\omega}^{\star})\right).$
(13)

Note that the right hand term represents the error, i.e. measurements vector
$\mathbf{y}$ subtracted from the measurement model
$\mathbf{h}(\boldsymbol{\omega}^{\star})$.

###  IV-E Factor graph formulation

A factor graph is a type of probabilistic graphical model for Bayesian
inference, which can exploit the conditional independence between variables in
the model and has been used widely for estimation problems in robotics [22].
In our setup we optimize the intra-modal fused measurements in between
keyframe at $t_{i}$ and $t_{j}$ for state estimation. Let, ${\cal I}_{ij}$,
${\cal L}_{ij}$ and ${\cal G}_{ij}$ be the IMU, lidar and GNSS measurements
between $t_{i}$ and $t_{j}$ and $\mathcal{K}$ be the set of all keyframes. If
the measurements are assumed to be conditionally independent then the current
state at keyframe $kf$, can be estimated as,

$\boldsymbol{x}_{kf}^{\star}=\operatorname*{arg\,min}_{\boldsymbol{x}_{kf}}\|\mathbf{r}_{0}\|^{2}_{\Sigma_{0}}+\sum_{(i,j)\in\mathcal{K}_{kf}}\left(\|\mathbf{r}_{{\cal
I}_{ij}}\|^{2}_{\Sigma_{{\cal I}_{ij}}}+\right.\\\ \left.\|\mathbf{r}_{{\cal
L}_{ij}}\|^{2}_{\Sigma_{{\cal L}_{ij}}}+\|\mathbf{r}_{{\cal
G}_{i}}\|^{2}_{\Sigma_{{\cal G}_{i}}}\right),$ (14)

where, each term is the squared residual error associated to a factor type,
weighted by the inverse of its covariance matrix. The notations:
$\mathbf{r}_{0}$, $\mathbf{r}_{{\cal I}_{ij}}$, $\mathbf{r}_{{\cal L}_{ij}}$,
$\mathbf{r}_{{\cal G}_{i}}$ represent residuals for state prior, pre-
integrated IMU, lidar odometry and GNSS factors respectively, described in the
following sections.

Figure 5: Factor graph structure with preintegrated IMU and lidar factors.
GNSS poses are also added as unary factors in the graph.

####  IV-E1 Prior factor

Prior factors are used to constrain the system with the $\mathtt{{W}}$ frame.
Assuming known prior
$\tensor[_{\mathtt{W}}]{\mathbf{T}}{{}_{\mathtt{{GNSS}}_{0}}}$, we can anchor
our system to the $\mathtt{W}$ frame by computing the residual between the
estimated state $\widehat{\boldsymbol{x}}_{0}$, and GNSS based state estimate
$\tensor[_{\mathtt{W}}]{\mathbf{T}}{{}_{\mathtt{{GNSS}}_{0}}}$, at $t_{0}$.
The residual becomes:

$\mathbf{r}_{0}=\left(\begin{array}[]{c}\check{\Pi}\left(\widehat{\mathbf{T}}_{0}^{-1}\mathbf{T}_{\mathtt{{GNSS}}_{0}}\right)\\\
\widehat{\mathbf{v}}_{0}-[0\;0\;0]^{\mathtt{T}}\\\
\widehat{\mathbf{\boldsymbol{\omega}}}_{0}-[0\;0\;0]^{\mathtt{T}}\\\
\widehat{\mathbf{b}}_{0}^{a}-\mathbf{b}_{{0}}^{a}\\\
\widehat{\mathbf{b}}_{0}^{\omega}-\mathbf{b}_{{0}}^{\omega}\end{array}\right),$
(15)

where, $\check{\Pi}$ is the lifting operator used to solve the optimization in
tangent space according to the “lift-solve-retract” principle presented in
[20].

####  IV-E2 IMU pre-integration factor

We use the standard formulation for IMU measurement pre-integration [20] to
create the IMU factor between consecutive nodes of the pose graph. The high
frequency pose update (at about 100Hz) is also used by the lidar odometry
module. The residual has the form (with a derivation taken from [20]):

$\mathbf{r}_{{\cal
I}_{ij}}=\left[\mathbf{r}^{\mathsf{T}}_{\Delta\mathbf{R}_{ij}},\mathbf{r}^{\mathsf{T}}_{\Delta\mathbf{t}_{ij}},\mathbf{r}^{\mathsf{T}}_{\Delta\mathbf{v}_{ij}},\mathbf{r}_{\mathbf{b}^{a}_{ij}},\mathbf{r}_{\mathbf{b}^{g}_{ij}}\right].$
(16)

####  IV-E3 Lidar odometry and between factor

##### Lidar odometry

For the lidar odometry we register the fused downsampled point cloud directly
to local submap using ikd-tree [15] leveraging its computational efficiency
when finding correspondences. ICP is performed between corresponding points to
infer relative motion, using a motion prior from IMU pre-integration. We
register lidar poses at ICP frequency ($\approx$ 2 Hz), so as to used to
maintain a smooth local submap. Since we directly register the raw point
clouds to the local submap without feature engineering it improves the
accuracy of the odometry system in feature-less environments. Also, as we are
not specifically doing SLAM, we update the local submap periodically. This
helps in managing a small memory footprint unlike maintaining a global map
(which would need large memory) required for the lidar odometry.

##### Between factor formulation

Based on our optimization window we may sample multiple lidar odometry
measurements between $t_{i}$ and $t_{j}$. However, to constrain the pose graph
with a single between factor within the smoothing window we use the first and
the last estimated poses from the lidar odometry module, denoted as,
$\widehat{\mathbf{T}}_{i+\epsilon}$ and $\widehat{\mathbf{T}}_{j-\epsilon}$
respectively. Note that we add the term $\epsilon$ to denote the lidar
odometry module processing time. Considering $\epsilon\to 0$, the residual has
the form:

$\mathbf{r}_{{\cal
L}_{ij}}=\check{\Pi}\left(({\mathbf{T}}_{i}^{-1}{\mathbf{T}}_{j})^{-1}\widehat{\mathbf{T}}_{i}^{-1}\widehat{\mathbf{T}}_{j}\right).$
(17)

####  IV-E4 GNSS prior factor

Lidar inertial odometry will drift over time. Hence, to constrain the state
estimation over long distances, we add only the
$\tensor[_{\mathtt{W}}]{\mathbf{t}}{{}_{\mathtt{{GNSS}}}}$ represented in UTM
co-ordinates as a prior on our factor graph formulation if the estimated
position covariance is larger than the GNSS position covariance adapted from
[13]. Note that our graph formulation will naturally produce state estimates
even when no GNSS information is available. Since, we only use the translation
components, the residual becomes:

$\mathbf{r}_{{\cal
G}_{i}}=|\tensor[_{\mathtt{W}}]{\mathbf{t}}{{}_{\mathtt{{GNSS}}_{i}}}-\tensor[_{\mathtt{W}}]{\mathbf{t}}{{}_{\mathtt{{\mathtt{{B}}}}_{i}}}|.$
(18)

For simplicity, we have used constant values for the covariances of IMU pre-
integration factor, lidar between factor and GNSS prior factor for our
experiments.

##  V Experimental Results

To demonstrate system performance, we conducted several experiments in real
world environments including urban and highway driving scenarios to collect
data using two vehicles with the sensor setup illustrated in Fig. 2.

###  V-A Dataset

All sensing data in our vehicle is synchronized using PTP (Precision Time
Protocol). For ground truth (GT) pose estimates, we use data from a GNSS
receiver.

We use the following notation to distinguish between different sensor
configurations – $\mathbf{L_{n_{l}}I_{n_{i}}G_{n_{g}}}$ (L: Lidar, I: IMU and
G: GNSS); with $\mathbf{n_{l}}$, $\mathbf{n_{i}}$ and $\mathbf{n_{g}}$:
indicating the corresponding number of sensor.

We collected 5 test sequences which is described in details in Table I. Lossy
signal conditions were emulated manually for our experiments by dropping
individual sensor signals randomly once for a period of 5-10 seconds in two
data logs. For the single modality experiments, we use the $F_{L}$ lidar and
its corresponding IMU in both the bus and the truck logs.

Datasets collection details for the experiments
---
Data | Vehicle | Scenario | Sensor setup | Length ($\mathrm{K}\mathrm{m}$) | Duration ($\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{s}$) | Signal loss
Seq-1 | Bus | Urban | ${L_{4}I_{4}G_{1}}$ | 1.643 | 333.37 | Yes
Seq-2† | Bus | Highway | ${L_{4}I_{4}G_{1}}$ | 1.948 | 168.64 | No
Seq-3† | Bus | Highway | ${L_{4}I_{4}G_{1}}$ | 1.577 | 118.56 | No
Seq-4 | Truck | Urban | ${L_{4}I_{4}G_{1}}$ | 0.453 | 78.91 | No
Seq-5 | Truck | Urban | ${L_{4}I_{4}G_{1}}$ | 1.267 | 151.37 | Yes
†Four front lidars and IMUs are used for this log sequence.

TABLE I:

We demonstrate the lidar fusion method and its results in signal dropout
scenario in the supplementary video.

###  V-B IMU fusion

All the IMUs from the integrated sensor housing (shown in Fig. 2) are oriented
along the $\mathtt{{B}}$ frame using the predefined calibration parameters.
Due to limited space availability, we only visualize and discuss the results
of the fusion of the angular velocity (Y-component) for Seq-1 and confirm that
a similar hypothesis holds true for the other IMU signal components in all the
sequences. As seen in Fig. 6, the angular velocity components look out-of-
phase, since they are not transformed to the $\mathtt{{B}}$ frame. Also,
around fifth second in Fig. 6, we observe that there was a signal loss for the
$F_{L}$ and $R_{R}$ IMU for around 10 seconds.

After applying the MIMU fusion technique, described in Sec. IV-D, with the
correct calibration parameters we get a fused MIMU signal in $\mathtt{{B}}$
frame, which is robust to signal loss as seen in Fig. 7. We also analyzed the
fused signal wrt GT and verify the quality of the fusion. The GT IMU signal is
logged from the GNSS system in $\mathtt{{B}}$ frame. As seen in Fig. 7 our
fused signal for the angular velocity (Y-component) closely corresponds to its
respective GT.

Figure 6: IMU raw signal of angular velocity (Y-component). Figure 7: IMU
fusion of angular velocity (Y-component).

RMSE of IMU Fusion wrt Ground Truth
---
Data | Average Fusion | MLE Fusion
Linear acceleration | Angular velocity | Linear acceleration | Angular velocity
Seq-1 | 2.868 | 0.004 | 1.716 | 0.003
Seq-2 | 2.219 | 0.003 | 1.814 | 0.003
Seq-3 | 2.186 | 0.003 | 1.792 | 0.003
Seq-4 | 2.041 | 0.008 | 1.818 | 0.007
Seq-5 | 2.263 | 0.008 | 1.820 | 0.007

TABLE II:

To provide more detailed analysis we summarize our RMSE error for the fused
signal wrt the GT signal and compare its performance to a fusion approach
based on averaging of the signals in Tab II. For the RMSE, we compared the
linear acceleration and the angular velocity components separately as,
$\operatorname{RMSE}=\sqrt{\frac{1}{N}\|I_{GT}-\widehat{I}\|^{2}_{F}}$.

Figure 8: Seq 1: Driving scenario with $F_{L}$ and $R_{R}$ lidar signal loss
(Bus)

Figure 9: Seq 2: Driving scenario with a roundabout (Bus)

Figure 10: Seq 4: Driving scenario with long turning radius (Truck)

Figure 11: Seq 5: Driving scenario with $F_{L}$ lidar and IMU signal loss
(Truck)

###  V-C State Estimation

The relative pose error (RPE) and absolute pose error (APE) metrics are used
to evaluate the local and global consistency respectively of the estimated
poses. RPE compares the relative pose along the estimated and the reference
trajectory between two timestamps $t_{i}$ and $t_{j}$, using the relative
transform over a specific distance traveled (10 $\mathrm{m}$) as,

$\displaystyle\operatorname{RPE}_{i,j}=({\mathbf{T}}_{GT_{i}}^{-1}{\mathbf{T}}_{GT_{j}})^{-1}(\widehat{{\mathbf{T}}}{{}_{i}^{-1}}\widehat{{\mathbf{T}}}{{}_{j}}).$
(19)

APE is based on the absolute relative pose difference between two poses
${\mathbf{T}}{{}_{GT_{i}}}$ and $\widehat{{\mathbf{T}}}{{}_{i}}$ at timestamp
$t_{i}$ and is compared for the whole trajectory of $N$ timestamps as,

$\displaystyle\operatorname{APE}_{N}=\sqrt{\frac{1}{N}\sum_{i=1}^{N}\|{\mathbf{T}}{{}_{GT_{i}}^{-1}}\widehat{\mathbf{T}}{{}_{i}}-I_{4\times
4}\|^{2}_{F}}.$ (20)

In our experiments to demonstrate the performance improvement of using
multiple sensors versus a single sensor when signal loss occurs we used the
$L_{4}I_{4}$ configuration, described in Sec. V-A. To illustrate the reduction
in drift when additionally fusing GNSS, we used the $L_{4}I_{4}G_{1}$
configuration.

As seen in Fig. 6, signal loss for few seconds occurs in Seq-1 (Fig. 11).
Fast-LIO2 drifts quickly at that point — but maintains consistency in the rest
of the trajectory. Meanwhile, the $L_{4}I_{4}$ configuration handles this
situation seamlessly using the other available sensor signals.

In the highway scenario of Seq-2 (Fig. 11), we see that again $L_{4}I_{4}$
configuration outperforms Fast-LIO2 as there are very few environmental
structures. Hence, expanding the lidar field-of-view using $L_{4}I_{4}$ gives
more constraints for the ICP matching, outperforming Fast-LIO2.

In Seq-5 (Fig. 11), we manually induced a data drop in the $F_{L}$ lidar and
IMU for 8 seconds. Fast-LIO2 again fails to recover from the situation,
whereas, $L_{4}I_{4}$ configuration performs well again.

Pose Error – RPE (Translation[$\mathrm{m}\text{/}$], Rotation
[$\mathrm{\SIUnitSymbolDegree}\text{/}$]) / APE
---
Data | $\mathbf{L_{4}I_{4}}$ | $\mathbf{L_{4}I_{4}G_{1}}$ | $\mathbf{Fast-LIO2}\textsuperscript{\textdagger}$
Seq-1 | 0.676, 1.248 / 4.992 |  0.115, 1.719 / 3.574 | 0.977, 1.458 / 11.137
Seq-2 | 0.419, 0.264 / 5.376 |  0.158, 0.804 / 4.342 | 0.673, 0.458 / 24.422
Seq-3 | 0.361, 0.228 / 3.479 |  0.221, 0.655 / 2.268 | 0.421, 0.298 / 13.770
Seq-4 | 0.835, 0.471 / 1.446 |  0.106, 0.719 / 0.785 | 0.824, 0.457 / 0.9840
Seq-5 | 1.265, 0.454 / 9.556 |  0.605, 1.479 / 1.645 | 26.815, 2.238 / 4613.7
†Fast-LIO2 $={L_{1}I_{1}}$ configuration.

TABLE III:

Table III, show RPE and APE metrics for each sequence with Fast-LIO2 acting as
the baseline. In all the scenarios, the $L_{4}I_{4}G_{1}$ configuration
outperform the other methods, as GNSS pose priors are added seamlessly in a
graph based optimization, which constrains the pose drift but suffers from
discontinuities; thus for visualization of the lidar map we use the
$L_{4}I_{4}$ lidar odometry frame. In most of the scenarios the RPE and APE
for $L_{4}I_{4}$ configuration is less than Fast-LIO2, since it encompasses a
wider field-of-view with more matching constraints for ICP. However, in Seq-4
(Fig. 11) Fast-LIO2 is performing slightly better as the vehicle is driven in
a structure rich environment with no signal dropouts. For Seq-3 results please
check our supplementary video.

##  VI Conclusion

We have presented a method for multi-lidar and MIMU state estimation in a
factor graph optimization framework and additionally demonstrated the use of
GNSS priors. The proposed method produces comparable performance to a state-
of-the-art lidar-inertial odometry system (Fast-LIO2) but outperforms it when
there are signal dropouts with an average improvement of 61% in relative
translation and 42% in relative rotational error. For future work we aim to
introduce radar as an additional modality to discard dynamic environmental
points for improving the odometry estimate.

##  VII Acknowledgements

This research has been jointly funded by the Swedish Foundation for Strategic
Research (SSF) and Scania. The authors would like to thank Ludvig af
Klinteberg for his valuable comments which improved the manuscript and Goksan
Isil for help with experimental data creation.

## References

  * [1] X. Zuo, P. Geneva, W. Lee, Y. Liu, and G. Huang, “LIC-Fusion: LiDAR-inertial-camera odometry,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2019, pp. 5848–5854.
  * [2] T. Shan, B. Englot, C. Ratti, and D. Rus, “LVI-SAM: Tightly-coupled lidar-visual-inertial odometry via smoothing and mapping,” in _IEEE International Conference on Robotics and Automation_ , 2021, pp. 5692–5698.
  * [3] D. Wisth, M. Camurri, and M. Fallon, “VILENS: Visual, inertial, lidar, and leg odometry for all-terrain legged robots,” _IEEE Transactions on Robotics_ , pp. 1–18, 2022.
  * [4] P. Sun, H. Kretzschmar, X. Dotiwalla, A. Chouard, V. Patnaik, P. Tsui, J. Guo, Y. Zhou, Y. Chai, B. Caine, _et al._ , “Scalability in perception for autonomous driving: Waymo open dataset,” in _IEEE Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 2446–2454.
  * [5] J. Geyer, Y. Kassahun, M. Mahmudi, X. Ricou, R. Durgesh, A. S. Chung, L. Hauswald, V. H. Pham, M. Mühlegg, S. Dorn, T. Fernandez, M. Jänicke, S. Mirashi, C. Savani, M. Sturm, O. Vorobiov, M. Oelker, S. Garreis, and P. Schuberth, “A2D2: Audi Autonomous Driving Dataset,” https://www.a2d2.audi, 2020.
  * [6] M. Palieri, B. Morrell, A. Thakur, K. Ebadi, J. Nash, A. Chatterjee, C. Kanellakis, L. Carlone, C. Guaragnella, and A. Agha-Mohammadi, “Locus: A multi-sensor lidar-centric solution for high-precision odometry and 3d mapping in real-time,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 421–428, 2020.
  * [7] J. Jiao, H. Ye, Y. Zhu, and M. Liu, “Robust odometry and mapping for multi-lidar systems with online extrinsic calibration,” _IEEE Transactions on Robotics_ , vol. 38, no. 1, pp. 351–371, 2021.
  * [8] I. Skog, J. O. Nilsson, P. Händel, and A. Nehorai, “Inertial sensor arrays, maximum likelihood, and cramér–rao bound,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 16, pp. 4218–4227, 2016.
  * [9] M. Faizullin and G. Ferrer, “Best axes composition: Multiple gyroscopes imu sensor fusion to reduce systematic error,” in _IEEE European Conference on Mobile Robots_ , 2021, pp. 1–7.
  * [10] J. Zhang and S. Singh, “Low-drift and real-time lidar odometry and mapping,” _Autonomous Robots_ , vol. 41, no. 2, pp. 401–416, Feb. 2017\.
  * [11] A. Segal, D. Haehnel, and S. Thrun, “Generalized-ICP,” in _Robotics: science and systems_ , vol. 2, no. 4. Seattle, WA, 2009, p. 435.
  * [12] T. Shan and B. Englot, “LeGO-LOAM: Lightweight and ground-optimized lidar odometry and mapping on variable terrain,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2018, pp. 4758–4765.
  * [13] T. Shan, B. Englot, D. Meyers, W. Wang, C. Ratti, and D. Rus, “Lio-sam: Tightly-coupled lidar inertial odometry via smoothing and mapping,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2020, pp. 5135–5142.
  * [14] J. Engel, T. Schöps, and D. Cremers, “LSD-SLAM: Large-scale direct monocular slam,” in _European Conference on computer vision_ , 2014, pp. 834–849.
  * [15] W. Xu, Y. Cai, D. He, J. Lin, and F. Zhang, “FAST-LIO2: Fast direct lidar-inertial odometry,” _IEEE Transactions on Robotics_ , vol. 38, no. 4, pp. 2053–2073, 2022.
  * [16] C. Bai, T. Xiao, Y. Chen, H. Wang, F. Zhang, and X. Gao, “Faster-LIO: Lightweight tightly coupled lidar-inertial odometry using parallel sparse incremental voxels,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 4861–4868, 2022.
  * [17] P. Furgale, “Representing robot pose: The good, the bad, and the ugly,” in _Workshop on Lessons Learned from Building Complex Systems, IEEE International Conference on Robotics and Automation_ , 2014.
  * [18] M. Pedley, “Tilt sensing using a three-axis accelerometer,” _Freescale semiconductor application note_ , vol. 1, pp. 2012–2013, 2013.
  * [19] D. W. Allan, “Statistics of atomic frequency standards,” _Proceedings of the IEEE_ , vol. 54, no. 2, pp. 221–230, 1966.
  * [20] C. Forster, L. Carlone, F. Dellaert, and D. Scaramuzza, “On-manifold preintegration for real-time visual-inertial odometry,” _IEEE Transactions on Robotics_ , vol. 33, no. 1, pp. 1–21, 2017.
  * [21] L. Kavan and J. Žára, “Spherical blend skinning: a real-time deformation of articulated models,” in _Symposium on Interactive 3D graphics and games_ , 2005, pp. 9–16.
  * [22] F. Dellaert and M. Kaess, “Factor graphs for robot perception,” _Foundations and Trends in Robotics_ , vol. 6, pp. 1–139, Aug. 2017.


#  Early-times Yang-Mills dynamics and the characterization of strongly
interacting matter with statistical learning

Matthew Heffernan  Department of Physics, McGill University, Montréal QC H3A
2T8, Canada  Charles Gale  Department of Physics, McGill University, Montréal
QC H3A 2T8, Canada  Sangyong Jeon  Department of Physics, McGill University,
Montréal QC H3A 2T8, Canada  Jean-François Paquet  Department of Physics and
Astronomy, Vanderbilt University, Nashville TN 37235  Department of
Mathematics, Vanderbilt University, Nashville TN 37235

###### Abstract

In ultrarelativistic heavy-ion collisions, a plasma of deconfined quarks and
gluons is formed within $1$ fm/c of the nuclei’s impact. The complex dynamics
of the collision before $\approx 1$ fm/c is often described with parametric
models, which affect the predictivity of calculations. In this work, we
perform a systematic analysis of LHC measurements from Pb-Pb collisions, by
combining an _ab-initio_ model of the early stage of the collisions with a
hydrodynamic model of the plasma. We obtain state-of-the-art constraints on
the shear and bulk viscosity of quark-gluon plasma. We mitigate the additional
cost of the ab-initio initial conditions by combining Bayesian model averaging
with transfer learning, allowing us to account for important theoretical
uncertainties in the hydrodynamics-to-hadron transition. We show that, despite
the apparent strong constraints on the shear viscosity, metrics that balance
the model’s predictivity with its degree of agreement with data do not prefer
a temperature-dependent specific shear viscosity over a constant value. We
validate the model by comparing with discriminating observables not used in
the calibration, finding excellent agreement.

_Introduction_ –

The fundamental degrees of freedom of Quantum Chromodynamics (QCD) – the
theory of the nuclear strong interaction – are quarks and gluons. In the
overlap region of two heavy nuclei colliding at ultrarelativistic energies, a
dense plasma of such quarks and gluons is formed [[See, forexample, ][,
andreferencestherein.]Busza:2018rrf]. The rapid expansion of this strongly-
interacting plasma can be described with relativistic viscous hydrodynamic [2,
3, 4, 5], providing a window into the transport properties of hot and dense
nuclear matter. This theoretical development brought the heavy-ion program
into a characterization phase.

Early analyses of measurements at the Relativistic Heavy Ion Collider
established that the specific shear viscosity $\eta/s$ (the coefficient of
shear viscosity over the entropy density) of quark-gluon plasma must be small,
of the order of $1/4\pi$ [6]. Since these early analyses, considerable effort
has been devoted to quantify the full temperature dependence of both shear and
bulk viscosity [2]. This is a formidable undertaking, since the quark-gluon
plasma is formed in heavy-ion collisions for less than $10^{-23}$ seconds (10
Fermi/$c$), evolving through multiple distinct stages between the impact of
the nuclei and the final measurement of particle showers by the detectors.
Furthermore, it was understood early-on that shortcomings in the theoretical
understanding of any stage of the collision translate into uncertainties in
the determination of the viscosities [7]. An important development has been
the use of Bayesian methods to quantify precisely the uncertainties on the
viscosities extracted from model-to-data comparisons. A number of Bayesian
analyses have now been performed [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
19, 20], comparing modern multistage models of heavy-ion collisions with broad
sets of measurements from the Relativistic Heavy Ion Collider and the Large
Hadron Collider.

A long-standing source of uncertainty in heavy-ion collisions is the early
stage of the collision, which provides the initial conditions for the
relativistic viscous fluid dynamic simulation. While multiple early-stage
models have been studied and compared in the literature, current Bayesian
analyses have mostly relied on a flexible parametric model of the impact of
the collision [21], sometimes augmented with more parameters [20] or
supplemented with simplified dynamical models of the early stage of the
collision [15, 17, 20]. In this work, we replace this parametric approach with
an ab-initio model of the early stage of the collision, the IP-Glasma model
[22, 23]. The foundation of the IP-Glasma model is to describe the nuclei,
before their collision, as two color glass condensates, which is a limiting
representation of the low-energy gluons in large nuclei accelerated at high
energies [24]. The collision of the two color glass condensates generates
gluon fields [25, 26] which are evolved in time with the classical Yang-Mills
equations. This model is known for its success in describing a wide range of
heavy-ion measurements with a limited set of parameters [27, 23, 2, 28, 29,
30, 31, 32].

Using these ab-initio IP-Glasma initial conditions, we show in this work that
our state-of-the-art model of heavy-ion collisions can describe a wide set of
heavy-ion measurements. We employ Bayesian techniques to constrain the
parameters in our approach, including shear and bulk specific viscosity
coefficients. We also illustrate the power of transfer learning to account for
uncertainty from modeling choices in simulations of heavy-ion collisions. This
paper is a companion paper to a larger write-up [33] that contains many of the
details omitted here for the sake of brevity.

_Model-data comparison_ – We consider a broad range of measurements from Pb-Pb
collisions with $\sqrt{s_{\textrm{NN}}}=2.76$ TeV center-of-mass energy per
nucleon pair at the Large Hadron Collider [33]: (i) the number of charged
hadrons per unit pseudorapidity $dN_{ch}/d\eta$ [34]; (ii) the number of
identified charged hadrons per unit rapidity $dN_{i}/dy$, $i\in\\{\pi,p,k\\}$
[35]; (iii) the transverse energy per unit pseudorapidity $dE_{T}/d\eta$ [36];
(iv) the mean transverse momenta of identified hadrons $\langle
p_{T}\rangle_{i}$, $i\in\\{\pi,K,p\\}$ [35]; (v) transverse momentum
fluctuations $\delta p_{T}/\langle p_{T}\rangle$ [37]; (vi) two-particle
radial Fourier coefficients [38]; (vii) four-particle transverse Fourier
coefficient $v_{2}\\{4\\}$ [39]; (viii) event plane correlators [40, 41]; (ix)
nonlinear response coefficients that quantify mixing between higher and lower-
order modes $\chi_{n,mk}$ [41]; (x) linear and non-linear flow modes [41].
These observables represent a diverse set of the available measurements that
can be calculated with reasonable statistical accuracy. Additional observables
are used to verify the model’s predictivity, as discussed below.

The Pb-Pb collisions are simulated in three distinct phases. The ab-initio IP-
Glasma model describes the first stage of the collision, from the nuclei’s
impact to the resulting production of color fields and their subsequent
initial evolution. The energy-momentum tensor of these gluons is then used as
initialization for relativistic viscous hydrodynamics [42, 43, 44]. Because we
focus on measurements made at midrapidity, we assume boost invariance along
the collision axis [45]. An equation of state based on to lattice quantum
chromodynamics calculations is used in the hydrodynamic phase [46, 47, 48].
The viscosities of the plasma are parametrized with the same flexible
parametrizations as used in Ref. [15, 49].111For the specific shear viscosity,
we use a piecewise linear function, with low temperature slope and a high
temperature slope separated by an single inflection point; the inflection
point is often expected to be a minimum, as observed in the temperature
dependence of other fluids [50], although we allow it to be either a minimum
or a maximum. The specific bulk viscosity is parametrized with a skewed Cauchy
distribution, allowing it to peak in the deconfinement region, and decrease at
different rate above and below this peak. All details are in Ref. [33].
Finally, at a chosen local energy density, the energy-momentum tensor of the
fluid is matched to hadronic momentum distributions, and relativistic hadronic
transport [51] is used for the third phase of the collision. There is
uncertainty in determining the hadronic momentum distributions at this
transition, and in consequence two different models are studied: one based on
Grad’s 14-moment approximation of the Boltzmann equation[52, 53, 54, 55, 56,
57, 58, 59], and the other using a Chapman-Enskog expansion in the relaxation
time approximation[60, 61, 62].

Two parameters are associated with the initial IP-Glasma phase: the first
relates the scale at which gluon saturation occurs to the density of color
charges, and the second is the hydrodynamic initialization time. Eight more
parameters are used to describe the temperature dependence of the viscosities
(see previous footnote), and one final parameter is used to define the energy
density for the transition between hydrodynamics and hadronic transport. The
smaller number of parameters in the first stage of the collision, at least
four fewer parameters than previous analysis, leads to a total of $11$
parameters. We will discuss later that the number of parameters can be further
reduced to $8$ by curtailing the flexibility of the shear viscosity
parametrization.

Model-to-data comparison is performed with Bayesian inference. To access the
model’s predictions for different value of viscosities and other parameters,
we build a fast surrogate using Gaussian process emulators. We define the
admissible range of the model parameters using Generalized Normal priors, and
we sample 350 sets of parameters with a Maximum Projection Latin hypercube
[63]. Model simulations are performed for all 350 parameter sets. Rather than
training a surrogate model for each observable of interest, a dimensional
reduction using principal component analysis is performed on the observables.
We then train Gaussian processes on 20 principal components, which capture
90.6% of the models information, and use these to emulate the model’s output
for all observables.

The above emulation is performed on calculations that use the Grad’s 14-moment
model to describe the transition between hydrodynamics and hadronic transport.
Because of the uncertainty in this transition, it is important to use more
than one model to compare with measurements. Repeating the entire procedure
for our second model (Chapman-Enskog) would have been inefficient, however,
given that only a single part of the model is changed. This is a general
challenge in studying theoretical uncertainties in heavy-ion collision
modelling: it is a high priority for the community to quantify model
uncertainties, but the computational costs are often prohibitive. We address
this problem with transfer learning, building upon Ref. [64]. We write the
model’s output as

$f_{\textrm{Chapman-Enskog}}(x)=\rho f_{\textrm{Grad}}(x)+\delta(x)$ (1)

where $\rho$ and $\delta$ are respectively a linear correlation term and a
discrepancy term between the source and target models. Constraints on $\rho$
and $\delta(x)$ are obtained by comparing with 50 design points evaluated with
the Chapman-Enskog model. This technique allows us to estimate, at a
substantially reduced computational cost, the uncertainty in the transition
from hydrodynamics to individual hadrons. We emphasize the potential of this
approach to study a wider range of theoretical uncertainties in the study of
heavy-ion collisions. Additional details of the implementation of transfer
learning and emulation can be found in Ref. [65].

We use a Gaussian likelihood function to compare models and measurements [65].
We evaluate the posterior $p(H|d,I)$, given by the product of the prior
$p(H,I)$ and the likelihood $p(d|H,I)$, normalized by the Bayes evidence
$p(d,I)$:

$p(H|d,I)=\frac{p(H,I)p(d|H,I)}{p(d,I)}.$ (2)

_Constraints on shear and bulk viscosity_ –

Figure 1: Bayes Model Averaged viscous posterior with shear viscosity allowed
to vary with temperature (top panels) and the Bayes Model Averaged viscous
posterior when the shear viscosity is constant with temperature (bottom
panels). Individual 90% credible intervals for the Grad and Chapman-Enskog
viscous corrections that underlie the averaging are shown in blue and red,
respectively.

The 90% credibility intervals for the shear and bulk viscosities, obtained
from the Bayesian calibration of our multistage model of heavy-ion collisions
to the LHC data sets listed in the previous section, are shown at the top of
Figure 1. Results are shown for both the Grad and the Chapman-Enskog model of
the hydrodynamics-to-hadrons transition. The results are similar for both
models. At low temperature ($T\approx 150$ MeV), large values of the specific
bulk viscosity $\zeta/s$ are disfavored; conversely, significant values of
bulk viscosity are favored between temperatures of 200 and 300 MeV. For the
specific shear viscosity $\eta/s$, the model-to-data comparison leads to
intermediate values $\approx~{}0.1$ to $0.2$ with a modest temperature
dependence. To obtain a single constraints from multiple models, we use
Bayesian model averaging [66, 49, 65], which yields the orange (dark) band
shown in Figure 1.

The weak temperature dependence of $\eta/s$ found in Figure 1 suggests that
good agreement with measurements could also be obtained using a constant value
of the specific shear viscosity. This question goes beyond temperature
dependence of $\eta/s$ and it directly related to the predictivity of the
model. A model with a temperature-dependent shear viscosity will always be
able to describe measurements as well or better than a model with a constant
$\eta/s$ (under the straightforward assumption that the prior of the
temperature-dependent $\eta/s$ includes the constant values). The more
important question is: given the additional model flexibility from the
temperature-dependent $\eta/s$, is there a statistically meaningful difference
in the model’s agreement with data? Metrics that balance model predictivity
and its degree of agreement with data are well known (see e.g. [67, 68]). In
this work, we use the Bayes factor, as in Ref. [15]. This is calculated by
taking the ratio of the Bayes evidences – the denominator $p(d,I)$ in Bayes’
theorem Eq. 2 – between two posterior distributions.

The Bayes factor between constant and temperature-dependent $\eta/s$ is $\ln
B=0.2\pm 2.4$ showing a negligible preference for temperature-dependent
$\eta/s$. This is the result with the Grad hydrodynamics-to-hadron model; the
result is very similar with the Chapman-Enskog RTA model, at $\ln B=1\pm 4$.
We note that Ref. [15] had also found similar Bayes factors with a constant
and temperature-dependent $\eta/s$.

While the degree of constraint on $\eta/s$ is significant, conclusive evidence
for temperature dependence remains elusive. In essence, the additional
complexity of the model with a temperature-dependence $\eta/s$ is not
justified by the corresponding small improvement in describing experimental
results. The posterior for the viscosities, assuming a constant $\eta/s$, are
shown in the bottom panels of Figure 1. Constraints on the bulk viscosity are
similar for both the constant and the temperature-dependent $\eta/$.

Figure 2: Postdiction of $\delta p_{T}/\langle p_{T}\rangle$ at Maximum a
Posteriori. Data from [37].

Confirming the results from the Bayes factor, we find a similar agreement with
measurements for both the constant and temperature-dependent value of $\eta/s$
shown in Figure 1. We show in Figure 2 the results for one of the measurements
used for the calibration, the transverse momentum fluctuation. The results for
a constant and a temperature-dependent $\eta/s$ are consistent within
statistical uncertainties. We verified that the same conclusion can be reached
for all measurements used in the calibration (see Ref. [65]).

Figure 3: Prediction of ATLAS event-plane correlations at Maximum a
Posteriori. Data [40] and calculations are shifted for clarity.

We now turn to comparison with measurements that were not used to constrain
the viscosities and other model parameters. We first compare with event-plane
correlations, which describe correlation between the symmetry axes of
different Fourier harmonics [40, 41]. As can be seen in Figure 3, we find good
agreement with four different event-plane correlators, with either model of
hydrodynamics-to-hadrons transition and with either a constant $\eta/s$ and a
temperature-dependent one.

Figure 4: Correlation between $v_{2}^{2}$ and $p_{T}$ at Maximum a Posteriori,
compared to data from a higher-energy collision. Note that data [69] are at
$\sqrt{s_{NN}}=5.02$ TeV while the MAP predictions are at $\sqrt{s_{NN}}=2.76$
TeV.

We further make a prediction for the correlation between momentum anisotropy
fluctuations and average transverse momentum fluctuations. This observable was
measurements for $\sqrt{s_{\textrm{NN}}}=5.02$ TeV Pb-Pb collisions by the
ALICE [69] and ATLAS [70] collaborations. We make a prediction for
$\sqrt{s_{\textrm{NN}}}=2.76$ TeV Pb-Pb collisions, which we superimpose with
the $\sqrt{s_{\textrm{NN}}}=5.02$ TeV Pb-Pb collisions by the ALICE [69]. If,
as expected, this observable is measured to have a very small
$\sqrt{s_{\textrm{NN}}}$ dependence, we expect excellent agreement of our
prediction with the measurement.

_Summary_ –

We compared a multistage model of ultrarelativistic nuclear collisions with a
large set of measurements from the Large Hadron Collider. By using an ab-
initio model for the early stage of the collision combined with Bayesian
inference techniques for the first time, we increased the model’s
predictivity, achieving excellent agreement with measurements used in the
model-to-data comparisons, as well as additional measurements not used in the
calibration. We found a relatively temperature-independent value for the shear
viscosity over entropy density ratio, and showed that a similar description of
measurements can be obtained with constant value for $\eta/s$. It is a
milestone that the resulting model, with a constant $\eta/s$ and ab-initio
initial condition model has eight parameters, approximately _half_ the number
used in other modern models of heavy-ion collisions, yet achieves a better
agreement with data.

We emphasize the importance of such a considerable reduction in the number of
parameters used in the multistage model. It is strong support for the overall
validity of the model to have an excellent agreement with data using a small
number of parameters. This in turn provides support for the constraints on the
shear and the bulk viscosity obtained in this work.

Building upon previous work [64], we further showed the benefits of transfer
learning methods to study model uncertainties, which we used to account for
hydrodynamics-to-transport uncertainty at a fraction of the cost of previous
studies. By exploiting this tool, we expect that a much broader range of
uncertainties can be explored in the near future, making possible a fuller
picture of the transport coefficients of Quantum Chromodynamics.

###### Acknowledgements.

We are happy to acknowledge useful exchanges with members of the JETSCAPE
Collaboration. In addition, we are grateful for useful discussions with S.
Bass, D. Everett, D. Liyanage, S. McDonald, N. Miro-Fortier, and M. Singh.
This work was supported in part by the Natural Sciences and Engineering
Research Council of Canada (NSERC), in part by the U.S. Department of Energy
Grant no. DE- FG-02-05ER41367, and in part by Vanderbilt University.
Computations were performed on the Narval, Béluga, Niagara, Graham, and Cedar
supercomputers managed by the Digital Research Alliance of Canada and its
regional partners, Calcul Québec, Compute Ontario, SciNet, SHARCNET, and
WestGrid.

## References

  * Busza  _et al._ [2018] W. Busza, K. Rajagopal, and W. van der Schee, Heavy Ion Collisions: The Big Picture, and the Big Questions, Ann. Rev. Nucl. Part. Sci. 68, 339 (2018), arXiv:1802.04801 [hep-ph] .
  * Gale  _et al._ [2013a] C. Gale, S. Jeon, and B. Schenke, Hydrodynamic Modeling of Heavy-Ion Collisions, Int. J. Mod. Phys. A28, 1340011 (2013a).
  * Heinz and Snellings [2013] U. Heinz and R. Snellings, Collective flow and viscosity in relativistic heavy-ion collisions, Ann. Rev. Nucl. Part. Sci. 63, 123 (2013), arXiv:1301.2826 [nucl-th] .
  * Derradi de Souza  _et al._ [2016] R. Derradi de Souza, T. Koide, and T. Kodama, Hydrodynamic Approaches in Relativistic Heavy Ion Reactions, Prog. Part. Nucl. Phys. 86, 35 (2016), arXiv:1506.03863 [nucl-th] .
  * Romatschke and Romatschke [2019] P. Romatschke and U. Romatschke,  _Relativistic Fluid Dynamics In and Out of Equilibrium_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2019) arXiv:1712.05815 [nucl-th] .
  * Romatschke and Romatschke [2007] P. Romatschke and U. Romatschke, Viscosity Information from Relativistic Nuclear Collisions: How Perfect is the Fluid Observed at RHIC?, Phys. Rev. Lett. 99, 172301 (2007), arXiv:0706.1522 [nucl-th] .
  * Luzum and Romatschke [2008] M. Luzum and P. Romatschke, Conformal Relativistic Viscous Hydrodynamics: Applications to RHIC results at s(NN)**(1/2) = 200-GeV, Phys. Rev. C 78, 034915 (2008), [Erratum: Phys.Rev.C 79, 039903(E) (2009)], arXiv:0804.4015 [nucl-th] .
  * Petersen  _et al._ [2011] H. Petersen, C. Coleman-Smith, S. A. Bass, and R. Wolpert, Constraining the initial state granularity with bulk observables in Au+Au collisions at $\sqrt{s_{NN}}=200$ GeV, J. Phys. G 38, 045102 (2011), arXiv:1012.4629 [nucl-th] .
  * Novak  _et al._ [2014] J. Novak, K. Novak, S. Pratt, J. Vredevoogd, C. Coleman-Smith, and R. Wolpert, Determining Fundamental Properties of Matter Created in Ultrarelativistic Heavy-Ion Collisions, Phys. Rev. C 89, 034917 (2014), arXiv:1303.5769 [nucl-th] .
  * Sangaline and Pratt [2016] E. Sangaline and S. Pratt, Toward a deeper understanding of how experiments constrain the underlying physics of heavy-ion collisions, Phys. Rev. C 93, 024908 (2016), arXiv:1508.07017 [nucl-th] .
  * Bernhard  _et al._ [2015] J. E. Bernhard, P. W. Marcy, C. E. Coleman-Smith, S. Huzurbazar, R. L. Wolpert, and S. A. Bass, Quantifying properties of hot and dense QCD matter through systematic model-to-data comparison, Phys. Rev. C 91, 054910 (2015), arXiv:1502.00339 [nucl-th] .
  * Bernhard  _et al._ [2016] J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu, and U. Heinz, Applying Bayesian parameter estimation to relativistic heavy-ion collisions: simultaneous characterization of the initial state and quark-gluon plasma medium, Phys. Rev. C 94, 024907 (2016), arXiv:1605.03954 [nucl-th] .
  * Moreland  _et al._ [2020] J. S. Moreland, J. E. Bernhard, and S. A. Bass, Bayesian calibration of a hybrid nuclear collision model using p-Pb and Pb-Pb data at energies available at the CERN Large Hadron Collider, Phys. Rev. C 101, 024911 (2020), arXiv:1808.02106 [nucl-th] .
  * Bernhard  _et al._ [2019] J. E. Bernhard, J. S. Moreland, and S. A. Bass, Bayesian estimation of the specific shear and bulk viscosity of quark–gluon plasma, Nature Phys. 15, 1113 (2019).
  * Everett  _et al._ [2021a] D. Everett _et al._ (JETSCAPE), Multisystem Bayesian constraints on the transport coefficients of QCD matter, Phys. Rev. C 103, 054904 (2021a), arXiv:2011.01430 [hep-ph] .
  * Nijs  _et al._ [2021a] G. Nijs, W. van der Schee, U. Gürsoy, and R. Snellings, Transverse Momentum Differential Global Analysis of Heavy-Ion Collisions, Phys. Rev. Lett. 126, 202301 (2021a), arXiv:2010.15130 [nucl-th] .
  * Nijs  _et al._ [2021b] G. Nijs, W. van der Schee, U. Gürsoy, and R. Snellings, Bayesian analysis of heavy ion collisions with the heavy ion computational framework Trajectum, Phys. Rev. C 103, 054909 (2021b), arXiv:2010.15134 [nucl-th] .
  * Nijs and van der Schee [2022] G. Nijs and W. van der Schee, Predictions and postdictions for relativistic lead and oxygen collisions with the computational simulation code Trajectum, Phys. Rev. C 106, 044903 (2022), arXiv:2110.13153 [nucl-th] .
  * Parkkila  _et al._ [2021] J. E. Parkkila, A. Onnerstad, and D. J. Kim, Bayesian estimation of the specific shear and bulk viscosity of the quark-gluon plasma with additional flow harmonic observables, Phys. Rev. C 104, 054904 (2021), arXiv:2106.05019 [hep-ph] .
  * Nijs and van der Schee [2023] G. Nijs and W. van der Schee, A generalized hydrodynamizing initial stage for Heavy Ion Collisions, (2023), arXiv:2304.06191 [nucl-th] .
  * Moreland  _et al._ [2015] J. S. Moreland, J. E. Bernhard, and S. A. Bass, Alternative ansatz to wounded nucleon and binary collision scaling in high-energy nuclear collisions, Phys. Rev. C 92, 011901 (2015), arXiv:1412.4708 [nucl-th] .
  * Schenke  _et al._ [2012a] B. Schenke, P. Tribedy, and R. Venugopalan, Event-by-event gluon multiplicity, energy density, and eccentricities in ultrarelativistic heavy-ion collisions, Phys. Rev. C 86, 034908 (2012a), arXiv:1206.6805 [hep-ph] .
  * Schenke  _et al._ [2012b] B. Schenke, P. Tribedy, and R. Venugopalan, Fluctuating Glasma initial conditions and flow in heavy ion collisions, Phys. Rev. Lett. 108, 252301 (2012b), arXiv:1202.6646 [nucl-th] .
  * McLerran and Venugopalan [1994] L. D. McLerran and R. Venugopalan, Computing quark and gluon distribution functions for very large nuclei, Phys. Rev. D 49, 2233 (1994), arXiv:hep-ph/9309289 .
  * Bartels  _et al._ [2002] J. Bartels, K. J. Golec-Biernat, and H. Kowalski, A modification of the saturation model: DGLAP evolution, Phys. Rev. D 66, 014001 (2002), arXiv:hep-ph/0203258 .
  * Kowalski and Teaney [2003] H. Kowalski and D. Teaney, An Impact parameter dipole saturation model, Phys. Rev. D 68, 114005 (2003), arXiv:hep-ph/0304189 .
  * Gale  _et al._ [2013b] C. Gale, S. Jeon, B. Schenke, P. Tribedy, and R. Venugopalan, Event-by-event anisotropic flow in heavy-ion collisions from combined Yang-Mills and viscous fluid dynamics, Phys. Rev. Lett. 110, 012302 (2013b), arXiv:1209.6330 [nucl-th] .
  * Ryu  _et al._ [2015] S. Ryu, J. F. Paquet, C. Shen, G. S. Denicol, B. Schenke, S. Jeon, and C. Gale, Importance of the Bulk Viscosity of QCD in Ultrarelativistic Heavy-Ion Collisions, Phys. Rev. Lett. 115, 132301 (2015), arXiv:1502.01675 [nucl-th] .
  * McDonald  _et al._ [2017] S. McDonald, C. Shen, F. Fillion-Gourdeau, S. Jeon, and C. Gale, Hydrodynamic predictions for Pb+Pb collisions at 5.02 TeV, Phys. Rev. C 95, 064913 (2017), arXiv:1609.02958 [hep-ph] .
  * Ryu  _et al._ [2018] S. Ryu, J.-F. Paquet, C. Shen, G. Denicol, B. Schenke, S. Jeon, and C. Gale, Effects of bulk viscosity and hadronic rescattering in heavy ion collisions at energies available at the BNL Relativistic Heavy Ion Collider and at the CERN Large Hadron Collider, Phys. Rev. C 97, 034910 (2018), arXiv:1704.04216 [nucl-th] .
  * Schenke  _et al._ [2020] B. Schenke, C. Shen, and P. Tribedy, Running the gamut of high energy nuclear collisions, Phys. Rev. C 102, 044905 (2020), arXiv:2005.14682 [nucl-th] .
  * Gale  _et al._ [2022] C. Gale, J.-F. Paquet, B. Schenke, and C. Shen, Multimessenger heavy-ion collision physics, Phys. Rev. C 105, 014909 (2022), arXiv:2106.11216 [nucl-th] .
  * Heffernan  _et al._ [2022] M. Heffernan, C. Gale, and J.-F. Paquet, Bayesian quantification of strongly-interacting matter with IP-Glasma initial conditions, In Preparation (2022).
  * Aamodt  _et al._ [2011a] K. Aamodt _et al._ (ALICE), Centrality dependence of the charged-particle multiplicity density at mid-rapidity in Pb-Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV, Phys. Rev. Lett. 106, 032301 (2011a), arXiv:1012.1657 [nucl-ex] .
  * Abelev  _et al._ [2013] B. Abelev _et al._ (ALICE), Centrality dependence of $\pi$, K, p production in Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV, Phys. Rev. C 88, 044910 (2013), arXiv:1303.0737 [hep-ex] .
  * Adam  _et al._ [2016a] J. Adam _et al._ (ALICE), Measurement of transverse energy at midrapidity in Pb-Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV, Phys. Rev. C 94, 034903 (2016a), arXiv:1603.04775 [nucl-ex] .
  * Abelev  _et al._ [2014] B. B. Abelev _et al._ (ALICE), Event-by-event mean ${p}_{\mathbf{T}}$ fluctuations in pp and Pb-Pb collisions at the LHC, Eur. Phys. J. C 74, 3077 (2014), arXiv:1407.5530 [nucl-ex] .
  * Aamodt  _et al._ [2011b] K. Aamodt _et al._ (ALICE), Higher harmonic anisotropic flow measurements of charged particles in Pb-Pb collisions at $\sqrt{s_{NN}}$=2.76 TeV, Phys. Rev. Lett. 107, 032301 (2011b), arXiv:1105.3865 [nucl-ex] .
  * Adam  _et al._ [2016b] J. Adam _et al._ (ALICE), Anisotropic flow of charged particles in Pb-Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV, Phys. Rev. Lett. 116, 132302 (2016b), arXiv:1602.01119 [nucl-ex] .
  * Aad  _et al._ [2014] G. Aad _et al._ (ATLAS), Measurement of event-plane correlations in $\sqrt{s_{NN}}=2.76$ TeV lead-lead collisions with the ATLAS detector, Phys. Rev. C 90, 024905 (2014), arXiv:1403.0489 [hep-ex] .
  * Acharya  _et al._ [2017] S. Acharya _et al._ (ALICE), Linear and non-linear flow modes in Pb-Pb collisions at $\sqrt{s_{NN}}=$ 2.76 TeV, Phys. Lett. B 773, 68 (2017), arXiv:1705.04377 [nucl-ex] .
  * Schenke  _et al._ [2010] B. Schenke, S. Jeon, and C. Gale, (3+1)D hydrodynamic simulation of relativistic heavy-ion collisions, Phys. Rev. C82, 014903 (2010).
  * Schenke  _et al._ [2011] B. Schenke, S. Jeon, and C. Gale, Elliptic and triangular flow in event-by-event (3+1)D viscous hydrodynamics, Phys. Rev. Lett. 106, 042301 (2011), arXiv:1009.3244 [hep-ph] .
  * Paquet  _et al._ [2016] J.-F. Paquet, C. Shen, G. S. Denicol, M. Luzum, B. Schenke, S. Jeon, and C. Gale, Production of photons in relativistic heavy-ion collisions, Phys. Rev. C 93, 044906 (2016), arXiv:1509.06738 [hep-ph] .
  * Bjorken [1983] J. D. Bjorken, Highly Relativistic Nucleus-Nucleus Collisions: The Central Rapidity Region, Phys. Rev. D 27, 140 (1983).
  * Bazavov  _et al._ [2014] A. Bazavov _et al._ (HotQCD), Equation of state in ( 2+1 )-flavor QCD, Phys. Rev. D 90, 094503 (2014), arXiv:1407.6387 [hep-lat] .
  * Bernhard [2018] J. E. Bernhard,  _Bayesian parameter estimation for relativistic heavy-ion collisions_ , Ph.D. thesis, Duke U. (2018), arXiv:1804.06469 [nucl-th] .
  * [48] eosmaker, https://github.com/mrhheffernan/eos_maker/tree/smash/.
  * Everett  _et al._ [2021b] D. Everett _et al._ (JETSCAPE), Phenomenological constraints on the transport properties of QCD matter with data-driven model averaging, Phys. Rev. Lett. 126, 242301 (2021b), arXiv:2010.03928 [hep-ph] .
  * Csernai  _et al._ [2006] L. P. Csernai, J. I. Kapusta, and L. D. McLerran, On the Strongly-Interacting Low-Viscosity Matter Created in Relativistic Nuclear Collisions, Phys. Rev. Lett. 97, 152303 (2006), arXiv:nucl-th/0604032 .
  * Weil  _et al._ [2016] J. Weil, V. Steinberg, J. Staudenmaier, L. G. Pang, D. Oliinychenko, J. Mohs, M. Kretz, T. Kehrenberg, A. Goldschmidt, B. Bäuchle, J. Auvinen, M. Attems, and H. Petersen, Particle production and equilibrium properties within a new hadron transport approach for heavy-ion collisions, Phys. Rev. C 94, 054905 (2016).
  * Grad [1949] H. Grad, On the kinetic theory of rarefied gases, Commun. Pure Appl. Math 2, 331 (1949).
  * Israel [1976] W. Israel, Nonstationary irreversible thermodynamics: A Causal relativistic theory, Annals Phys. 100, 310 (1976).
  * Israel and Stewart [1979] W. Israel and J. M. Stewart, Transient relativistic thermodynamics and kinetic theory, Annals Phys. 118, 341 (1979).
  * Teaney [2003] D. Teaney, The Effects of viscosity on spectra, elliptic flow, and HBT radii, Phys. Rev. C 68, 034913 (2003), arXiv:nucl-th/0301099 .
  * Dusling  _et al._ [2010] K. Dusling, G. D. Moore, and D. Teaney, Radiative energy loss and v(2) spectra for viscous hydrodynamics, Phys. Rev. C 81, 034907 (2010), arXiv:0909.0754 [nucl-th] .
  * Monnai and Hirano [2009] A. Monnai and T. Hirano, Effects of Bulk Viscosity at Freezeout, Phys. Rev. C 80, 054906 (2009), arXiv:0903.4436 [nucl-th] .
  * Dusling and Schäfer [2012] K. Dusling and T. Schäfer, Bulk viscosity, particle spectra and flow in heavy-ion collisions, Phys. Rev. C 85, 044909 (2012), arXiv:1109.5181 [hep-ph] .
  * Denicol  _et al._ [2012] G. S. Denicol, H. Niemi, E. Molnar, and D. H. Rischke, Derivation of transient relativistic fluid dynamics from the Boltzmann equation, Phys. Rev. D 85, 114047 (2012), [Erratum: Phys.Rev.D 91, 039902(E) (2015)], arXiv:1202.4551 [nucl-th] .
  * Chapman and Cowling [1990] S. Chapman and T. G. Cowling,  _The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases_ (Cambridge university press, 1990).
  * Anderson and Witting [1974] J. Anderson and H. Witting, A relativistic relaxation-time model for the boltzmann equation, Physica 74, 466 (1974).
  * Jaiswal  _et al._ [2014] A. Jaiswal, R. Ryblewski, and M. Strickland, Transport coefficients for bulk viscous evolution in the relaxation time approximation, Phys. Rev. C 90, 044908 (2014), arXiv:1407.7231 [hep-ph] .
  * Joseph  _et al._ [2015] V. R. Joseph, E. Gul, and S. Ba, Maximum projection designs for computer experiments, Biometrika 102, 371 (2015).
  * Liyanage  _et al._ [2022] D. Liyanage, Y. Ji, D. Everett, M. Heffernan, U. Heinz, S. Mak, and J.-F. Paquet, Efficient emulation of relativistic heavy ion collisions with transfer learning, Phys. Rev. C 105, 034910 (2022), arXiv:2201.07302 [nucl-th] .
  * Heffernan  _et al._ [2023] M. R. Heffernan, C. Gale, S. Jeon, and J.-F. Paquet, Bayesian quantification of strongly-interacting matter with color glass condensate initial conditions, (2023), arXiv:2302.09478 [nucl-th] .
  * Hoeting  _et al._ [1999] J. A. Hoeting, D. Madigan, A. E. Raftery, and C. T. Volinsky, Bayesian model averaging: a tutorial, Statistical science , 382 (1999).
  * Liddle [2004] A. R. Liddle, How many cosmological parameters?, Mon. Not. Roy. Astron. Soc. 351, L49 (2004), arXiv:astro-ph/0401198 .
  * Liddle [2007] A. R. Liddle, Information criteria for astrophysical model selection, Mon. Not. Roy. Astron. Soc. 377, L74 (2007), arXiv:astro-ph/0701113 .
  * Acharya  _et al._ [2022] S. Acharya _et al._ (ALICE), Characterizing the initial conditions of heavy-ion collisions at the LHC with mean transverse momentum and anisotropic flow correlations, Phys. Lett. B 834, 137393 (2022), arXiv:2111.06106 [nucl-ex] .
  * Aad  _et al._ [2023] G. Aad _et al._ (ATLAS), Correlations between flow and transverse momentum in Xe+Xe and Pb+Pb collisions at the LHC with the ATLAS detector: A probe of the heavy-ion initial state and nuclear deformation, Phys. Rev. C 107, 054910 (2023), arXiv:2205.00039 [nucl-ex] .