# Deep learning-based surrogate model for 3-D patient-specific computational
fluid dynamics

Pan Du  Xiaozhi Zhu  Jian-Xun Wang  Aerospace and Mechanical Engineering,
University of Notre Dame, Notre Dame, IN  Applied and Computational
Mathematics and Statistics, University of Notre Dame, Notre Dame, IN
<EMAIL_ADDRESS>
###### Abstract

Optimization and uncertainty quantification have been playing an increasingly
important role in computational hemodynamics. However, existing methods based
on principled modeling and classic numerical techniques have faced significant
challenges, particularly when it comes to complex 3D patient-specific shapes
in the real world. First, it is notoriously challenging to parameterize the
input space of arbitrarily complex 3-D geometries. Second, the process often
involves massive forward simulations, which are extremely computationally
demanding or even infeasible. We propose a novel deep learning surrogate
modeling solution to address these challenges and enable rapid hemodynamic
predictions. Specifically, a statistical generative model for 3-D patient-
specific shapes is developed based on a small set of baseline patient-specific
geometries. An unsupervised shape correspondence solution is used to enable
geometric morphing and scalable shape synthesis statistically. Moreover, a
simulation routine is developed for automatic data generation by automatic
meshing, boundary setting, simulation, and post-processing. An efficient
supervised learning solution is proposed to map the geometric inputs to the
hemodynamics predictions in latent spaces. Numerical studies on aortic flows
are conducted to demonstrate the effectiveness and merit of the proposed
techniques.

###### keywords:

hemodynamics , machine learning , CFD , neural networks , cardiovascular
biomechanics

††journal: Journal of The Royal Society Interface

##  1 Introduction

Cardiovascular disease (CVD) remains the first leading cause of death and
morbidity in the U.S., posing a major healthcare concern [1]. Although medical
imaging techniques, e.g., computed tomography (CT) and magnetic resonance
imaging (MRI), have enabled the acquisition of exquisite anatomical
information and revolutionized cardiovascular medicine [2], they are usually
not able to provide hemodynamic information (e.g., flow field, pressure
losses) on their own, which is believed to be more important to CVD diagnosis
and prognosis [3]. Computational models based on physical principles of
cardiovascular systems, combined with medical imaging, enable the derivation
of functional information inaccessible by medical images alone. In the past
decade, image-based computational modeling has become a paradigm in
cardiovascular research [4, 5, 6] and is pioneering new clinical applications
[7, 8, 9]. However, the first-principle physics-based models are often
computationally expensive since they involve solving a large-scale discrete
partial differential equation (PDE) system using numerical techniques, e.g.,
finite volume or finite element methods (FVM/FEM). Particularly, when studying
dynamics of blood flows (i.e., hemodynamics), computational fluid dynamics
(CFD) or/and fluid-structure interaction (FSI) models are needed, which are
enormously time-consuming (requiring supercomputing clusters) and prone to
numerical challenges (requiring significant domain expertise). These
roadblocks have primarily limited image-based CFD modeling in clinical
applications that require timely feedback for further therapeutic assessment
and treatment planning. For example, in the case of Heartflow (a medical
company for precision heart care), image data must be sent outside the
hospital to be processed by in-house supercomputers and results sent back to
the hospital. Likewise, the high computational cost prohibits many-query
simulations for, e.g., uncertainty quantification, parameter estimation,
design optimization, etc., which are becoming increasingly vital in advancing
the utility of cardiovascular modeling to practical applications [10, 11, 12,
13, 14, 15].

As a more efficient alternative to the principled CFD/FSI models, reduced-
order modeling (ROM) has been an area of intense investigation for years. One
class of ROM is constructed by simplifying the physics of cardiovascular
systems, such as lumped parameter (LP) models, one-dimensional (1-D) models,
or hybrid 3-D/1-D/LP models [16, 17, 18], where the physical dimensions are
reduced to different levels. For example, LP models typically treat the
hemodynamic system as an electric circuit, where the relations among pressure
drop, viscous effects, and flow velocity are modeled as simplified ordinary
differential equations (ODEs). Although LP/1D models are very cheap, they only
can predict global information (integral quantities) instead of local ones
(e.g., spatiotemporal fields of velocity or wall shear stresses), while the
latter is more useful to advance cardiovascular research/healthcare. To
preserve detailed hemodynamics information, the other class of ROM is built by
projecting the full-order governing PDEs (e.g., Navier-Stokes equations) onto
a reduced subspace spanned by a group of basis functions, such as proper
orthogonal decomposition (POD) modes, known as POD-Galerkin projection, which
has been developed for facilitating hemodynamics analysis in many 3-D patient-
specific configurations [19, 20, 21, 22]. However, the projection-based ROM is
often less stable in parametric settings and highly code-intrusive, posing
great challenges to leveraging legacy CFD/FSI solvers [23, 24, 25]. With the
increasing data availability and recent advances in machine learning (ML),
there has been growing interest in developing non-intrusive data-driven
approaches for surrogate modeling of cardiovascular systems. In general, data-
driven surrogate models aim to learn the mapping between modeling inputs
(e.g., geometry, inlet/outlet, material properties, etc.) and computed outputs
(e.g., resolved flow field, pressure contour, wall shear stress distribution,
etc.) based on full-fidelity CFD/FSI simulation data. A well-trained surrogate
model can predict detailed hemodynamics information rapidly, supporting real-
time or many-query applications. For example, collocation-based polynomial
chaos expansion (PCE) has been used to build non-intrusive surrogates of 3-D
blood flow simulations for UQ tasks [11, 10]. Gao et al. [14, 15] developed a
bi-fidelity surrogate model that leverages the efficiency and accuracy of the
low- and high-fidelity CFD simulations, respectively. Deep neural network
(DNN) has become a popular surrogate modeling approach, renowned for their
universal functional approximation capability for high-dimensional system [26,
27]. Most recently, DNN-based surrogate models have been developed for
cardiovascular applications and shown great potential [28, 29, 30, 31, 32, 33,
34, 35].

However, ML-based surrogate modeling for complex patient-specific hemodynamics
still faces significant challenges that existing techniques struggle to
address simultaneously. First, the parameterization of the geometric space
spanned by irregular 3-D patient-specific shapes is very challenging. As a
result, most existing ML-based surrogate hemodynamic models focus on
inlet/outlet parameterization or studying idealized geometries, which can be
easily parameterized with a few descriptors (e.g., radius), and only a few
studies dealt with 3-D patient-specific cases [34, 35]. Second, a well-
performed ML surrogate model often requires a large number of training
samples. However, patient-specific geometries are commonly obtained by
segmentation from medical images (e.g., CT and MRA), which are often
cumbersome and time-consuming, leading to a prohibitive cost for data
generation. To address the issue, it is necessary to effectively synthesize
numerous new samples from a small set of patient-specific geometries obtained
by image segmentation.

In this paper, we try to fill the gaps and develop an ML-based surrogate
modeling framework for rapid predictions of comprehensive hemodynamics
information in 3-D patient-specific aortic geometries. In contrast to many ML-
based surrogates that predict global hemodynamics information such as pressure
loss, FFR, and mean wall shear stress (WSS), this work focuses on predicting
local blood flow information (e.g., velocity, pressure, and WSS fields), given
irregular 3-D patient-specific shapes. The contributions of the current paper
are summarized as follows. First, we developed a method to parameterize 3-D
patient-specific aorta geometries based on stochastic shape modeling (SSM).
Second, a generative model is proposed to synthesize a massive amount of 3-D
geometries from a small patient-specific dataset. Third, an automatic
simulation data generation routine is built, enabling automatic meshing,
prescribing boundary conditions, solving, and post-processing. Lastly, an
efficient ML-based surrogate is developed by learning the functional map from
geometries to hemodynamics in latent spaces, and the performance of using DNN
and bi-fidelity technique is compared. The rest of the paper is organized as
follows. The proposed framework is introduced in Section 2. Numerical results
of the surrogate modeling are presented and discussed in Section 3, and
Section 4 concludes the paper.

##  2 Methodology

The pipeline of the proposed ML-based surrogate modeling framework is
summarized by Figure 1:

Figure 1: The schematic of the proposed ML-based surrogate modeling of
patient-specific aortic flows, with several components, including dataset
generation (green box), network training (orange box), and DNN
inference/prediction (blue box).

(i) a small set of 3-D aorta geometries are reconstructed from anatomical
images acquired via CTA/MRI scans on real patients; (ii) the geometric space
spanned by the group of real patient-specific geometries is parameterized
based on statistical shape modeling (SSM), where the variation of aorta
surfaces can be described mathematically; (iii) a large number of virtual 3-D
aorta shapes are synthesized within the latent geometric space; (iv) a python-
based simulation routine is developed for automating meshing, equation
solving, and post-processing of aortic CFD simulations on the generated
ensemble of 3-D synthetic shapes, and a large dataset of hemodynamics (e.g.,
pressure, velocity, and WSS fields) are obtained; (v) The 3-D vascular shape
and hemodynamic fields of interest (FOI) of each aorta sample are encoded into
latent spaces using Principal component analysis (PCA) or proper orthogonal
decomposition (POD), and then a deep neural network (DNN) is constructed to
learn the non-linear mapping within the latent space. Once the network is
sufficiently trained, the ML surrogate model will be able to predict the fully
resolved 3-D hemodynamic flow information within a split second, given a new
aorta shape.

###  2.1 Patient-specific aorta geometry parameterization

####  2.1.1 Aorta dataset pre-processing

Eight 3-D aorta geometries are reconstructed from MRI scans of patients with
Coarctation of the aorta (COA) conditions by the German Heart Institute Berlin
[36]. Prior to the parameterization of the variation of aorta shapes, one
crucial step is to build correspondence among them, which is very challenging
since the deformation across different patients is non-isometric, especially
aortas with a complex topological structure (e.g. many branches). To reduce
complexity, this work is focused on 3D aortic geometries with homeomorphisum.
In particular, several branches of the original thoracic aorta geometries,
including left subclavian artery (LSA), left carotid artery (LCA), right
subclavian artery (RSA) and right carotid artery (RCA), are trimmed off,
simplifying their structure to a single channel. Later, a point cloud of 4000
vertices is sampled from each aorta using voronoi clustering implemented via
the pyacvd module in python, which is loosely based on Approximated Centroidal
Voronoi Diagrams (ACVD) [37]. To regularize the position and the orientation
of the aorta samples in the 3-D Cartesian coordinate system, rigid iterative
closest points (ICP) algorithm [38] is adopted to align every aorta geometry
to a template selected from the aorta dataset.

####  2.1.2 Statistical shape modeling and virtual aorta synthesis

With the aligned real aorta geometries, we establish an SSM framework to
parameterize the complex 3-D patient-specific geometric input space and
synthesize enormous virtual aorta geometries. The procedure is conducted in
three phases: (1) building correspondence; (2) projection to latent space; (3)
synthesizing virtual aortas.

Building correspondence. Let $\\{A^{i}\\}_{i=1,...,n}$ denote the set of $n$
original patient-specific aorta shapes, and each shape sample
$A^{i}\in\mathbb{R}^{3\times p}$ is represented by a point cloud of surfaces
vertices $A^{i}=\\{\mathbf{a}^{i}_{j}\\}_{j=1}^{p}$, where $p$ is the total
number of vertices. Given a template geometry
$A^{S}\in\\{A_{i}\\}_{i=1,...,n}$ selected from the set
$\\{A^{i}\\}_{i=1,...,n}$, we aim to build correspondence by deforming the
template $A^{S}$ to all the other aorta shapes $A^{T}\in\\{A_{i}\\}_{i\neq
s}$, which are referred to as source geometry and target geometries,
respectively. The diffeomprphism function
$A^{T}=\Phi(A^{S}):\mathbb{R}^{3\times p}\to\mathbb{R}^{3\times p}$,
describing the deformation from the template $A^{S}$ to each of the rest aorta
samples $\\{A_{i}\\}_{i\neq S}$, can be solved based on the controlled-point-
based large deformation diffeomorphic metric mapping (LDDMM) algorithm [39].
In LDDMM, the deformation is described as a “current” of surface, which is the
flux of a 3-D velocity vector field over the template
$\\{\mathbf{a}^{S}_{i}\\}_{i=1}^{p}$,

$\bm{v}(\mathbf{a}_{i})=\sum_{j=1}^{r}K(\mathbf{a}_{i},\mathbf{q}_{j})\cdot\bm{\mu}_{j},$
(1)

where $\\{\mathbf{q}_{j}\\}_{1}^{r}$ and $\\{\bm{\mu}_{k}\\}_{1}^{r}$
represent a set of $r$ control points and corresponding momentum vectors,
which are defined on the source geometry $A^{s}$; K is a Gaussian kernel
$K(\mathbf{a},\mathbf{q})=exp(-||\mathbf{a}-\mathbf{q}||_{L_{2}}/\sigma^{2})$
with $\sigma$ controlling the typical width of the deformation. The
deformation is performed along the pseudo time axis $\tau$, and the deformed
point cloud locations at $\tau$ are denoted by $A^{s}(\tau)$, where
$A^{s}(t=0)$ represents the initial state (i.e., the point cloud
$\\{a^{s}_{j}\\}_{j=1}^{p}$ of template $A_{s}$). The point deformations are
governed by the following ODEs,

$\dot{A}^{s}(\tau)=\mathbf{v}\big{(}A^{s}(\tau),\tau\big{)},$ (2)

which can be solved using the second-order Runge-Kutta method, and meanwhile
the control points and momentum vectors are updated based on the Hamiltonian
dynamics,

$\displaystyle\dot{\mathbf{q}}(\tau)=\mathbf{K}\Big{(}\mathbf{q}(\tau),\mathbf{q}(\tau)\Big{)}\bm{\mu}(\tau),$
(3a)
$\displaystyle\dot{\bm{\mu}}(\tau)=-\frac{1}{2}\nabla_{\mathbf{q}}{\mathbf{K}\Big{(}\mathbf{q}(\tau),\mathbf{q}(\tau)\Big{)}\bm{\mu}(\tau)^{\top}\bm{\mu}},$
(3b)

where $\mathbf{K}$ represents the kernel matrix and
$\big{(}\mathbf{K}(\mathbf{q}(\tau),\mathbf{q}(\tau))\big{)}_{i,j}=K(\mathbf{q}_{i}(\tau),\mathbf{q}_{j}(\tau))$.
We force the deformed template $\tilde{A}^{S}_{\tau^{e}}$ at the final step
(i.e., $\tau=\tau^{e}$) to match the target geometry $A^{T}$ using the
steepest gradient descent optimization algorithm. To define the loss function,
a distance metric, evaluating the difference between the deformed template
$\tilde{A}^{S}_{\tau^{e}}$ and the target geometry $A_{T}$, is required. As
the aorta geometries are composed of triangular meshes with grid centers
$\\{\mathbf{c}_{i}\\}_{i=1}^{m}$ and edge normal vectors
$\\{\mathbf{n}_{i}\\}_{i=1}^{m}$, the varifold distance metric $d$ is defined
as,

$d\bigg{(}\Big{\\{}(\mathbf{n}^{S}_{i},\mathbf{c}^{S}_{i})\Big{\\}}_{i=1}^{m^{S}},\Big{\\{}(\mathbf{n}^{T}_{j},\mathbf{c}^{T}_{j})\Big{\\}}_{j=1}^{m^{T}}\bigg{)}=\sum_{i}^{m^{S}}\sum_{j}^{m^{T}}K(\mathbf{c}^{S}_{i},\mathbf{c}^{T}_{j})\cdot\frac{\Big{(}(\mathbf{n}^{S}_{i})^{\top}\cdot\mathbf{n}^{T}_{j}\Big{)}^{2}}{||\mathbf{n}^{S}_{i}||\cdot||\mathbf{n}^{T}_{j}||}.$
(4)

The metric is calculable regardless of whether point-to-point correspondence
exists between two geometries. Namely, the total amount of grids $m^{S}$ and
$m^{T}$ for source and target geometries can be different. Once the
diffeomorphism functions are solved, all target geometries can be approximated
by the deformations from the template, establishing a correspondence for the
original shape dataset. The LDDMM algorithm has been implemented in
Deformetrica, an open-source software designed for statistical analysis of
meshes [40].

Projection to latent space. To enable a concise parameterization of the input
geometric space, PCA is applied to encode the 3-D irregular shape onto the
latent space. Specifically,
$\\{(x^{i}_{j},y^{i}_{j},z^{i}_{j})\\}_{j=1,..,p}$, denoting the Cartesian
coordinates of the point cloud on the aorta sample $\tilde{A}_{i}$, was
flattened into a vector
$\tilde{\boldsymbol{V}}^{i}=(x^{i}_{1},y^{i}_{1},z^{i}_{1},...,x^{i}_{p},y^{i}_{p},z^{i}_{p})$.
The entire 3-D shape dataset can be squeezed into a matrix
$\boldsymbol{A}=({\tilde{\boldsymbol{V}}_{1}}^{\top},...,{\tilde{\boldsymbol{V}}_{N}}^{\top})^{\top}$.
Subsequently, PCA acting on $A$ yields the eigenvectors
$\\{\boldsymbol{W}_{i}\\}_{i=1,...,z}$ and eigenvalues
$\\{\lambda_{i}\\}_{i=1,...,z}$, where $z$ is the truncation integer
controlling the number of primary components needed. Each aorta shape can be
decomposed as follow,

$\tilde{\boldsymbol{A}_{i}}\approx\boldsymbol{A}_{mean}+\sum_{j=1}^{z}\alpha^{i}_{j}\sqrt{\lambda_{j}}\boldsymbol{W}_{j},$
(5)

where the coefficient vector
$\boldsymbol{\alpha}_{i}=(\alpha^{i}_{1},...,\alpha^{i}_{z})$ represents the
encoded shape, which, in other words, is the projection of sample $A_{i}$ in
the latent space. Let $F$ and $F^{-1}$ denote the PCA operator and its
inverse, respectively. The reconstruction error $E_{pca}$ is defined to
evaluate the accuracy of the PCA transform,

$E_{pca}=\frac{||F^{-1}(F(\boldsymbol{A}))-\boldsymbol{A}||_{L2}}{||\boldsymbol{A}||_{L2}}.$
(6)

Synthesizing virtual aorta geometries. Motivated by the need of sufficient
training samples for data-driven ML surrogate modeling, we propose a shape
synthesizing method for producing a large amount of synthetic 3-D aortas from
a small set of original aorta shapes. The general idea is to sample from the
latent geometric space spanned by the original aorta shapes and then
synthesize new aortas by decoding the sampled latent shape vectors back to the
Cartesian space, as illustrated in Figure 2, where the geometric space is
spanned by three original aorta shapes.

Figure 2: Illustration of latent space distribution and visualization of
generated synthetic aortas ( ) from 3 original aorta shapes ( ), where PC1 and
PC2 represent the first and second primary components of the latent geometric
space, respectively.

To synthesize new geometries within the geometric space, we propose two
different methods: (1) random linear interpolation and (2) uniform PCA
sampling. The random linear interpolation approach linearly combines the $n$
original samples in the latent space with randomly assigned weights to
generate new samples $\boldsymbol{A}^{*}$,

$\boldsymbol{A}^{*}=F^{-1}(\sum_{i=1}^{n}\omega_{i}\boldsymbol{\alpha_{i}})=\boldsymbol{A}_{mean}+\sum_{i=1}^{n}\sum_{j=1}^{z}\omega_{i}\alpha^{i}_{j}\sqrt{\lambda_{j}}\boldsymbol{W}_{j},$
(7)

where $\boldsymbol{\omega}=[\omega_{1},\cdots,\omega_{n}]^{T}$ is a randomly
generated weight vector that satisfies $\sum_{i=1}^{z}\omega_{i}=1$. As an
alternative, uniform PCA sampling approach generates new samples by uniformly
perturbing the PCA coefficients within the ranges bounded by the original
geometries,

$\boldsymbol{A}^{*}=\boldsymbol{A}_{mean}+\sum_{j=1}^{z}\tilde{\alpha}^{*}_{j}\sqrt{\lambda_{j}}\boldsymbol{W}_{j},$
(8)

where
$\boldsymbol{\tilde{\alpha}^{*}}=[\tilde{\alpha}^{*}_{1},\cdots,\tilde{\alpha}^{*}_{z}]^{T}$
is a randomly generated vector, and each dimension $\tilde{\alpha}^{*}_{j}$ is
sampled from a uniform distribution
U$\big{[}\min{(\alpha_{j}^{i}|_{i=1}^{n})},\max{(\alpha_{j}^{i}|_{i=1}^{n})}\big{]}$.
Figure 3 plots a few generated samples within the 2-D latent space defined by
its first (PC1) and second (PC2) primary components and the corresponding
synthetic geometries. The synthetic geometries (blue) are scattered within the
triangular region bounded by three original geometries (orange).

###  2.2 Automatic CFD simulation routine

To systematically generate massive CFD simulation data, we develop a python
routine to automate tedious simulation procedures, including geometry pre-
processing, meshing, boundary setting, simulation, and post-processing. As the
first step, pre-processing of the aorta geometries is needed in order to be
suitable for CFD simulations. First, since the edges at the entry and the exit
of the raw aorta channel could be jaggy, the entry is slightly extended to
construct a circular rim at the inlet. Similarly, the exit of the aorta is
also extended, which also avoids inverse flows at the outlet. Second, the
aorta surface is smoothed to reduce bumpiness while still preserving its
geometric characteristics. Third, caps are added to seal nonphysical holes on
the aorta geometry. Once the geometries pre-processing is done, 3-D
unstructured triangular meshes are generated accordingly and then converted
into msh format suitable for subsequent CFD simulations. These operations are
implemented based on the Vascular Modeling Toolkit (VMTK), which is an open-
source module specialized in 3D reconstruction, geometric analysis and mesh
generation. A python wrapper is implemented to automate the entire procedure
in parallel for the ensemble setting.

Subsequently, a large-size ensemble of CFD simulations are performed to obtain
the corresponding FOIs for each sample of the generated dataset, based on the
open-source CFD platform, OpenFOAM. As handling a massive set of cases can
easily be cumbersome and time-consuming, a python subroutine coupled with
OpenFOAM is established, such that the procedure of each simulation, including
meshing, setting up cases, specifying boundary conditions, solving, and post-
processing, is automatically executed and repeated through all samples in the
dataset. This process can be conducted in parallel as well. During the post-
processing step, the simulated field data on CFD meshes are projected onto the
surface or volumetric mesh grids where the correspondences are established. In
particular, surface pressure and wall shear stress are projected onto surface
vertices and velocity data are projected onto the volumetric nodes via the
nearest-neighbor interpolation. The volumetric nodes are obtained by uniformly
sampling points between the centerline and the surface vertices. Finally, a
large number of input geometries and corresponding flow solutions are obtained
as the labeled dataset for surrogate modeling.

###  2.3 ML-based surrogate model

In order to build a fast forward map from the geometric space to the solution
space of interest, we propose a supervised deep learning solution based on PCA
encoding-decoding and fully-connected deep neural networks.

####  2.3.1 Encoding-decoding formulation.

Considering the significantly high dimensions of 3-D shapes and solution
fields in their discrete mesh forms, we adopt an encoding-decoding strategy to
reduce the learning complexity. Specifically, the input geometries and output
hemodynamic solution fields are encoded into low-dimensional latent spaces. To
promote training efficiency, the encoder-decoder construction is decoupled
with DNN-based surrogate training, and PCA (or POD) instead of DNN-based
autoencoder is used for the encoding-decoding process. To balance the
reconstruction accuracy and learning complexity, 100% of the PCA total
variation is preserved for the input geometries, and over 95% is preserved for
all output solution fields, resulting in reconstruction errors less than 1%.
The relationship between the encoded shape and solution fields is learned by a
multilayer perceptron (MLP), where Rectified Linear Unit (Relu) is used as the
activation function.

####  2.3.2 Learning architecture optimization.

The DNN-based forward map is learned from the labeled dataset using the Adam
optimizer, where the learning rate is adaptively changed based on the
estimates of the moments. To achieve the best learning performance, we
optimize the MLP architecture and other learning hyperparameters, including
the number of layers, the number of neurons of each layer, batch size, and
initial learning rate, using a Bayesian optimization algorithm. Specifically,
RAY-tune [41], an open-source python module for scalable hyperparameter
tuning, is adopted here to achieve this goal. The general idea is to construct
a posterior distribution of functions based on Gaussian processes within the
prescribed ranges of the hyperparameters to be optimized. During the
optimization, the posterior of hyperparameters will be updated given more and
more training data, leading to the best learning configuration. The Bayesian
optimization is superior to a brutal grid search method and can reach the
optimized state in fewer iterations. During the optimization, the Asynchronous
Successive Halving Algorithm (ASHA) is utilized to aggressively terminate non-
ideal trials in advance in the interest of saving time, and more details can
be found in [42].

####  2.3.3 Evaluation of performance

To evaluate the performance of the ML-based surrogate model, relative mean
square error (RMSE) is used to calculate the difference between prediction and
label, which is defined as,

$RMSE=\frac{\sum_{i}^{N}(x_{i}-y_{i})^{2}}{\sum_{i}^{N}{y_{i}}^{2}}\times
100\%,$ (9)

where $\boldsymbol{X}=(x_{1},...,x_{N})$ and
$\boldsymbol{Y}=(y_{1},...,y_{N})$ are predictions and labels of FOI on $N$
vertices/grids, respectively. For scalar fields such as pressure, RMSE is
calculated directly upon the pressure field. For vector fields like wall shear
stress and velocity, an array flattened from the $N\times 3$ vector-matrix
containing three directional components (x,y and z-direction) and the array of
the vector magnitude are calculated to estimate the RMSE, referred to as
vector RMSE and magnitude RMSE, respectively. When the error is evaluated upon
a test dataset with multiple samples, the average of the RMSE of each sample
is calculated.

##  3 Results and Discussion

###  3.1 Statistical shape analysis and synthesis

We first register the eight original patient-specific aorta geometries by
solving seven diffeomorphisms from a randomly chosen template to the rest of
seven geometries. An example of such diffeomorphism is demonstrated in Figure
3. The template $A_{1}$ is firstly sampled with 4000 points $P_{1}$ and so is
the target geometry $A_{2}$ with $P_{2}$. The shape deformation $A(\tau)$ is
solved with the boundary condition of $A(\tau=0)=P_{1}$ and $A(\tau=1)=P_{2}$.
Several intermediate states of the function $A$ are shown in the middle, among
which the point clouds near the left side resemble $P_{1}$ and, likewise,
point clouds near the right side resemble $P_{2}$. In addition, the
intermediate states can be used to “interpolate” between 3-D irregular,
showing geometric morphing from aorta $A_{1}$ to aorta $A_{2}$.

Figure 3: Building correspondence: Step①, sample vertices from surface; Step②,
solve the diffeomorphism function; Step③, reconstruct intermediate aorta

(a) Uniform PCA sampling

(b) Random interpolation sampling

Figure 4: Latent space distribution of synthetic ( ) and original ( ) samples
for different synthesizing method: (a) uniform PCA sampling, (b) random
interpolation sampling

Based on the set of eight original aorta geometries, a large number of virtual
aortas are synthesized by sampling the geometric space spanned by the original
shapes. Both the uniform PCA sampling and random interpolation methods are
implemented, and the sample distributions of them are plotted in Figure 4 to
compare their coverage of the input space. The 1,000 synthetic shapes are
scattered in the coordinate plane of their PCA latent space with the first
(PC1) and second (PC2) primary components. The synthetic samples obtained by
the uniform PCA sampling are much more spread out over the coordinate plane
than those from the random interpolation method, indicating a better coverage
of the input space and hence adopted for the data generation in this work.

Figure 5: Latent space distribution (top) and visualization (bottom) of
synthetic ( ) and original ( ) aorta samples

Part of the synthetic geometries and the eight original geometries are shown
in Figure 5. The uniform sampling method successfully explores the input
topological space and the synthetic geometries, though virtual, possess
geometric characteristics of biological shapes of real aortas and fully
displays the variation of the aorta shapes.

###  3.2 Surrogate CFD modeling using deep learning

Once the aorta geometry dataset is built, ensemble CFD simulations are
conducted to generate labeled hemodynamics data, which will be used to train
an ML-based surrogate.

Numerical setup. In this work, we made the following assumptions: 1) the blood
flow is Newtonian and governed by the steady-state incompressible Navier-
Stokes equations; 2) the vessel walls are rigid; 3) the inlet velocity has a
parabolic profile with a maximum magnitude of 1 m/s, which is close to the
peak blood flow velocity during one cardiac cycle in a human body. A zero-
gradient boundary condition is applied at the inlet for pressure and a
constant reference pressure is set at the outlet. The OpenFOAM, an open-source
C++ library for FVM, is used for CFD simulations. In particular, the Semi-
Implicit Method for Pressure Linked Equations (SIMPLE) algorithms [43] were
used for solving the incompressible Navier-Stokes equations, and the Rhie and
Chow interpolation with collocated grids was employed to prevent the
pressure–velocity decoupling [44]. All CFD simulations were run in parallel (4
CPUs per case), and each simulation takes around $1.2\times 10^{3}$ s to reach
the convergence. Three decoupled MLPs are constructed to predict the velocity,
surface pressure, and WSS fields separately. The network architectures and
learning parameters are optimized using RAY-tune module, which is given in
Table 2. The DNNs are trained on an Nvidia RTX A6000 GPU with a minimum of
2,000 epochs to ensure convergence.

Learning performance study. As a demonstration, we generated a dataset of
$N_{tot}=1,000$ aorta geometries and CFD simulated labels for training and
validation study. The dataset is divided into a training set and testing set
with the ratio of $r=N_{train}/(N_{train}+N_{test})$, where $N_{train}$ and
$N_{test}$ denote the size of training set and test set, respectively. The
learning-prediction performance of the trained DNN surrogates with different
training-testing ratios are reported in Table 1, where the averaged MSEs of
the DNN-predicted surface pressure, velocity, and WSS fields are compared.
Considering the randomness of the DNN initialization, the training processes
for all DNNs are repeated by ten times with different random initializations,
and the averaged prediction errors over these trials are reported for a more
rigorous assessment.

$r=\frac{N_{train}}{N_{train}+N_{test}}$ | $\epsilon_{P}$ | $\epsilon_{V}$ | $\epsilon_{VM}$ | $\epsilon_{WSS}$ | $\epsilon_{WSSM}$
---|---|---|---|---|---
r=60% | 1.83% | 0.735% | 0.504% | 1.88% | 1.33%
r=70% | 1.81% | 0.726% | 0.5% | 1.86% | 1.31%
r=80% | 1.67% | 0.734% | 0.505% | 1.86% | 1.32%
r=90% | 1.48% | 0.64% | 0.434% | 1.74% | 1.26%
Table 1: DNN prediction errors (averaged MSE) with different training to
testing size ratios for hemodynamics fields of interest, including surface
pressure ($\epsilon_{P}$), velocity vector ($\epsilon_{V}$), velocity
magnitude ($\epsilon_{VM}$), wall shear stress ($\epsilon_{W}$), wall shear
stress magnitude ($\epsilon_{WM}$).

It can be seen that the relative errors of the NN surrogate prediction are
very low once the DNNs are sufficiently trained. For velocity field
prediction, the prediction error is less than $1\%$ even only $60\%$ of the
dataset is used for training. Compared to the velocity, the errors of surface
field predictions such as surface pressure and WSS fields are slightly higher,
since they are more sensitive to the shape of vessel walls and the boundary
layer conditions, which are not easy to capture well. Also, for vector-valued
fluid quantities like velocity and WSS, the magnitude errors are smaller than
the vector errors by about $30\%$. Viewing the columns for each flow quantity,
the prediction accuracy increases as the training size grows, which is as
expected. For example, the pressure error climbs from $1.48\%$ at $r=90\%$ to
$1.83\%$ at $r=60\%$ as $r$ grows. In general, the prediction of the FOIs is
accurate, indicating the great capability of the proposed DNN-based surrogate
model.

Predicting local hemodynamic information. In contrast to traditional surrogate
model (e.g., GP models) or reduced-space model (e.g., LP/1D models), the
proposed ML-based surrogate model is able to rapidly predict local hemodynamic
information such as blood flow pattern, surface pressure and WSS distributions
over the vessel walls. For illustration, we randomly selected four test aorta
geometries, and the hemodynamic fields predicted by the DNN surrogate
($r=80\%$) are compared with the CFD reference (referred to as ground truth),
as shown in Figures 6 and 7. Besides, the absolute error contours are also
plotted out by substracting the ground truth from the prediction. For
pressures, the surrogate model yields very close surface distribution to the
ground truth, yet slight differences exist at the stenosis of the descending
aorta. From the error contours, this discrepancy is very mild on samples 1, 2
and 5, but notable on sample 3, where the severity of stenosis is much higher
than other samples. This indicates that the surrogate model is slightly less
accurate in capturing flow features with large pressure gradients, since more
flow vortices are induced by sharp changes of surface curvatures complicating
the flow physics. Surrogate prediction of the WSS distribution is in a good
agreement with the ground truth. The pattern and low/high WSS regions can be
accurately captured. However, the surrogate model tends to smooth out the
small WSS fluctuations compared to the CFD reference. It appears that high-
frequency components of the output are partially trimmed off when
interpolating among neighboring training samples based on the non-linear
relations learned by the neural network. Similar to pressure, notable
inconsistency is also observed near the stenosis on the descending aorta for
sample 3. The velocity contours on the cross-section and internal velocity
vector fields are visualized in Figure 7. Again, the predicted velocity cross-
section contours agree with the CFD ground truth very well and the overall
velocity vector fields from the neural network and CFD are almost identical
for all samples, though the error contours of planar velocity show slight
prediction errors near the stenosis region for Sample 3.

Figure 6: ML surrogate predictions of pressure (top) and wall shear stress
(bottom) are compared with the CFD ground truth: prediction (1st and 4th row),
CFD ground truth (2nd and 5th row), prediction error (3rd and 6th row) Figure
7: Cross-section velocity magnitude (top) and velocity vector filed (bottom)
comparison of prediction and CFD ground truth: prediction (1st and 4th row),
ground truth (2nd and 5th row), error (3rd row)

Training cost and prediction speedup. The mappings between shape input and
FOIs, including pressure, velocity and WSS, are learned using separated neural
networks, and each takes about two minutes to be fully trained on an Nvidia
RTX-3090 GPU. The network inference for solution field prediction takes less
than 0.0012 seconds. In contrast, a full-resolution CFD simulation needs about
20 minutes to reach the steady state, not counting the overhead for mesh
generation and case setup, which could also be time-consuming. As a result,
the fully trained DNN surrogate model has $100,000$x speedup over the full CFD
simulation in the online inference phase. We should note that a fair
comparison should also consider the offline cost of constructing the
surrogate, which includes the cost for network training and label generation.
As for network training, the cost is negligible compared to that of a CFD
simulation, thanks to the light DNN structure defined in the latent space.
However, the cost for data generation depends on the number of training labels
required, which is determined by the desired trade-off between the efficiency
and accuracy of the surrogate. For example, supposing 500 CFD labels are used
for network training, the CFD data generation cost is $500\times 20$ minutes,
which is significantly greater than that of a single CFD simulation.
Nonetheless, the surrogate model is built for many-query applications (e.g.,
optimization, uncertainty quantification, inverse modeling, etc.), and the
huge speedup of single model evaluation can be leveraged when a large number
of model evaluations are required as the offline cost is paid off at once. For
example, Markov chain Monte Carlo (MCMC), as the gold standard Bayesian
inference method for forward and inverse UQ, usually requires hundreds of
thousands of model evaluations to reach convergence, which well justifies the
data generation at the cost of hundreds of full-order CFD evaluations.
Moreover, it is worth noting that the training dataset size is application-
dependent and the data generation cost can be significantly reduced if one
allocates the computational budget wisely. A comprehensive study on the use of
the proposed surrogate model in specific UQ or optimization applications is
out of the scope of this work and will be conducted in the future.

###  3.3 Comparison of DNN-based and bi-fidelity surrogate models

Lastly, we investigate how well the proposed DNN-based surrogate model
performs compared to other state-of-art surrogate models that can provide
full-field hemodynamic predictions. To enable the full-field prediction
capability with the same resolution as that of full-order CFD simulations, a
multi-fidelity strategy is often adopted for surrogate modeling. The multi-
fidelity paradigm is pioneered by Kennedy & O’Hagan [45] decades ago,
originally from a statistical point of view (e.g., GP-based multi-model
approach). Very recently, multi-fidelity strategy has been applied for
surrogate modeling in inference and uncertainty quantification (UQ) problems
in cardiovascular biomechanics problems [46, 47, 48, 14, 15]. In particular,
Gao et al. [14] developed a bi-fidelity (BF) surrogate modeling approach for
3-D patient-specific hemodynamic simulations, which has been demonstrated
effective for both forward and inverse UQ problems [14, 15]. The theoretical
backbone of the BF surrogate is based on the multi-fidelity paradigm, which
leverages the accuracy of high-fidelity CFD solutions and the efficiency of
low-fidelity CFD solutions, largely reducing the total computation cost
compared with conventional full-order CFD simulations. Here we compared the
proposed DNN-based surrogate with the state-of-the-art BF surrogate model,
which are built with the same amount of training labels. It is worth noting
that the training overhead and prediction cost of the BF surrogate are always
higher than those of the DNN-based surrogate with the same amount of training
data (i.e., high-fidelity CFD labels), because low-fidelity CFD simulations
have to be performed for both training and inference. Specifically, to build
the BF surrogate model, in addition to high-fidelity (HF) CFD training data,
another 1,000 low-fidelity (LF) CFD simulations are conducted on low-
resolution meshes. As the low-res meshes are very coarse and the convergence
criterion is set to be large, the LF simulation cost is much lower than that
of a HF simulation. The LF data are used to explore the input geometric space
to determine important points, where HF CFD simulations are performed for
training. More details of the BF surrogate construction can be found in [14].

(a) Pressure

(b) Velocity

(c) Wall shear stress

Figure 8: Comparison of neural network(NN) and bi-fidelity model(BF): Y axis,
averaged mean square error (AMSE); X axis, trainset size (TS)

To find out how the prediction errors change with the training size, the
training set size (TS) changes from 300 to 900, and the testing set is a fixed
group of 100 aortas randomly generated. The comparison of performance is shown
in Figure 8. Both surrogate models have reasonably good predictions, but the
proposed DNN-based surrogate model outperforms the BF surrogate model for all
FOIs under all different training scenarios. Both surrogate models reach their
best accuracy at $TS=900$ for all FOIs. The prediction error of the BF model
quickly increases as the training sample size decreases, while this trend is
less notable for the DNN-based model, particularly for velocity and WSS
predictions, indicating a more robust behavior in the small data regime. For
both models, the learning performance for the 3-D velocity vector field is
slightly better than that of surface field outputs, such as pressure and WSS,
due to the fact that they are more sensitive to the local changes of the input
geometries cross difference patients. Moreover, the BF prediction cost is
higher than the DNN-based surrogate, since each BF prediction requires an LF
CFD simulation. Therefore, the proposed DNN-based surrogate model is superior
to the BF model in terms of accuracy and efficiency, though both models can
predict all FOIs reasonably well.

##  4 Conclusion

In this study, we developed a statistical shape modeling (SSM) approach to
parameterize the geometric space spanned by a set of 3-D patient aorta
geometries and systematically synthesize a large number of virtual aorta
geometries. The results show that the proposed SSM method successfully
established correspondence among aorta geometries from different patients and
the uniform PCA sampling algorithm sufficiently explored the input space,
generating abundant synthetic aorta shapes that fully cover the geometric
variation. Subsequently, an ML-based surrogate model is proposed to predict
comprehensive local hemodynamic information (e.g., surface pressure, velocity
and WSS fields), where an encoding-decoding structure is built based on PCA
and multiple MLPs are constructed to learn the non-linear mappings between the
input geometries and flow fields of interest within the latent spaces. The PCA
encoder and decoders managed to reduce the dimension of the mesh-based 3-D
aorta shapes and field outputs to a large extent and maintained a
reconstruction error as low as $1\%$. The MLPs are well optimized in terms of
their hyperparameters and trained on different sizes of sample sets. The
testing results prove that the proposed surrogate model is able to yield
hemodynamic field predictions in a good agreement with the ground truth
obtained from full-order CFD simulations, yet the accuracy is slightly lower
near the stenosis of the descending aorta, where flow features are more
complicated due to sharp pressure gradients. The accuracy of the proposed
model can be improved by increasing the size of training set. The proposed
DNN-based surrogate is compared with the state-of-the-art bi-fidelity (BF)
surrogate model, and the DNN surrogate outperforms the BF surrogate in terms
of both efficiency and accuracy for all FOIs under different training
scenarios. Furthermore, the proposed model exhibits a robust behavior when the
number of samples drops whereas the BF model easily loses its accuracy,
especially for surface pressure predictions. In conclusion, the proposed DNN-
based surrogate model shows a great capability of approximating the non-
linearity of the CFD simulations for cardiovascular flows. Considering
extremely low costs of DNN for inference, the trained DNN-based surrogate is
able to provide high-resolution local hemodynamic information in split of
second, showing a great potential for applications requires massive model
queries, such model inference, optimization, and uncertainty quantification
problems. Although the current work focus on aortic flows, the proposed
methods can be applied to hemodynamic modeling for other anatomies in general.

## Acknowledgment

The authors would like to acknowledge the funds from National Science
Foundation under award numbers CMMI-1934300 and OAC-2047127, and startup funds
from the College of Engineering at University of Notre Dame in supporting this
study.

## Compliance with Ethical Standards

Conflict of Interest: The authors declare that they have no conflict of
interest.

##  Appendix A Optimized MLP structures and hyperparameters

FOI | DNN Layers | In & Out features | batch size | learning rate
---|---|---|---|---
| Linear Layer | (7,64) |  |
| Relu | - |  |
| Linear Layer | (64,8) |  |
Pressure | Relu | - | 8 | $1\times 10^{-3}$
| Linear Layer | (8,256) |  |
| Relu | - |  |
| Linear Layer | (265,10) |  |
| Linear Layer | (7,64) |  |
| Relu | - |  |
| Linear Layer | (64,512) |  |
Velocity | Relu | - | 8 | $4.07\times 10^{-4}$
| Linear Layer | (512,256) |  |
| Relu | - |  |
| Linear Layer | (265,10) |  |
| Linear Layer | (7,512) |  |
| Relu | - |  |
| Linear Layer | (512,16) |  |
| Relu | - |  |
| Linear Layer | (16,512) |  |
| Relu | - |  |
WSS | Linear Layer | (512,32) | 8 | $2.30\times 10^{-4}$
| Relu | - |  |
| Linear Layer | (32,16) |  |
| Relu | - |  |
| Linear Layer | (16,256) |  |
| Relu | - |  |
| Linear Layer | (256,10) |  |
Table 2: The optimized MLP structure and hyperparameters based on Bayesian
tuning

## References

  * [1] D. Mozaffarian, E. J. Benjamin, A. S. Go, D. K. Arnett, M. J. Blaha, M. Cushman, S. De Ferranti, J.-P. Després, H. J. Fullerton, V. J. Howard, et al., Executive summary: heart disease and stroke statistics—2015 update: a report from the american heart association, Circulation 131 (4) (2015) 434–441.
  * [2] K. S. Nayak, J.-F. Nielsen, M. A. Bernstein, M. Markl, P. D. Gatehouse, R. M. Botnar, D. Saloner, C. Lorenz, H. Wen, B. S. Hu, et al., Cardiovascular magnetic resonance phase contrast imaging, Journal of Cardiovascular Magnetic Resonance 17 (1) (2015) 1–26.
  * [3] C. A. Taylor, C. Figueroa, Patient-specific modeling of cardiovascular mechanics, Annual review of biomedical engineering 11 (2009) 109–134.
  * [4] D. A. Steinman, C. A. Taylor, Flow imaging and computing: large artery hemodynamics, Annals of biomedical engineering 33 (12) (2005) 1704–1709.
  * [5] J. Xiang, L. Antiga, N. Varble, K. V. Snyder, E. I. Levy, A. H. Siddiqui, H. Meng, Aview: an image-based clinical computational tool for intracranial aneurysm flow visualization and clinical management, Annals of biomedical engineering 44 (4) (2016) 1085–1096.
  * [6] A. Quarteroni, A. Manzoni, C. Vergara, The cardiovascular system: mathematical modelling, numerical algorithms and clinical applications, Acta Numerica 26 (2017) 365–590.
  * [7] C. A. Taylor, T. A. Fonte, J. K. Min, Computational fluid dynamics applied to cardiac computed tomography for noninvasive quantification of fractional flow reserve: scientific basis, Journal of the American College of Cardiology 61 (22) (2013) 2233–2241.
  * [8] K. M. Chinnaiyan, T. Akasaka, T. Amano, J. J. Bax, P. Blanke, B. De Bruyne, T. Kawasaki, J. Leipsic, H. Matsuo, Y. Morino, et al., Rationale, design and goals of the heartflow assessing diagnostic value of non-invasive ffrct in coronary care (advance) registry, Journal of Cardiovascular Computed Tomography 11 (1) (2017) 62–67.
  * [9] T. W. Hokken, J. M. Ribeiro, P. P. De Jaegere, N. M. Van Mieghem, Precision medicine in interventional cardiology, Interventional Cardiology Review 15 (2020).
  * [10] S. Sankaran, A. L. Marsden, A stochastic collocation method for uncertainty quantification and propagation in cardiovascular simulations, Journal of biomechanical engineering 133 (3) (2011).
  * [11] S. Sankaran, A. L. Marsden, The impact of uncertainty on shape optimization of idealized bypass graft models in unsteady flow, Physics of fluids 22 (12) (2010) 121902.
  * [12] J. S. Tran, D. E. Schiavazzi, A. B. Ramachandra, A. M. Kahn, A. L. Marsden, Automated tuning for parameter identification and uncertainty quantification in multi-scale coronary simulations, Computers & fluids 142 (2017) 128–138.
  * [13] D. Schiavazzi, A. Doostan, G. Iaccarino, A. Marsden, A generalized multi-resolution expansion for uncertainty propagation with application to cardiovascular modeling, Computer methods in applied mechanics and engineering 314 (2017) 196–221.
  * [14] H. Gao, J.-X. Wang, A bi-fidelity ensemble kalman method for PDE-constrained inverse problems in computational mechanics, Computational Mechanics 67 (4) (2021) 1115–1131.
  * [15] H. Gao, X. Zhu, J.-X. Wang, A bi-fidelity surrogate modeling approach for uncertainty propagation in three-dimensional hemodynamic simulations, Computer Methods in Applied Mechanics and Engineering 366 (2020) 113047.
  * [16] M. Mirramezani, S. L. Diamond, H. I. Litt, S. C. Shadden, Reduced order models for transstenotic pressure drop in the coronary arteries, Journal of biomechanical engineering 141 (3) (2019) 031005.
  * [17] M. Mirramezani, S. C. Shadden, A distributed lumped parameter model of blood flow, Annals of Biomedical Engineering 48 (12) (2020) 2870–2886.
  * [18] M. R. Pfaller, J. Pham, A. Verma, N. M. Wilson, D. W. Parker, W. Yang, A. L. Marsden, Automated generation of 0d and 1d reduced-order models of patient-specific blood flow, arXiv preprint arXiv:2111.04878 (2021).
  * [19] A. Manzoni, A. Quarteroni, G. Rozza, Model reduction techniques for fast blood flow simulation in parametrized geometries, International journal for numerical methods in biomedical engineering 28 (6-7) (2012) 604–625.
  * [20] T. Lassila, A. Manzoni, A. Quarteroni, G. Rozza, A reduced computational and geometrical framework for inverse problems in hemodynamics, International journal for numerical methods in biomedical engineering 29 (7) (2013) 741–776.
  * [21] F. Ballarin, G. Rozza, Pod–galerkin monolithic reduced order models for parametrized fluid-structure interaction problems, International Journal for Numerical Methods in Fluids 82 (12) (2016) 1010–1034.
  * [22] F. Ballarin, E. Faggiano, S. Ippolito, A. Manzoni, A. Quarteroni, G. Rozza, R. Scrofani, Fast simulations of patient-specific haemodynamics of coronary artery bypass grafts based on a pod–galerkin method and a vascular shape parametrization, Journal of Computational Physics 315 (2016) 609–628.
  * [23] T. Lassila, A. Manzoni, A. Quarteroni, G. Rozza, Model order reduction in fluid dynamics: challenges and perspectives, Reduced Order Methods for modeling and computational reduction (2014) 235–273.
  * [24] S. Chaturantabut, D. C. Sorensen, Nonlinear model reduction via discrete empirical interpolation, SIAM Journal on Scientific Computing 32 (5) (2010) 2737–2764.
  * [25] H. Gao, J.-X. Wang, M. J. Zahr, Non-intrusive model reduction of large-scale, nonlinear dynamical systems using deep learning, Physica D: Nonlinear Phenomena 412 (2020) 132614.
  * [26] F. Scarselli, A. C. Tsoi, Universal approximation using feedforward neural networks: A survey of some existing methods, and some new results, Neural networks 11 (1) (1998) 15–37.
  * [27] M. Hutzenthaler, A. Jentzen, T. Kruse, T. Anh Nguyen, P. von Wurstemberger, Overcoming the curse of dimensionality in the numerical approximation of semilinear parabolic partial differential equations, Proceedings of the Royal Society A 476 (2244) (2020) 20190630.
  * [28] L. Liang, M. Liu, C. Martin, J. A. Elefteriades, W. Sun, A machine learning approach to investigate the relationship between shape features and numerically predicted risk of ascending aortic aneurysm, Biomechanics and modeling in mechanobiology 16 (5) (2017) 1519–1533.
  * [29] L. Liang, M. Liu, C. Martin, W. Sun, A deep learning approach to estimate stress distribution: a fast and accurate surrogate of finite-element analysis, Journal of The Royal Society Interface 15 (138) (2018) 20170844.
  * [30] A. Madani, A. Bakhaty, J. Kim, Y. Mubarak, M. R. Mofrad, Bridging finite element and machine learning modeling: stress prediction of arterial walls in atherosclerosis, Journal of biomechanical engineering 141 (8) (2019).
  * [31] A. Balu, S. Nallagonda, F. Xu, A. Krishnamurthy, M.-C. Hsu, S. Sarkar, A deep learning framework for design and analysis of surgical bioprosthetic heart valves, Scientific reports 9 (1) (2019) 1–12.
  * [32] F. Sahli Costabal, Y. Yang, P. Perdikaris, D. E. Hurtado, E. Kuhl, Physics-informed neural networks for cardiac activation mapping, Frontiers in Physics 8 (2020) 42.
  * [33] G. Kissas, Y. Yang, E. Hwuang, W. R. Witschey, J. A. Detre, P. Perdikaris, Machine learning in cardiovascular flows modeling: Predicting arterial blood pressure from non-invasive 4d flow mri data using physics-informed neural networks, Computer Methods in Applied Mechanics and Engineering 358 (2020) 112623\.
  * [34] L. Liang, W. Mao, W. Sun, A feasibility study of deep learning for predicting hemodynamics of human thoracic aorta, Journal of biomechanics 99 (2020) 109544\.
  * [35] G. Li, H. Wang, M. Zhang, S. Tupin, A. Qiao, Y. Liu, M. Ohta, H. Anzai, Prediction of 3d cardiovascular hemodynamics before and after coronary artery bypass surgery via deep learning, Communications biology 4 (1) (2021) 1–12.
  * [36] J. Brüning, F. Hellmeier, P. Yevtushenko, T. Kühne, L. Goubergrits, Uncertainty quantification for non-invasive assessment of pressure drop across a coarctation of the aorta using cfd, Cardiovascular engineering and technology 9 (4) (2018) 582–596.
  * [37] S. Valette, J. Chassery, Approximated centroidal voronoi diagrams for uniform polygonal mesh coarsening, Computer Graphics Forum 23 (3) (2004) 381–389.
  * [38] S. Bouaziz, A. Tagliasacchi, M. Pauly, Sparse iterative closest point, in: Computer graphics forum, Vol. 32, Wiley Online Library, 2013, pp. 113–123.
  * [39] M. F. Beg, M. I. Miller, A. Trouvé, L. Younes, Computing large deformation metric mappings via geodesic flows of diffeomorphisms, International journal of computer vision 61 (2) (2005) 139–157.
  * [40] S. Durrleman, M. Prastawa, N. Charon, J. R. Korenberg, S. Joshi, G. Gerig, A. Trouvé, Morphometry of anatomical shape complexes with dense deformations and sparse parameters, NeuroImage 101 (2014) 35–49.
  * [41] R. Liaw, E. Liang, R. Nishihara, P. Moritz, J. E. Gonzalez, I. Stoica, Tune: A research platform for distributed model selection and training, arXiv preprint arXiv:1807.05118 (2018).
  * [42] L. Li, et al., A system for massively parallel hyperparameter tuning, Proceedings of Machine Learning and Systems 2 (2020) 230–246.
  * [43] R. H. Pletcher, J. C. Tannehill, D. Anderson, Computational fluid mechanics and heat transfer, CRC press, 2012.
  * [44] C. M. Rhie, W.-L. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, AIAA journal 21 (11) (1983) 1525–1532.
  * [45] A. O. M. C. Kennedy, Bayesian calibration of computer models, Journal of the Royal Statistical Society: Series B (Statistical Methodology) (63) (2001) 425–464.
  * [46] J. Biehler, M. W. Gee, W. A. Wall, Towards efficient uncertainty quantification in complex and large-scale biomechanical problems based on a bayesian multi-fidelity scheme, Biomechanics and modeling in mechanobiology 14 (3) (2015) 489–513.
  * [47] J. Biehler, S. Kehl, M. W. Gee, F. Schmies, J. Pelisek, A. Maier, C. Reeps, H.-H. Eckstein, W. A. Wall, Probabilistic noninvasive prediction of wall properties of abdominal aortic aneurysms using bayesian regression, Biomechanics and modeling in mechanobiology 16 (1) (2017) 45–61.
  * [48] C. M. Fleeter, G. Geraci, D. E. Schiavazzi, A. M. Kahn, A. L. Marsden, Multilevel and multifidelity uncertainty quantification for cardiovascular hemodynamics, Computer methods in applied mechanics and engineering 365 (2020) 113030.


# The geodesic flow of the BGPP metric
is Liouville integrable

Andrzej J. Maciejewski
Janusz Gil Institute of Astronomy, University of Zielona Góra,
ul. Licealna 9, 65-417, Zielona Góra, Poland
e-mail<EMAIL_ADDRESS>
Maria Przybylska
Institute of Physics, University of Zielona Góra,
Licealna 9, 65–417, Zielona Góra, Poland
e-mail<EMAIL_ADDRESS>
and
Galliano Valent
LPMP: Institut de Physique Mathématique de Provence
13100 Aix en Provence, France
e-mail<EMAIL_ADDRESS>
###### Abstract

We prove that the geodesics equations corresponding to the BGPP metric are
integrable in the Liouville sense. The $\mathrm{SO}(3,\mathbb{R})$ symmetry of
the model allows to reduce the system from four to two degrees of freedom.
Moreover, solutions of the reduced system and its degenerations can be solved
explicitly or reduced to a certain quadrature. In degenerated cases BGPP
metric coincides with the Eguchi-Hanson metric and for this case the mentioned
quadrature can be calculated explicitly in terms of elliptic integrals.

##  1 Introduction

When looking for gravitational instantons was more fashionable, the articles
describing new hyperkähler metrics flourished. Among these metrics, two
Euclidean triaxial metrics, with the Lie algebra $\mathfrak{su}(2)$ of
infinitesimal isometries, were impressive. The first one, the BGPP metric, was
found by Belinski, Gibbons, Page, and Pope in 1978, see [5]. The second one,
called the Atiyah-Hitchin (A-H) metric, was introduced in [1] and its relation
with monopole scattering was developed in [2].

The BGPP metric is still of importance due to its role of a seed in the
development of geometries in higher dimensions: in [6] Einstein 5D geometries
with a negative cosmological constant, in [7] Einstein-Maxwell metrics for
$D\geq 6$ and in [16] a supersymmetric black-hole in 5D.

Although both these metrics are triaxial, they are somewhat different in
complexity. Indeed, components of the BGPP metric are given by simple
algebraic functions, while in the case of A-H they are given by elliptic
integrals, and it seems that they cannot be reduced to a simpler form.

A second difference is that the first metric is not _complete_ while the
second one, which describes the relative motion of two magnetic monopoles,
happens to be complete.

Both metrics exhibit biaxial limits: Eguchi-Hanson for the BGPP metric and
Taub-NUT for the A-H metric [12].

The BGPP metric lies in the multicentre family as proved in [13], [14],
characterized by at least one tri-holomorphic Killing vector, while A-H metric
does not belong to this family.

Considering only the multicentre metrics, the analysis of the integrability of
their geodesic flow was developed in [8] for Taub-NUT and the two-centre
metrics while in [20] a more general analysis led to seven metrics mostly
related with Bianchi A metrics as shown in [21].

Despite these results and further work in [14] and in [21], the integrability
status of the BGPP metric remained unknown.

The aim of this article is to study the integrability of their geodesic flow.
In Section 2 we review the known integrability results for the Multicentre
metrics in order to have a better insight of the BGPP metric. In Section 3 we
lay down the Hamiltonian formulation of the geodesic flow of this metric and
show the integrability in the Liouville sense of this flow. In Section 4 we
analyse the corresponding reduced five-dimensional system which is a Poisson
system with one Casimir function and on a symplectic leaf of this Poisson
structure this system is Hamiltonian with two degrees of freedom and appears
to be integrable in the Liouville sense. Moreover, in Section 5 we describe
its solutions expressible by the Jacobi elliptic functions
$\operatorname{sn}$, $\operatorname{cn}$ and $\operatorname{dn}$. Section 6 is
devoted to degenerate cases of the BGPP metric which appear to coincide with
metrics I and II from paper of Eguchi and Hanson [10]. The corresponding
geodesics flows are integrable as well as the corresponding reduced systems
are integrable and their solutions are expressible either by elliptic or by
trigonometric functions. Some conclusions are presented in Section 7.

##  2 Multicentre versus BGPP metric

The Multicentre metrics are four dimensional Euclidean metrics. Taking for
coordinates $(\sigma,\,x^{i}),\ i=1,2,3$ one has

$g=\frac{1}{V}(d\sigma+\omega)^{2}+V\,\gamma_{0},$ (1)

where $V$ and $\omega$ are independent of $\sigma$ and
$\gamma_{0}=d\vec{x}\cdot d\vec{x}$, where
$\vec{x}=\left(x^{1},x^{2},x^{3}\right)$, is the flat metric in
${\mathbb{R}}^{3}$. It follows that all these metrics share the Killing vector
${\cal K}=\partial_{\sigma}$. Taking

$e_{0}=\frac{1}{\sqrt{V}}(d\sigma+\omega),\quad\quad\quad\quad
e_{i}=\sqrt{V}\,dx^{i},\quad i=1,2,3$

one can construct a triplet of two-forms given by

$\Omega_{i}=e_{0}\wedge e_{i}+\frac{1}{2}\epsilon_{ijk}\,e_{j}\wedge
e_{k}=\frac{1}{2}(J_{i})_{\mu\nu}\,dx^{\mu}\wedge dx^{\nu},\quad\quad
i=1,2,3,$ (2)

with the quaternionic multiplication table

$J_{i}\ J_{j}=-\delta_{ij}\,Id+\epsilon_{ijk}\,J_{k}.$

The closedness of these two-forms follow from

$d\omega=\mathop{\smash{*}}_{\gamma_{0}}\,dV\quad\quad\Longrightarrow\quad\quad\Delta
V=0.$ (3)

By a theorem of Hitchin [15] this construction produces a family of
hyperkähler hence Ricci-flat metrics characterized by tri-holomorphic complex
structures:

$\mathop{\smash{\cal L}}_{{\cal K}}\,J_{i}=0,\quad\quad i=1,2,3.$ (4)

The quest for integrability relies on the existence of a (symmetric) Killing-
Stäckel (K-S) tensor $S_{\mu\nu}$ such that

$\nabla_{(\mu}\,S_{\nu\rho)}=0,$

which gives a quadratic integral

${\cal S}=S^{\mu\nu}\,P_{\mu}\,P_{\nu},\quad\quad\quad\quad\\{H,{\cal
S}\\}=0.$

The first results establishing the Liouville integrability of the flow for the
Multicentre metrics began with the work of Gibbons and Manton [12] for the
Taub-NUT metric for which

$V=1+\frac{m}{r},$

and was generalized in [8] to the two-centre metric for which

$V=v_{0}+\frac{m_{1}}{|\vec{r}+\vec{c}|}+\frac{m_{2}}{|\vec{r}-\vec{c}|}.$

In the first case the construction of a generalized Runge-Lenz vector allowed
for an elegant geometric derivation of the geodesics while in the second case
an educated guess of a set of coordinates separating the Hamilton-Jacobi
equation led to the Killing-Stackel tensor.

A more systematic analysis in [20] led to the following result: let it be
supposed that $V$ has one Killing vector i.e. $K^{i}\,\partial_{i}V=0$. There
are 8 choices of coordinates separating the Hamilton-Jacobi in $\gamma_{0}$
(among the 11 possible choices) which lift up to separating coordinates for
the Hamilton-Jacobi equation for $g$, producing in each case a quadratic
Killing-Stackel tensor. The Liouville integrability stems from the set of
integrals

$H,\quad\quad P_{\sigma},\quad\quad K=K^{i}\,P_{i},\quad\quad{\cal S},$

which include two linear ones and two quadratic ones. This is no longer the
case for the BGPP metric.

The BGPP metric [5] is given by

$g=f(t)^{2}\,dt^{2}+a(t)^{2}\,\sigma_{1}^{2}+b(t)^{2}\,\sigma_{2}^{2}+c(t)^{2}\,\sigma_{3}^{2},\quad\quad
d\sigma_{i}=\frac{1}{2}\,\epsilon_{ijk}\,\sigma_{j}\wedge\sigma_{k}.$ (5)

where $\sigma_{i}$ are the Maurer-Cartan 1-forms for the group
$G=\mathrm{SU}(2)$. Using Euler angles, we have

$\begin{array}[]{lrr}\sigma_{1}=&-\sin\psi\,d\theta&+\sin\theta\,\cos\psi\,d\phi,\\\
\sigma_{2}=&+\cos\psi\,d\theta&+\sin\theta\,\sin\psi\,d\phi,\\\
\sigma_{3}=&d\psi&+\cos\theta\,d\phi.\end{array}$

For the BGPP metric one has

$f(t)^{2}=\frac{1}{4A(t)B(t)C(t)},\,\,\,\,a(t)^{2}=\frac{B(t)C(t)}{A(t)},\,\,\,\,b(t)^{2}=\frac{C(t)A(t)}{B(t)},\,\,\,\,c(t)^{2}=\frac{A(t)B(t)}{C(t)},$
(6)

with

$A(t)=\sqrt{t-t_{1}},\quad\quad B(t)=\sqrt{t-t_{2}},\quad\quad
C(t)=\sqrt{t-t_{3}},$ (7)

compare with formulas (12) and (13) in [5]. The real three parameters
$t_{1},t_{2},t_{3}$ are non-negative and can be chosen arbitrarily. The
coordinate $t$ is restricted to the half line $t>t_{m}$, where $t_{m}$ is the
biggest parameter $t_{i}$.

There is the algebra $\mathfrak{g}=\mathfrak{su}(2)$ of infinitesimal
isometries

$\begin{array}[]{l}\displaystyle\widehat{\cal
L}_{1}=-\sin\phi\,\partial_{\theta}+\frac{\cos\phi}{\sin\theta}\,\partial_{\psi}-\frac{\cos\phi}{\tan\theta}\,\partial_{\phi},\\\\[11.38109pt]
\displaystyle\widehat{\cal
L}_{2}=-\cos\phi\,\partial_{\theta}-\frac{\sin\phi}{\sin\theta}\,\partial_{\psi}+\frac{\cos\phi}{\tan\theta}\,\partial_{\phi},\\\\[11.38109pt]
\widehat{\cal
L}_{3}=\partial_{\phi},\end{array}\quad\quad\quad\quad[\widehat{\cal
L}_{i},\widehat{\cal L}_{j}]=\epsilon_{ijk}\,\widehat{\cal L}_{k}.$ (8)

Using the co-frames

$e_{0}=f(t)dt,\quad e_{1}=a(t)\sigma_{1},\quad e_{2}=b(t)\sigma_{2},\quad
e_{3}=c(t)\sigma_{3},$

it is easy to check that the self-dual 2-forms

$\Omega_{i}=e_{0}\wedge e_{i}+\frac{1}{2}\,\epsilon_{ijk}\,e_{j}\wedge
e_{k}\quad\quad i=1,2,3$ (9)

are closed since we have

$\Omega_{1}=d(A(t)\,\sigma_{1}),\quad\quad\Omega_{2}=d(B(t)\,\sigma_{2}),\quad\quad\Omega_{3}=d(C(t)\,\sigma_{3}).$
(10)

Since the Killing vector $\widehat{\cal L}_{3}=\partial_{\phi}$ is tri-
holomorphic the BGPP metric belongs to the Multicentre family. Indeed, taking
$d\sigma=d\phi$ it was shown in [13] that the potential is given by

$V=\frac{1}{a(t)^{2}\sin^{2}\theta\cos^{2}\psi+b(t)^{2}\sin^{2}\theta\sin^{2}\psi+c(t)^{2}\cos^{2}\theta},$
(11)

and one-form $\omega$ is

$\omega=V\left[\left(b(t)^{2}-a(t)^{2}\right)\sin\theta\sin\psi\cos\psi
d\theta+c(t)^{2}\cos\theta d\psi\right].$ (12)

In these settings coordinates $\vec{x}=\left(x^{1},x^{2},x^{3}\right)$ are

$x^{1}=A\sin\theta\cos\psi,\quad\quad x^{2}=B\sin\theta\sin\psi,\quad\quad
x^{3}=C\cos\theta.$

The problem is that this potential has no symmetry at all. This explains why,
despite the analyses developed in [21] and [14], the integrability status of
its geodesic flow remained unknown for 45 years…

##  3 Hamiltonian formulation of the geodesic flow

The Lagrange function corresponding to the BGPP metric reads

$L=\frac{1}{2}\Big{(}f(t)^{2}\,\dot{t}^{2}+a(t)^{2}\,s_{1}^{2}+b(t)^{2}\,s_{2}^{2}+c(t)^{2}\,s_{3}^{2}\Big{)},$
(13)

where

$\begin{array}[]{l}s_{1}=-\sin\psi\,\dot{\theta}+\sin\theta\,\cos\psi\,\dot{\phi},\\\\[11.38109pt]
s_{2}=+\cos\psi\,\dot{\theta}+\sin\theta\,\sin\psi\,\dot{\phi},\\\\[11.38109pt]
s_{3}=\dot{\psi}+\cos\theta\,\dot{\phi},\end{array}$ (14)

and the dot over a symbol denotes the derivation with respect to parameter
$\lambda$ along geodesic, while the prime denotes the derivation with respect
to coordinate $t$. The canonical momenta conjugated to coordinates
$t,\theta,\phi,\psi$ are the following:

$\begin{array}[]{l}\displaystyle P_{t}=\frac{\partial
L}{\partial\dot{t}}=f(t)^{2}\,\dot{t},\quad\quad P_{\theta}=\frac{\partial
L}{\partial\dot{\theta}}=-a(t)^{2}\,s_{1}\,\sin\psi+b(t)^{2}\,s_{2}\,\cos\psi,\\\\[11.38109pt]
\displaystyle P_{\phi}=\frac{\partial
L}{\partial\dot{\phi}}=a(t)^{2}\,s_{1}\,\sin\theta\,\cos\psi+b(t)^{2}\,s_{2}\,\sin\theta\,\sin\psi+c(t)^{2}\,s_{3}\,\cos\theta,\,\,\,P_{\psi}=\frac{\partial
L}{\partial\dot{\psi}}=c(t)^{2}\,s_{3},\end{array}$ (15)

and Hamilton function takes the form

$\begin{split}H&=\frac{\left(\csc\theta\cos\psi(P_{\phi}-P_{\psi}\cos\theta)-P_{\theta}\sin\psi\right)^{2}}{2a(t)^{2}}+\frac{\left(P_{\theta}\cos\psi+\sin\psi(P_{\phi}\csc\theta-
P_{\psi}\cot\theta)\right)^{2}}{2b(t)^{2}}\\\
&+\frac{P_{\psi}^{2}}{2c(t)^{2}}+\frac{P_{t}^{2}}{2f(t)^{2}}.\end{split}$ (16)

The corresponding Hamilton equations are quite lengthy but one can simplify
them introducing the following variables

$\begin{array}[]{l}\displaystyle
M_{1}=-\sin\psi\,P_{\theta}+\frac{\cos\psi}{\sin\theta}\,P_{\phi}-\frac{\cos\psi}{\tan\theta}\,P_{\psi}=a(t)^{2}\,s_{1},\\\\[11.38109pt]
\displaystyle
M_{2}=+\cos\psi\,P_{\theta}+\frac{\sin\psi}{\sin\theta}\,P_{\phi}-\frac{\sin\psi}{\tan\theta}\,P_{\psi}=b(t)^{2}\,s_{2},\\\\[11.38109pt]
M_{3}=P_{\psi}=c(t)^{2}\,s_{3}.\end{array}$ (17)

The Poisson brackets of these functions are

$\\{M_{i},M_{j}\\}=\epsilon_{ijk}\,M_{k}.$ (18)

In the non-canonical variables $(t,P_{t},M_{1},M_{2},M_{3},\phi,\theta,\psi)$
the Hamiltonian (16) has particularly simple form

$H=\frac{1}{2}\left(\frac{P_{t}^{2}}{f(t)^{2}}+\frac{M_{1}^{2}}{a(t)^{2}}+\frac{M_{2}^{2}}{b(t)^{2}}+\frac{M_{3}^{2}}{c(t)^{2}}\right).$
(19)

However, the symplectic form is no longer constant. In these non-canonical
variables the $8\times 8$ Poisson tensor reads

$J=\begin{bmatrix}J_{2}&0&0\\\ 0&\widehat{M}&-K^{T}\\\ 0&K&0\end{bmatrix},$
(20)

where

$\widehat{M}=\begin{bmatrix}0&M_{3}&-M_{2}\\\ -M_{3}&0&M_{1}\\\
M_{2}&-M_{1}&0\end{bmatrix},\,\,K=\begin{bmatrix}\csc\theta\cos\psi&\csc\theta\sin\psi&0\\\
-\sin\psi&\cos\psi&0\\\
-\cot\theta\cos\psi&-\cot\theta\sin\psi&1\end{bmatrix},\,\,J_{2}=\begin{bmatrix}0&1\\\
-1&0\end{bmatrix}$ (21)

and the explicit form of Hamilton equations corresponding to this structure is
the following

$\begin{split}\dot{t}=&\frac{1}{f(t)^{2}}P_{t},\\\
\dot{P}_{t}=&\frac{f^{\prime}(t)}{f^{3}(t)}P_{t}^{2}+\frac{a^{\prime}(t)}{a^{3}(t)}M_{1}^{2}+\frac{b^{\prime}(t)}{b^{3}(t)}M_{2}^{2}+\frac{c^{\prime}(t)}{c^{3}(t)}M_{3}^{2},\\\
\dot{M}_{1}=&\left(\frac{1}{b(t)^{2}}-\frac{1}{c(t)^{2}}\right)M_{2}\,M_{3}=\frac{(t_{3}-t_{2})}{A(t)B(t)C(t)}M_{2}\,M_{3},\\\
\dot{M}_{2}=&\left(\frac{1}{c(t)^{2}}-\frac{1}{a(t)^{2}}\right)M_{3}\,M_{1}=\frac{(t_{1}-t_{3})}{A(t)B(t)C(t)}M_{3}\,M_{1},\\\
\dot{M}_{3}=&\left(\frac{1}{a(t)^{2}}-\frac{1}{b(t)^{2}}\right)M_{1}\,M_{2}=\frac{(t_{2}-t_{1})}{A(t)B(t)C(t)}M_{1}\,M_{2},\\\
\dot{\phi}=&\csc\theta\left(\frac{M_{1}\cos\psi}{a(t)^{2}}+\frac{M_{2}\sin\psi}{b(t)^{2}}\right),\\\
\dot{\theta}=&\frac{M_{2}\cos\psi}{b(t)^{2}}-\frac{M_{1}\sin\psi}{a(t)^{2}},\\\
\dot{\psi}=&-\frac{M_{1}\cot\theta\cos\psi}{a(t)^{2}}-\frac{M_{2}\cot\theta\sin\psi}{b(t)^{2}}+\frac{M_{3}}{c(t)^{2}}.\end{split}$
(22)

This system is Hamiltonian with four degrees of freedom and for its
integrability four independent commuting first integrals are necessary

###### Theorem 1

The geodesic flow (22) of the BGPP Hamiltonian is integrable in Liouville
sense.

Proof: Apart from the Hamiltonian we have to exhibit 3 integrals in involution
and algebraically independent. The first one follows from the cyclic variable
$\phi$ and it is the corresponding momentum $P_{\phi}$ and additionally exist
two others which have the following forms

$\begin{split}&P_{\phi}=\sin\theta(M_{1}\cos\psi+M_{2}\sin\psi)+M_{3}\cos\theta,\\\
&{\cal C}=M_{1}^{2}+M_{2}^{2}+M_{3}^{2},\\\ &{\cal
I}=t_{1}M_{1}^{2}+t_{2}M_{2}^{2}+t_{3}M_{3}^{2},\end{split}$ (23)

and $H$, $P_{\phi}$, ${\cal C}$ and ${\cal I}$ are independent and pairwise
commuting with the standard Poisson bracket. $\quad\Box$

##  4 Reduced system of the geodesic flow

Let us notice that system (22) contains a closed subsystem

$\begin{array}[]{ll}\displaystyle\quad\quad\dot{t}=\frac{1}{f(t)^{2}}P_{t},\\\\[11.38109pt]
\displaystyle\quad\quad\dot{P}_{t}=\frac{f^{\prime}(t)}{f^{3}(t)}P_{t}^{2}+\frac{a^{\prime}(t)}{a^{3}(t)}M_{1}^{2}+\frac{b^{\prime}(t)}{b^{3}(t)}M_{2}^{2}+\frac{c^{\prime}(t)}{c^{3}(t)}M_{3}^{2},\\\\[11.38109pt]
\displaystyle\quad\quad\dot{M}_{1}=\left(\frac{1}{b(t)^{2}}-\frac{1}{c(t)^{2}}\right)M_{2}\,M_{3}=4f(t)^{2}(t_{3}-t_{2})M_{2}\,M_{3},\\\\[11.38109pt]
\displaystyle\quad\quad\dot{M}_{2}=\left(\frac{1}{c(t)^{2}}-\frac{1}{a(t)^{2}}\right)M_{3}\,M_{1}=4f(t)^{2}(t_{1}-t_{3})M_{3}\,M_{1},\\\\[11.38109pt]
\displaystyle\quad\quad\dot{M}_{3}=\left(\frac{1}{a(t)^{2}}-\frac{1}{b(t)^{2}}\right)M_{1}\,M_{2}=4f(t)^{2}(t_{2}-t_{1})M_{1}\,M_{2}.\end{array}$
(24)

These equations are of the form

$\frac{\mathrm{d}X_{i}}{\mathrm{d}\lambda}=\sum_{j=1}^{5}P_{ij}\frac{\partial
H}{\partial X_{j}},$ (25)

where $X=[t,P_{t},M_{1},M_{2},M_{3}]^{T}$, function $H$ is given by (19), and

$P=[P_{ij}]=\begin{bmatrix}J_{2}&0\\\ 0&\widehat{M}\\\ \end{bmatrix}.$ (26)

So, they are Hamiltonian with Poisson structure given by $P$. This Poisson
tensor defines the Poisson bracket. Namely, for any two functions $F(X)$ and
$G(X)$ we set

$\\{F,G\\}(X)=\frac{\partial F^{T}}{\partial X}P\frac{\partial G}{\partial
X}.$ (27)

Clearly, bracket $\\{\cdot,\cdot\\}$ is antisymmetric and satisfies the
Leibniz rule and the Jacobi identity. However, it is degenerated. The Casimir
function

${\cal C}=M_{1}^{2}+M_{2}^{2}+M_{3}^{2},$ (28)

commute with any function $F(X)$: $\\{{\cal C},F\\}(X)=0$. On a symplectic
leaf ${\cal C}(X)=\alpha>0$ this system becomes Hamiltonian with two degrees
of freedom, thus for its integrability one more functionally independent first
integrals is necessary. This reduction from four to two degrees of freedom is
possible thanks to the $\mathrm{SO}(3,\mathbb{R})$ symmetry of the original
system.

We have the following.

###### Theorem 2

The reduced system corresponding to the geodesic flow of the BGPP Hamiltonian
is integrable in Liouville sense on its symplectic leaf.

Proof: Apart from the Casimir function (28) and Hamiltonian (19) Poisson
system (24) possesses additional first integral ${\cal I}$ defined in (23).
Functions ${\cal C}$, $H$ and ${\cal I}$ are functionally independent and
their Poisson brackets (27) vanish. Moreover, on a level ${\cal
C}(X)=\alpha>0$ the differentials $dH(X)$ are $d{\cal I}(X)$ are linearly
independent. Hence, restriction of $H$ and ${\cal I}$ to this level gives a
Liouville integrable system. $\quad\Box$

##  5 Integration of the reduced system

If the system is integrable, then one would like to give on explicit formulae
for its solutions. However, generally it is a highly non-trivial task. In
particular, it is true for the considered system.

At first, let us notice that the last three equations in the reduced system
(24) have the form of the Euler equations of a rigid body whose moments of
inertia are controlled by dynamical variables $(t,P_{t})$. Our idea how to
effectively perform the integration is similar to that applied in integration
of the Kepler problem: on a given fixed level of first integrals we introduce
implicitly new independent variable $\tau$, then we give solutions as function
of $\tau$. Clearly the last step of this approach has to be a formula that
defines implicitly or explicitly $\tau(\lambda)$. In this section we assume
that $0\leq t_{1}<t_{2}<t_{3}$.

On the common level sets of the first integrals ${\cal C}=m^{2}$, ${\cal
I}=n^{2}$ and $H=e$ we have

$\left(\frac{dt}{d\lambda}\right)^{2}=\frac{P_{t}^{2}}{f(t)^{4}}=\frac{2e}{f(t)^{2}}-4(A(t)^{2}\,M_{1}^{2}+B(t)^{2}\,M_{2}^{2}+C(t)^{2}\,M_{3}^{2})=\frac{2e}{f(t)^{2}}-4m^{2}\,t+4n^{2},$

leading to

$\left(\frac{dt}{d\lambda}\right)^{2}={S}(t)\quad\quad\quad\quad{S}(t)=4\left(n^{2}-m^{2}t+2e\sqrt{(t-t_{1})(t-t_{2})(t-t_{3})}\right),$
(29)

that, in principle, enables to express $t$ as a function of $\lambda$.
However, we do not know if it is possible to find it in an analytic way. In
fact, notice that the genus of the algebraic curve $y^{2}={S}(t)$ is
generically equal to 3.

If we introduce the new independent variable $\tau$ such that

$\frac{\mathrm{d}\tau}{\mathrm{d}\lambda}=4f(t)^{2},$ (30)

then the last three equations of system (24) can be written as

$\frac{dM_{1}}{d\tau}=(t_{3}-t_{2})M_{2}\,M_{3},\quad\quad\frac{dM_{2}}{d\tau}=(t_{1}-t_{3})M_{3}\,M_{1},\quad\quad\frac{dM_{3}}{d\tau}=(t_{2}-t_{1})M_{1}\,M_{2},$
(31)

in which we can recognise Euler equations for a rigid body.

Notice that equations (30) and (29) imply

$\frac{\mathrm{d}\tau}{\mathrm{d}t}=\frac{4f(t)^{2}}{\frac{\mathrm{d}t}{\mathrm{d}\lambda}}=\frac{4f(t)^{2}}{\sqrt{{S}(t)}},$
(32)

that gives function $\tau=\tau(t)$ as a quadrature

$\tau(t)=\int_{t_{3}}^{t}\,\frac{du}{\sqrt{(u-t_{1})(u-t_{2})(u-t_{3})\,{S}(u)}}.$
(33)

The main difference between the Euler equations and equations (31) is that
time is replaced by the variable $\tau$, given by the Abelian integral (33).

Following calculations of explicit solutions for an asymmetric top given e.g.
in [18, Chap. II, Sec. 102] or [17, par. 37] on a common level of first
integrals ${\cal C}=m^{2}$, ${\cal I}=n^{2}$ depending on values of $n^{2}$
and $m^{2}$ explicit solutions are following:
in case I when $t_{3}>\tfrac{n^{2}}{m^{2}}>t_{2}>t_{1}$:

$\begin{split}&M_{1}=\mp\sqrt{\frac{m^{2}t_{3}-n^{2}}{t_{3}-t_{1}}}\,{\rm
cn}\,(\sigma,k),\,\,M_{2}=\pm\sqrt{\frac{m^{2}t_{3}-n^{2}}{t_{3}-t_{2}}}\,{\rm
sn}\,(\sigma,k),\,\,M_{3}=\sqrt{\frac{n^{2}-m^{2}t_{1}}{t_{3}-t_{1}}}\,{\rm
dn}\,(\sigma,k),\\\ &\sigma=\sqrt{(t_{3}-t_{2})(n^{2}-m^{2}t_{1})}\,\tau,\quad
k^{2}=\frac{(t_{2}-t_{1})}{(t_{3}-t_{2})}\frac{(m^{2}t_{3}-n^{2})}{(n^{2}-m^{2}t_{1})}<1;\end{split}$

in case II when $t_{2}>\tfrac{n^{2}}{m^{2}}>t_{1}$:

$\begin{split}&M_{1}=\sqrt{\frac{m^{2}t_{3}-n^{2}}{t_{3}-t_{1}}}\,{\rm
dn}\,(\sigma,k),\,\,M_{2}=\pm\sqrt{\frac{n^{2}-m^{2}t_{1}}{t_{2}-t_{1}}}\,{\rm
sn}\,(\sigma,k),\,\,M_{3}=\mp\sqrt{\frac{n^{2}-m^{2}t_{1}}{t_{3}-t_{1}}}\,{\rm
cn}\,(\sigma,k),\\\ &\sigma=\sqrt{(t_{2}-t_{1})(m^{2}t_{3}-n^{2})}\,\tau,\quad
k^{2}=\frac{(t_{3}-t_{2})}{(t_{2}-t_{1})}\frac{(n^{2}-m^{2}t_{1})}{(m^{2}t_{3}-n^{2})}<1;\end{split}$

in case III when equality $\tfrac{n^{2}}{m^{2}}=t_{2}$ holds

$\begin{split}&M_{1}=m\sqrt{\frac{t_{3}-t_{2}}{t_{3}-t_{1}}}\frac{1}{\cosh\sigma},\quad
M_{2}=m\tanh\sigma,\quad
M_{3}=-m\sqrt{\frac{t_{2}-t_{1}}{t_{3}-t_{1}}}\frac{1}{\cosh\sigma},\\\
&\sigma=m\sqrt{(t_{3}-t_{2})(t_{2}-t_{1})}\,\tau.\end{split}$

Let us notice that $\tau$ is expressible by quadrature of $t$ given in (33)
thus solutions are no longer double-periodic functions of $t$ as it was for a
free rigid body.

Here we remark that our solution is nor fully satisfactory. We did not express
analytically $t$ as a function of $\tau$ and because of this we did not
express $\tau$ as a function of $\lambda$.

##  6 Eguchi-Hanson limit

In the previous section we considered generic case when constants
$(t_{1},t_{2},t_{3})$ are pairwise different. If $t_{2}=t_{3}$, or
$t_{1}=t_{2}$, the either $M_{1}$, or $M_{3}$ is a first integral.
Interpreting equations for $M_{i}$ as Euler equations, these cases correspond
to an axially symmetric body. In both cases $M_{i}(\tau)$ are given by
trigonometric functions. In the paper [5] the metrics corresponding these
limiting cases are called the Eguchi-Hanson metric I and metric II,
respectivly. Both of them are discussed in the second part of the paper. See
also [10] where these metric have been introduced. As it was pointed in these
references the first metric is not complete. So, here we will focus on the
second case only, and we will call it the Eguchi-Hanson metric, or the Eguchi-
Hanson limit.

It is convenient to assume here that $t_{3}<t_{2}<t_{1}<t$. Then the Eguchi-
Hanson limit we are going to investigate is given by $t_{1}=t_{2}$. After the
change of variable $t\,\to\,\rho=\sqrt{t-t_{3}}$ the BGPP metric (5) becomes

$g=\frac{\rho\,d\rho^{2}}{\rho^{2}-\gamma^{2}}+\rho(\sigma_{1}^{2}+\sigma_{2}^{2})+\frac{\rho^{2}-\gamma^{2}}{\rho}\sigma_{3}^{2},\quad\quad\gamma^{2}=t_{1}-t_{3}.$
(34)

Hence, we have $\rho>\gamma$ and in this expression we recognize the Eguchi-
Hanson metric derived in [10]

$g_{\mathrm{II}}=\frac{4r^{4}}{r^{4}-a^{4}}dr^{2}+r^{2}(\sigma_{1}^{2}+\sigma_{2})^{2}+\frac{r^{4}-a^{4}}{r^{2}}\sigma_{3}^{2},$

when we substitute $r^{2}=\rho$ and $a^{4}=\gamma^{2}$. Despite the apparent
singularity $\rho=\gamma$ (a “bolt” according to Gibbons and Hawking) let us
mention that this metric is complete and defined on the cotangent bundle of
the two-sphere [11]. Its geodesic flow at the classical and quantum levels was
first discussed in [19]. However, the author, precluding elliptic integrals,
did not give the explicit flow formulae in general. These can be conveniently
derived within our approach.

The Hamiltonian (19) now reads

$H=\frac{1}{2}\left(\frac{\rho^{2}-\gamma^{2}}{\rho}\,P_{\rho}^{2}+\frac{M_{1}^{2}+M_{2}^{2}}{\rho}+\frac{\rho}{\rho^{2}-\gamma^{2}}M_{3}^{2}\right),$
(35)

and equations (24) transforms into

$\begin{array}[]{l}\displaystyle\dot{\rho}=\frac{(\rho^{2}-\gamma^{2})P_{\rho}}{\rho},\\\
\displaystyle\dot{P}_{\rho}=\frac{1}{2}\left(-\frac{P_{\rho}^{2}(\rho^{2}+\gamma^{2})}{\rho^{2}}+\frac{M_{1}^{2}+M_{2}^{2}}{\rho^{2}}+\frac{M_{3}^{2}(\rho^{2}+\gamma^{2})}{(\rho^{2}-\gamma^{2})^{2}}\right),\\\
\displaystyle\dot{M}_{1}=-\frac{\gamma^{2}\,M_{2}M_{3}}{\rho(\rho^{2}-\gamma^{2})},\\\
\displaystyle\dot{M}_{2}=+\frac{\gamma^{2}\,M_{3}M_{1}}{\rho(\rho^{2}-\gamma^{2})},\\\
\displaystyle\dot{M}_{3}=0.\end{array}$ (36)

This system has first integrals

$H,\quad\quad M_{3},\quad\quad M_{1}^{2}+M_{2}^{2},$

i.e. one of them is linear, and two others are quadratic in variables. For the
BGPP metric the reduced system has three first integrals quadratic in
variables.

On a common set of first integrals $M_{3}=m_{3}$,
$M_{1}^{2}+M_{2}^{2}=\mu^{2}$ and $H=e$ we have

$\left(\frac{d\rho}{d\lambda}\right)^{2}=\frac{R(\rho)}{\rho^{2}},\quad\mbox{where}\quad
R(\rho)=2e\rho^{3}-(\mu^{2}+m_{3}^{2})\rho^{2}-2e\gamma^{2}\rho+\gamma^{2}\mu^{2}.$
(37)

Let us notice that the discriminant of the polynomial $R(\rho)$ is always
positive

${\rm
disc}\,R(\rho)=4\gamma^{2}\Big{[}(\mu^{4}-4\gamma^{2}e^{2})^{2}+m_{3}^{2}\Big{(}\gamma^{2}e^{2}(20\mu^{2}+m_{3}^{2})+\mu^{2}(3\mu^{4}+m_{3}^{4}+3\mu^{2}m_{3}^{2})\Big{)}\Big{]},$
(38)

except for one case which will be considered later. Thus we have three real
roots $\rho_{1}<\rho_{2}<\rho_{3}$ and we put
$R(\rho)=2e(\rho-\rho_{1})(\rho-\rho_{2})(\rho-\rho_{3})$. Since
$R(\pm\gamma)<0$ and $R(0)>0$ we conclude to

$-\gamma<\rho_{1}<0<\rho_{2}<\gamma<\rho_{3}.$

Notice that equation (37) defines domains of possible motion. Curve
$y^{2}=R(\rho)/\rho^{2}$ in half plane has two connected components, see Fig.
1.

Figure 1: Phase curves $y^{2}=R(\rho)/\rho^{2}$

For the first component $0<\rho\leq\rho_{1}$ and for the second
$\rho\geq\rho_{3}$. Since $\rho>\gamma$ we have to consider only the case
$\rho\geq\rho_{3}$.

Now we introduce the variable $\tau$ such that

$\frac{d\tau}{d\lambda}=\frac{1}{\rho(\rho^{2}-\gamma^{2})},$ (39)

which gives

$\frac{dM_{1}}{d\tau}=-\omega\,M_{2},\quad\quad\quad\quad\frac{dM_{2}}{d\tau}=\omega\,M_{1},\quad\quad\quad\quad\omega=\gamma^{2}\,m_{3},$
(40)

and is easily integrated to

$M_{1}=\mu\sin(\omega\tau+\varphi_{0}),\quad
M_{2}=-\mu\cos(\omega\tau+\varphi_{0}),$

where $\mu$ and $\varphi_{0}$ are integration constants.

Equations (37) and (39) imply that

$\left(\frac{d\rho}{d\tau}\right)^{2}=(\rho^{2}-\gamma^{2})^{2}R(\rho).$ (41)

Integrating

$d\tau=\frac{d\rho}{(\rho^{2}-\gamma^{2})\sqrt{R(\rho)}}$ (42)

with assumption $\rho\geq\rho_{3}$ and using formulae 3.137 no 8 from [10] we
obtain

$\tau=\frac{1}{\gamma\sqrt{2e}\sqrt{\rho_{3}-\rho_{1}}(\gamma^{2}-\rho_{1}^{2})}\Big{[}2\gamma\,F(\sigma,k)-(\gamma+\rho_{1})\Pi(\sigma,n_{-},k)-(\gamma-\rho_{1})\Pi(\sigma,n_{+},k)\Big{]},$
(43)

where $F(z,k)$ and $\Pi(z,n,k)$ are the incomplete elliptic integrals of the
first and the third kind defined in 8.111 no 2-4 from [10]. In the above
formula

$\sigma=\arcsin\sqrt{\frac{\rho_{3}-\rho_{1}}{\rho-\rho_{1}}},\quad\quad
n_{\pm}=\frac{\rho_{1}\pm\gamma}{\rho_{1}-\rho_{3}},\quad\quad
k=\sqrt{\frac{\rho_{2}-\rho_{1}}{\rho_{3}-\rho_{1}}}.$

Discriminant (38) is positive except for $m_{3}=0$ and $\mu^{2}=2e\gamma$ in
which case $R(\rho)$ has a double root

$R(\rho)=(\rho-\gamma)^{2}(\rho+\gamma).$

For this case we obtain

$\tau=\frac{1}{4\sqrt{2e}\gamma^{2}}\left[\frac{(\gamma-3\rho)}{(\rho-\gamma)\sqrt{\rho+\gamma}}+\frac{3}{2\sqrt{2\gamma}}\ln\left(\frac{\sqrt{\rho+\gamma}+\sqrt{2\gamma}}{\sqrt{\rho+\gamma}-\sqrt{2\gamma}}\right)\right],\quad\quad\rho>\gamma.$
(44)

##  7 Conclusion

In this article we analyze the integrability of the geodesic flow
corresponding to the BGPP metric in generic and non-generic degenerated cases.
In all these cases the corresponding Hamiltonian systems are integrable in the
Liouville sense. The $\mathrm{SO}(3,\mathbb{R})$ symmetry of the model allows
to reduce the system to two degrees of freedom. After introducing variables
$M_{1},M_{2},M_{3}$ by formulae (17) obtained system splits into two
subsystems and the reduced system also appear to be integrable in all cases.
Moreover, the reduced system can be solved explicitly as functions of a new
variable $\tau$ that is related to $t$ by a quadrature. When the BGPP metric
reduces to the Eguchi-Hanson metric it is possible to calculate explicitly the
mentioned quadrature. An important and interesting question remains open: is
the geodesic flow of the A-H metric integrable?

Acknowledgments: We are greatly indebted to Professor S. Shevchishin for the
determination of the genus of the algebraic curve $y^{2}=S(t)$ in Section 5.

This research has been partially founded by The National Science Center of
Poland under Grant No. 2020/39/D/ST1/01632. For the purpose of Open Access,
the authors have applied a CC-BY public copyright license to any Author
Accepted Manuscript (AAM) version arising from this submission.

## References

  * [1] M. F. Atiyah and N. J. Hitchin, Phys. Lett., A 107 (1985) 21-25.
  * [2] M. F. Atiyah and N. J. Hitchin, The Geometry and Dynamics of Magnetic Monopoles, Princeton University Press (2014).
  * [3] M. Audin, Spinning tops, Cambridge Studies in Advanced Mathematics 51, Cambridge University Press, Cambridge (1996).
  * [4] L. Bates and R. Montgomery, Commun. Math. Phys., 118 (1988) 635-640.
  * [5] V. Belinskii, G. W. Gibbons, D. N. Page and C. N. Pope, Phys. Lett., B 76 (1978) 433-435.
  * [6] D. Dold, Class. Quantum Grav., 35(9), (2018) 095012.
  * [7] A. M. Ghezelbash, Class. Quantum Grav., 39 (7), (2022) 075012.
  * [8] G. W. Gibbons and P. J. Ruback, Comm. Math. Phys., 115 (1968) 267-300.
  * [9] I. S Gradshteyn and I. M. Ryzhik, Table of integrals, series, and products, Eighth edition, Elsevier/Academic Press, Amsterdam (2015).
  * [10] T. Eguchi and A. J. Hanson, Phys. Lett., 74B (1978) 249-251.
  * [11] T. Eguchi and A. J. Hanson, Ann. Phys., 120 (1979) 82-106.
  * [12] G. W. Gibbons and N. S. Manton, Nucl. Phys. B274 (1986) 183-224.
  * [13] G. W. Gibbons, D. Olivier, G. Valent and P. J. Ruback, Nucl. Phys. B 296 (1988) 676-696.
  * [14] G. W. Gibbons, Class. Quantum Grav., 20 (2003) 4401-4408.
  * [15] N. Hitchin, NATO Advanced Study Institute, Les Presses de l’Université de Montréal, Montréal (Québec) (1987).
  * [16] D. Katona and J. Lucietti, Commun. Math. Phys., 399 (2023)
  * [17] L. D. Landau and E. M. Lifshitz, Mechanics,Vol. 1 of Course of Theoretical Physics, Third Edition, Butterworth-Heinenann, Oxford (1997).
  * [18] A. P. Markeev, Teoreticheskaya mekhanika, 2nd ed., R&C Dynamics, Institute of Computer Science, Izhevsk (1999), (In Russian).
  * [19] S. Mignemi, J. Math. Phys., 32 (1991) 3047-3054.
  * [20] G. Valent, Comm. Math. Phys., 244 (2004) 571-594.
  * [21] G. Valent and H. Ben Yahia, Class. Quantum Grav., 24 (2007) 255-276.
  * [22] M. P. Wojtkowski, Bull. Amer. Math. Soc., 18 (1988) 179-183.