# System size dependence of energy loss and correlations of heavy mesons at
LHC energies

Jiaxing Zhao, Joerg Aichelin, Pol Bernard Gossiaux, and Klaus Werner
SUBATECH, Nantes University, IMT Atlantique, IN2P3/CNRS,
4 rue Alfred Kastler, 44307 Nantes cedex 3, France

###### Abstract

We study the system size dependence of heavy quark (HQ) observables at a
center of mass energy of $\sqrt{s_{\rm NN}}=5.02$ TeV to explore whether it
can provide further constraints on the physical processes which are involved:
energy loss of HQs in the quark gluon plasma (QGP), hadronization and hadronic
rescattering. We use the EPOS4HQ approach to study p-p and 0-10% central O-O,
Ar-Ar, Kr-Kr and Pb-Pb reactions and investigate in detail the momentum change
of heavy quarks from creation until their detection as part of a hadron as
well as the enhancement of heavy baryon production at low $p_{T}$. We
investigate furthermore the origin and the system size dependence of the
azimuthal correlations between the heavy quark $Q$ and the heavy antiquark
$\bar{Q}$ and how one can bypass the problem that correlations are washed out
due to the combinatorial background. We conclude that a systematic study of
the system size dependence of the momentum loss allows to separate the
momentum loss due to the passage through the QGP from the momentum change due
to hadronization and the $Q\bar{Q}$ correlations allow to gain inside into the
different pQCD processes in which the heavy quarks are created.

###### pacs:

12.38Mh

##  I Introduction

Theory predicts that at high temperature and/or high density a plasma of
quarks and gluons (QGP) is formed in which partons are not confined anymore.
One of the main goals of heavy-ion reactions at ultra-relativistc energies is
to study this new form of deconfined matter. In the last years ample evidence
has been gathered that this state exists. The observation of strangeness
enhancement [1], quark scaling of the elliptic flow [2], jet quenching [3],
and quarkonium suppression [4] points in this direction. It is the task of
future experiments and the challenge for theory to study the properties of the
QGP, not only at the moment when it hadronizes into mesons and baryons but
also before, during the expansion of the QGP when the temperature is
considerably higher than the hadronization temperature of about $T_{c}\approx
160\rm\ MeV$. Because the multiplicity of mesons, which contain only light
valence quarks, is well described by a statistical model with $T=T_{c}$ and
the spectra of light hadrons are perturbed by hadronic rescattering, the
efforts concentrated in the last year on the study of penetrating probes,
which include jets and heavy flavour hadrons. Heavy flavour quarks have the
advantage that they are created only in the initial hard collisions between
projectile and target nucleons because their mass $m_{Q}$ is large compared to
QGP temperatures and that, due to $m_{Q}>\Lambda_{QCD}$, their production can
be described by perturbative QCD (pQCD).

In the last years many theoretical approaches have been advanced to describe
the passage of heavy quarks through the medium, the subsequent hadronization
and hadronic rescattering. They include PHSD [5, 6], in which the heavy quark
physics is embedded in an approach, which described the light hadrons as well.
It is based on the dynamical quasi particle model (DQPM), which respects the
equation of state of strongly interacting matter. The Catania model [7, 8, 9]
is also based on a DQPM approach and describes the expanding medium by a
Boltzmann equation. The LBT model [10] solves as well a Boltzmann equation for
the heavy quarks including elastic and inelastic collisions but the medium is
described by viscous hydrodynamics. Other models like TAMU[11, 12], Duke[13],
and Torino [14, 15] replace the Boltzmann equation by a Langevin equation to
describe the dynamics of heavy quarks and model the expanding QGP by ideal or
viscous hydrodynamics, fitted to different key observables of light hadrons.

Recently we advanced the EPOS4HQ model, which is based on the EPOS4 approach,
[16, 17, 18, 19], to describe the physics of heavy quarks. EPOS4 has been very
successful in describing a multitude of light hadron observables from RHIC to
LHC energies [18]. The agreement between EPOS4 and the experimental
observables indicates that the time evolution of the QGP, including
hadronization and hadronic rescattering, is well reproduced. EPOS4 creates
heavy quarks as well but they do not interact with the QGP. EPOS4HQ embeds the
heavy quark physics in the EPOS4 approach. There the heavy quarks, created in
the EPOS4 interaction points, interact with the QGP partons by elastic and
inelastic collisions, described by pQCD elastic and inelastic cross sections.
Heavy quarks, which have passed a QGP, hadronize either by coalescence or by
fragmentation. The relative weight of both processes depends on the transverse
momentum, $p_{T}$, of the heavy quark. A detailed description of EPOS4HQ,
including the scattering with the plasma constituents and hadronization, can
be found in [20]. There is also shown that EPOS4HQ describes all the presently
available data quite well.

In [21] the EPOS4HQ approach has been applied to pp collisions and it has been
demonstrated that it describes equally well the heavy hadron production in pp
collision, assuming that in pp collisions as well as, in heavy-ion collisions
the same criterion for the formation of a QGP applies, see Ref. [17].

The agreement of the EPOS4HQ approach with p-p as well as with Pb-Pb data
motivates to study systematically the system size dependence of several key
observables. We start out with a short overview over the heavy quark
production processes (section II) and of the QGP in EPOS4 (section III). The
momentum difference between the $c$ and $\bar{c}$-quark at creation and
between the finally observed $D$-mesons as well as the momentum loss of HQ in
a QGP is shown in section IV. This is followed by a study of the system size
dependence of the heavy baryon enhancement (section V) as well as of the
change of the azimuthal correlation between a $D$ and a $\bar{D}$ meson,
created in the same elementary collision, due to the QGP (section VI).
Finally, in section VII, we present our conclusions.

The calculations are based on version EPOS4.0.1.s9.

##  II heavy quark production and evolution in EPOS4

In EPOS4HQ the heavy quarks are produced together with light quarks and
gluons. Energy and momentum is conserved at each vertex. The interaction
vertices are given by EPOS4. Therefore the kinetic variables as well as the
initial heavy quark coordinates are identical in EPOS4 and EPOS4HQ. The
difference is, as said, that in EPOS4HQ we follow the heavy quarks through the
QGP where they may interact with the QGP partons and that heavy quarks, which
pass a QGP, may hadronize not only by fragmentation but also by coalescence.

Figure 1: Different possibilities to create heavy flavor: (a) flavor
excitation (space-like cascade (SLC)), (b) flavor creation (Born process), (c)
gluon splitting (time-like cascade (TLC)).

In EPOS4 $Q\bar{Q}$ pairs can be produced in different processes, which are
shown in Fig. 1. Before the partons from projectile and target interact in a
hard process they develop a space like DGLAP cascade (SLC), displayed in the
left panel of Fig. 1. In the DGLAP process gluons can be emitted, which can
split into a heavy quark-antiquark pair. This is a next to leading order
process. In the hard process, Fig. 1, middle, collisions of gluons or light
quark-antiquark pairs can create a $Q\bar{Q}$ pair as well. This is the
leading order (Born) process. The gluons, created in the hard process, are
time like (TLC) and may therefore, in a next to leading order process,
disintegrate into a $Q\bar{Q}$ pair. This process is shown in the right panel
of Fig. 1. As we will see later, the correlations between the $Q$ and the
$\bar{Q}$ are rather different in the three processes.

In EPOS4HQ heavy quarks interact with the partons of the QGP, formed by the
light core partons and gluons. We include both, elastic [22] and radiative
[23] collisions. To calculate the energy exchange of the heavy quark with the
medium we first calculate the interaction rate for elastic and inelastic
collisions. If an interaction takes place we select randomly the momentum of
the scattering partner, a light quark or a gluon, from their corresponding
thermal distribution. The thermal distribution is determined by the local
temperature and the local mean parton velocity, at the position where the
heavy quark is localized. The parton momentum is then boosted into the
computational frame, the projectile/target center of mass system in which the
heavy quark momentum is defined. The scattering cross sections of the heavy
quark with gluons and light quarks are calculated by pQCD matrix elements with
a running coupling constant. The kinematic variables of the final state
particles are drawn randomly from the cross section integrated over energy and
momentum conserving delta functions.

When the QGP has expanded to the critical value of the local temperature or
energy density ($T_{\rm FO}=167\ \rm MeV$ or $\epsilon_{\rm FO}=0.57\ \rm
GeV/fm^{3}$), the system hadronizes and the heavy quarks, which are localized
in the QGP, are converted into a heavy flavor (HF) hadron. In the EPOS4HQ
framework, there are two ways in which the heavy quark can hadronize, either
by fragmentation or by coalescence. The momentum distribution of a HF hadron,
which is produced via the coalescence process, can be calculated via,

$\displaystyle{dN\over d^{3}{\bf P}}$ $\displaystyle=$ $\displaystyle\sum
g_{H}\int\prod_{i=1}^{k}{d^{3}p_{i}\over(2\pi)^{3}E_{i}}f_{i}({\bf p}_{i})$
(1) $\displaystyle\times$ $\displaystyle W_{H}({\bf p}_{1},..,{\bf
p}_{i})\,\delta^{(3)}\left({\bf P}-\sum_{i=1}^{k}{\bf p}_{i}\right),$

where $g_{H}$ is the degeneracy factor of color and spin. $k$ equals 2(3) for
mesons (baryons). ${\bf P}$ and ${\bf p}_{i}$ are the momenta of heavy flavor
hadrons and the quark momenta, respectively. The delta function conserves the
momentum. $f_{1}({\bf p}_{1})=(2\pi)^{3}\delta^{(3)}\left({\bf p}_{c}-{\bf
p}_{1}\right)$ is the normalized momentum space distribution of the heavy
quark and $f_{i}({\bf p}_{i})$ for $i>1$ is the momentum space distribution of
the light quarks at the hadronization point. $W_{H}({\bf p}_{1},..,{\bf
p}_{i})$ is the momentum distribution of the quarks bound in a heavy hadron
$H$ in momentum space, corresponding to the Wigner density distribution after
integration over coordinate space. It can be constructed from the hadron wave
function. The hadron wave function can be approximated by a three-dimensional
harmonic oscillator state with the same root mean square radius. For the
ground state of charmed mesons, the Wigner density in the center-of-mass (CM)
frame can be expressed as,

$\displaystyle W(p_{r})$ $\displaystyle=$
$\displaystyle{(2\sqrt{\pi}\sigma)^{3}}e^{-\sigma^{2}{p_{r}}^{2}},$ (2)

where $p_{r}={|E_{2}{\bf p}_{1}-E_{1}{\bf p}_{2}|/(E_{1}+E_{2})}$ is the
relative momentum between the two constituent quarks in their CM frame. It is
normalized, $\int W(p_{r})d^{3}{\bf p}_{r}/(2\pi)^{3}=1$. $E_{1}({\bf p}_{1})$
and $E_{2}({\bf p}_{2})$ are the energies and momenta of the quark or
antiquark in the CM frame, respectively. The width $\sigma$ in the Wigner
density is controlled by the root-mean-radius. If heavy quarks hadronize by
fragmentation we apply the HQET-based fragmentation function [24, 25]. The
probability for coalescence, $P_{coal}$, and of fragmentation, $1-P_{coal}$,
depends on the momentum of the heavy quark. It is displayed in [20].
Fragmentation and coalescence act differently: Whereas the momentum of the
$D$-meson formed by coalescence is usually larger than that of the $c$-quark
due to the momentum of the light quark, which is picked up, the momentum of
the $D$-meson formed by fragmentation is smaller.

Figure 2: Temperature profile of one event in the transverse plane at the
initial time $\tau_{0}$ in p-p, O-O, Ar-Ar, Kr-Kr, and Pb-Pb collisions at
$\sqrt{s_{\rm NN}}=5.02$ TeV and at $z=0$.

##  III Hot QCD medium in EPOS4

EPOS4 and EPOS4HQ use event by event hydrodynamics. This means that in every
simulation the interaction points and the parton productions are different.
The parton energy distribution is the initial condition for the hydrodynamic
evolution of QGP. Typical temperature profiles of a single event in the plane
perpendicular to the beam direction for $z=0$ at the time $\tau_{0}$
($\tau_{0}$=0.4, 1.024, 1.078, 1.167, 1.457 fm/c for p-p, 0-10% O-O, Ar-Ar,
Kr-Kr, and Pb-Pb collisions) are shown in Fig. 2 for central O-O, Ar-Ar, Kr-
Kr, and Pb-Pb events. $x$ is the direction of the impact parameter. We observe
strong fluctuations of the temperature whose scale is given in GeV. Large
areas are above the critical temperature of 0.167 GeV above which a QGP
exists. So almost all heavy quarks, in O-O as well as in Pb-Pb collisions,
pass a QGP. The hydrodynamic evolution of the QGP stops when the critical
energy density of $\epsilon_{\rm FO}$ = 0.57 GeV/fm3 is reached. Then the
system hadronizes, as described in [17].

##  IV heavy quark momentum loss

Figure 3: $p_{T}$ spectra of charm quarks at creation, before hadronization as
well as of $D^{0}$ after hadronization and UrQMD in p-p, O-O, Ar-Ar, Kr-Kr,
and Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV and central rapidity
$|y|<0.5$. Here only charm quarks are included, which are finally part of
$D^{0}$ mesons.

Figure 4: Ratio of the momentum spectrum of charm quarks before hadronization
and that of charm quarks at creation(top left) as well as that of final
$D^{0}$ meson and charm quarks before hadronization (top right) in p-p, O-O,
Ar-Ar, Kr-Kr, and Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV and central
rapidity $|y|<0.5$. The bottom figure displays the ratio of the spectra of
$D$-mesons after and before hadronic rescattering.

On the way to the detector heavy quarks, which are later converted into heavy
mesons, change their momentum by three processes: a) by interaction with the
QGP partons, b) by hadronization, when the heavy quark becomes a heavy hadron
and c) by hadronic rescattering.

In Fig. 3 we display the transverse momentum distribution of heavy quarks and
heavy mesons created at midrapidity ($|y|<0.5$) in minimum bias (MB) for p-p
and central collisions for different symmetric heavy-ion collisions. We select
here only those heavy quarks which hadronize to $D^{0}$-mesons. Therefore the
norm of the spectra for $c$-quarks and $D$-mesons is identical. The black line
shows the distribution of the charm quarks when they are created in the pQCD
processes, the blue line that at the moment before they hadronize into
$D^{0}$-mesons. The red and green lines show the transverse momentum spectrum
of the $D$-mesons immediately after hadronization and after hadronic
rescattering, described by UrQMD, respectively.

We see, first of all, that the red and green lines are almost identical, so
the hadronic rescattering does not modify the spectra considerably. We observe
furthermore that the difference between the black and the blue lines increases
with the size of the colliding system, as we expect, because the spatial
extension of the QGP increases with the system size and we expect a larger
momentum loss if system gets larger. It is remarkable that already in p-p
collisions the black and blue line differ, indicating the formation of a QGP,
as we have pointed out recently [21]. Hadronization changes the spectra in a
non negligible way. For $c$-quarks with a $p_{T}>5$ GeV hadronization by
fragmentation is the dominant process and so we expect that the $p_{T}$ of the
$D$-mesons is smaller that that of the c-quark before hadronization. This
shift is seen in the spectra. At low $p_{T}$ coalescence dominates and
therefore we expect a momentum gain in the hadronization process, as also seen
in the spectra.

To elucidate the physics from a different point of view we display in Fig. 4
ratios of the $p_{T}$ spectra for different projectile-target combinations.
Top left we display the ratio of the $c$-quark spectra just before
hadronization and at creation. Top right shows the ratio of the $D$-meson
$p_{T}$ distribution after hadronization and the $p_{T}$ distribution of the
charm quark before hadronization, hence the change of the spectrum due to
hadronization. The bottom figure shows the ratio of the $D$-meson spectra
after and before UrQMD and consequently the influence of the hadronic
rescattering.

For $p_{T}>5$ GeV the ratio before hadronization/at creation depends strongly
on the system size. For Pb-Pb the spectrum before hadronization is suppressed
by more than a factor of 5 but even for a system as small as O-O the
suppression is close to a factor of 2, so clearly visible. At low $p_{T}$ the
heavy quark rescattering leads to a partial thermalization of the heavy quarks
with the plasma partons. This enhances the hadronization $p_{T}$ spectra
around the thermal equilibrium value of $p_{T}$. The ratio of the $p_{T}$
spectrum of the $D$-meson immediately after hadronization and of the $c$-quark
before hadronization, displayed on the right hand side, shows that the
$D$-meson spectrum is suppressed as compared to that of the $c$-quarks. For
$p_{T}>5$ GeV this suppression of the order of 2 is almost independent of the
size of the colliding system. This confirms that hadronization is a local
process, which happens when the expanding QGP has reached the critical
temperature and which is therefore almost independent of the system size. This
shift is a superposition between the momentum increase due to hadronization by
coalescence and the momentum loss due to hadronization by fragmentation. Fig.
4 bottom, shows the ratio of the $p_{T}$ spectrum after and before hadronic
rescattering, which is close to one for $p_{T}>$ 5 GeV.

Figure 5: System size dependence of the ratio of $p_{T}$ integrated transverse
momentum spectra for $p_{T}>7\;\rm GeV/c$. On the left we display the ratio
before hadronization/at creation, on the right the ratio after UrQMD/at
creation for $\sqrt{s_{\rm NN}}=5.02$ TeV and at central rapidity $|y|<0.5$.

In Fig. 5 we compare integrated $p_{T}$ spectra (for $p_{T}>7$ GeV) as a
function of $A^{1/3}$ to quantify more their system size dependence. On the
left hand side we present the ratio before hadronization/at creation, on the
right hand side before UrQMD/at creation. As shown, the final spectrum equals
the spectrum before UrQMD. We see again from the left figure that larger
systems shift the high momentum particles more to lower $p_{T}$ values than
smaller systems. In Pb-Pb 90% of the heavy quarks with an initial $p_{T}>7$
GeV are shifted to momenta below 7 GeV. In p-p collisions this is only the
case for 20%. For smaller systems the fraction of $c$-quarks shifted to lower
$p_{T}$ is roughly proportional to $A^{1/3}$, and so proportional to the
average length, which the heavy quarks have to travel in the QGP. If we
compare the ratio after UrQMD/at creation on the left hand side we find a
curve which is characteristic for the relative fraction of the momentum loss
due to collisions in the QGP and due to hadronization. This curve can be
normalized to the pp values and compared to the experimental finding.
Therefore, a systematic study of the system size dependence offers the
possibility to separate experimentally the momentum loss of heavy quarks in
the QGP (which increases with the system size) and that due to hadronization
(which is system size independent and dominates at small systems) just by
comparing the integrated $p_{T}$ spectra of $D$ mesons. This would also allow
to calibrate better the transport models which do not agree on the relative
fraction of both momentum losses [26, 20].

The microscopic description of the time evolution of heavy quarks in EPOS4HQ
allows to follow the heavy quarks individually until their hadronization. We
can therefore calculate the momentum transfer, which a heavy quark suffers
while it is moving in the QGP, as well as the momentum difference between the
$c$-quark before hadronization and the $D$-meson after hadronization.

Figure 6: Transverse momentum change in the hot medium (left) and during
hadronization (right) in p-p, O-O, Ar-Ar, Kr-Kr, and Pb-Pb collisions at
$\sqrt{s_{\rm NN}}=5.02$ TeV and central rapidity $|y|<0.5$. $\Delta p_{T}$ is
in GeV.

For $\sqrt{s_{\rm NN}}=$ 200 GeV Au-Au collisions, this momentum loss has
already been calculated using another hydrodynamical model [27]. There a
momentum shift of approx 7 GeV was observed for an initial $c$-quark momentum
of 20 GeV.

In Fig. 6 we present the results for $\Delta p_{T}$ (top row), as well as
those for $\Delta p_{T}/p_{T}$ (bottom row). On the left hand side we see, for
large $p_{T}$, an increase of the momentum loss in the QGP with $p_{T}$.
Because for all systems $\Delta p_{T}/p_{T}$ is almost constant at high
$p_{T}$, the momentum change of a heavy quark in a Pb-Pb collisions with
$p_{T}$= 20 GeV is almost a factor of four larger than for a quark at 5 GeV.
This is a consequence of the energy dependence of the cross sections and the
kinematic of the reaction. Although for all systems the momentum shift as a
function of $p_{T}$ is almost linear, the slope is different for the different
systems due to different lengths, which the heavy quarks travel in the QGP. At
low $p_{T}$, where the heavy quarks are slow and stay for a long time in
contact with the QGP particles, we see that the $c$-quarks are in the process
to come to an equilibrium with the QGP, as already observed in [28]. Therefore
the heavy quarks may gain momentum.

Also the momentum shift due to hadronization, displayed at the right hand
side, depends on $p_{T}$. Here we observe roughly a factor of two between
$p_{T}$= 10 GeV and $p_{T}$= 20 GeV: The fragmentation function scales with
the momentum fraction of the $D$-meson with respect to the momentum of the
$c$-quark. Hence the momentum shift of a $D$-meson produced by a $p_{T}$ = 10
GeV $c$-quark should be roughly half of that produced by a $p_{T}$ = 20 GeV
$c$-quark, what is indeed observed. The momentum shift due to hadronization is
system size independent due to the local character of hadronization. At low
$p_{T}$ the momentum change is positive due to the coalescence process.

To explore the relative importance of the physical processes to which the
finally observed $D$-mesons (here are all charmed hadrons before UrQMD) are
sensitive to, we divide both ratios, the ratio of the momentum shift of the
$c$-quark due to interactions with the QGP and the momentum shift caused by
hadronization. The result is displayed in Fig. 7. Let us start with the high
$p_{T}$ part. Since both, momentum shift in the plasma and momentum shift due
to hadronization, increase with $p_{T}$, the ratio is rather flat. For light
projectile-target combination the momentum shift due to hadronization is
larger than that due to interactions with the QGP. The final $D$-meson $p_{T}$
spectrum is there more sensitive to the fragmentation process (and its correct
description in the medium) then to the energy loss of the $c$-quark in the QGP
and is therefore little informative about the $c$-quark interactions with the
QGP. For larger systems this ratio may decrease by a factor of three and
therefore the final spectra become more sensitive to the $c$-quark
interactions with the QGP.

Figure 7: Ratio of the transverse momentum shift in the QGP and due to
hadronization for p-p, O-O, Ar-Ar, Kr-Kr, and Pb-Pb collisions at
$\sqrt{s_{\rm NN}}=5.02$ TeV and at central rapidity $|y|<0.5$. The two
horizontal lines serve to guide the eye.

##  V Multiplicity of Heavy Baryons

Figure 8: Yield ratio of $\Lambda_{c}/D^{0}$, $\Xi_{c}^{0}/D^{0}$, and
$\Omega_{c}/D^{0}$ in p-p, O-O, Ar-Ar, Kr-Kr, and Pb-Pb collisions at
$\sqrt{s_{\rm NN}}=5.02$ TeV and central rapidity $|y|<0.5$.

The enhanced production of heavy baryons at low $p_{T}$ as compared to the
value expected from $e^{+}e^{-}$ fragmentation functions was one of the
surprising results in p-p collisions. Such an enhancement has also been
observed in Pb-Pb and there it is attributed the the hadronization of heavy
quarks at the end of the expansion phase of the QGP, where heavy quarks may
form heavy baryons by coalescence. A smooth system size dependence of this
enhancement may therefore verify the conjecture that also in very small
systems like p-p a QGP is formed, where heavy baryons can be produced by
coalescence. This conjecture is based on the successful description of
different observables in p-p collision in the EPOS4HQ approach [21].

The system size dependence of the enhancement of heavy baryons at low $p_{T}$
is shown in Fig. 8, which displays from left to right the ratio of the $p_{T}$
spectra of $\Lambda_{c}/D^{0},\Xi_{c}/D^{0}$ and $\Omega_{c}/D^{0}$ for p-p,
O-O, Ar-Ar, Kr-Kr and Pb-Pb collisions. We observe that this enhancement is
rather similar for all systems even if a detailed analysis shows that the
maximum of the enhancement shifts from $p_{T}\approx 3$ GeV for p-p to
$p_{T}\approx$4 GeV for Pb-Pb. Also the enhancement at the maximum increases
by 40% from p-p to Pb-Pb. If this smooth system size dependence, and
explicitly the small difference between p-p and O-O, is verified by
experiments it is hard to imagine that a p-p collision proceeds differently
than a heavy ion collision.

##  VI Correlations

Another source of information, not only on the interaction of $c$-quarks with
the QGP but also on the $c$-quark production mechanism itself, are
correlations among heavy hadrons. They have been measured in p-p collisions at
forward rapidity by the LHCb collaboration [29] and been successfully
described by EPOS4HQ [21, 19].

Figure 9: The correlation between the rapidity $y$ of $c$ and $\bar{c}$ at
creation separated for the different production processes in Pb-Pb collisions
at $\sqrt{s_{\rm NN}}=5.02$ TeV.  Figure 10: The distribution of the total
rapidity of $c$ and $\bar{c}$ for the different production processes in 0-10%
Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV. All curves are normalized to
1.

As discussed in section II, in collisions among hadrons a heavy flavour quark-
antiquark ($Q\bar{Q}$) pair may be produced in three different ways, shown in
Fig. 1. The correlation between $Q$ and $\bar{Q}$ is rather different for the
three different pQCD creation processes.

###  VI.1 Correlations in rapidity

The rapidity correlation between $Q$ and $\bar{Q}$ for the three production
processes as well as for all three processes together is displayed in Fig. 9.

We see that the correlations for flavour excitation and gluon splitting are
centered around $y_{1}=y_{2}$. For flavor excitation the correlations are
rather weak and therefore the $Q$ and $\bar{Q}$, created at the same vertex,
may have rather different rapidities. This means that in experiments a large
rapidity interval is necessary to register both heavy quarks. For gluon
splitting the distribution is considerably narrower because there is no hard
interaction which may cause a strong rapidity shift. For the flavour creating
process, however, the distribution has a maximum at a finite rapidity
difference between the heavy quarks. Parton distribution functions increase
strongly with decreasing $x=p_{parton}/p_{proton}$ values. Therefore
collisions with different $x$ values for projectile and target are more
probable than symmetric collisions. This is also true for collisions with a
sufficient $\sqrt{s}$ value to produce a $Q\bar{Q}$ pair.

In Fig. 10 we display the distributions of the total rapidity, $y_{1}+y_{2}$,
of the $Q\bar{Q}$ pairs. They are normalized to one. All peak at $y_{total}=0$
but the distribution caused by the flavor creation is much broader.

###  VI.2 Correlation in transverse momentum

The correlation between the absolute value of the transverse momentum of the
$Q$ and $\bar{Q}$ of those pairs, where both quarks have a rapidity $|y|<1$,
is displayed for central Pb-Pb collisions in Fig. 11, as well separated for
the three different processes. We display the initial distribution (top) and
the distribution before hadronization, after the passage through the QGP
(bottom).

Due to momentum conservation we expect that initially quarks, coming from the
flavour creation process, have very similar transverse momenta and this is
indeed seen in this LO process where the $c\bar{c}$ pair is emitted nearly
back-to-back. For flavour excitation, where one of the heavy quarks gets a
kick due to the hard interaction, we expect a transverse momentum difference
between the heavy quarks, what is also observed.

Figure 11: The correlation between $p_{T}$ of $c$ and $\bar{c}$ at creation
(top) and before hadronization (bottom) of different processes in Pb-Pb
collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV and at central rapidity $|y|<1$.
Figure 12: The distribution of the transverse momentum difference $|\Delta
p_{T}|$ of charm quarks with $p_{T}>3~{}\rm GeV/c$ and from different
processes in 0-10% Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV. All
curves are normalized to 1.On the left side we display the distribution at the
moment of the production of the $Q\bar{Q}$ pair, on the right side when the
heavy quarks hadronize.

Gluon splitting shows a rather weak correlation, which peaks at
$p_{T,1}\approx p_{T,2}$. The dominating process for the creation of a
$Q\bar{Q}$ pair is flavour excitation.

If we compare the top and the bottom figure we see that the correlations are
considerably weakened during the travel of the heavy quarks through the QGP,
due to their elastic and inelastic collisions. In addition, the transverse
momentum is lowered. The $p_{T,1}\approx p_{T,2}$ correlation of the pairs
created by flavor creation has disappeared and for all creation processes most
of the pairs have quite different $p_{T,1}$ and $p_{T,2}$ values. If we
investigate smaller systems or more peripheral reactions we expect
correlations intermediate between the top and the bottom figure.

This figure shows as well that, if we demand that the $p_{T}$ of both heavy
quarks is larger than 4 GeV, we select pairs which are created by flavour
excitation and not by flavour creation.

To quantify the influence of the QGP on the relative momentum between the
heavy quarks $|\Delta p_{T}|$ = $|p_{T1}-p_{T2}|$ we display in Fig. 12 this
quantity for Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV for the
different production mechanisms, on the left hand side at production, on the
right hand side before hadronization. Here both quarks have to have a $p_{T}>$
3 GeV and the distributions are normalized to one. We see that initially
$Q\bar{Q}$ pairs created by flavour creation have a small momentum difference
whereas those produced in flavor excitation have on the average a considerably
larger $\Delta p_{T}$. Those from gluon splitting are in between. After
traversing the QGP all distributions are rather similar and even that for
flavor excitation increases up to $\Delta p_{T}\to 0$ (see left hand side).
Therefore a $\Delta p_{T}$ cut cannot serve as possibility to separate the
three production processes in central Pb-Pb collisions.

###  VI.3 Azimuthal correlations

Figure 13: Correlations of $D\bar{D}$ (left) and from different processes
(middle and right) in p-p, Ar-Ar, Kr-Kr, and Pb-Pb collisions at $\sqrt{s_{\rm
NN}}=5.02$ TeV and central rapidity $|y|<1$. Here the correlations are charm
quarks or $D$ mesons from the same vertex.  Figure 14: Same as Fig. 13 but for
a wider rapidity range, $|y|<4$.

Next we study the azimuthal correlations between the transverse momenta of the
heavy quarks for the different production mechanisms of the $Q\bar{Q}$ pairs.
These are displayed in Fig. 13 for pairs in which both quarks have $|y|<1$. On
the left hand side we display the distribution of the difference of the
azimuthal angle, $\Delta\phi$, between the $Q$ and $\bar{Q}$ momenta,
separated for the three different production processes, for HQ which both have
$p_{T}>$ 3 GeV. In the middle column we display this correlation of the HQ
after hadronization and in the right column we display the sum of all
production processes for three different $p_{T}$ cuts for both HQ. From top to
bottom we display these correlations for p-p and central Ar-Ar, Kr-Kr and Pb-
Pb collisions.

In the top row, p-p collisions, we see that each of the production processes
creates a different correlation, which is preserved during the hadronization
process. $Q\bar{Q}$s created by flavour creation show an enhancement for back
to back production, as expected by the two body kinematics, which is at the
origin of this type of production. Gluon splitting, on the contrary, produces
an enhancement around $\Delta\phi=0$ because the direction of the momentum of
both quarks is close to that of the emitting gluon. The correlation between
the transverse momenta of those pairs, which are produced by flavour
excitation, is weak due to the kick, which one of the heavy quarks obtains in
the hard process. Summing up the different correlations the correlation
function looks smoother. For $p_{T}<1$ GeV there are practically no strong
azimuthal correlations, which show up when we demand $p_{T}>3$ GeV. There we
see a strong enhancement at $\Delta\phi=0$ and $\Delta\phi=\pi$ coming from
flavour creation and gluon splitting. We can therefore conclude that the
experimentally observed forward-backward enhancement of the correlation
function is due to the three different production processes [21].

When we increase the system size, displayed in the lower rows of Fig. 13 for
Ar-Ar, Kr-Kr and Pb-Pb, we observe that the initial correlations are rather
similar. The passage through the QGP changes the correlations in a different
way. The correlation due to gluon splitting, where both HQ have a similar
$p_{T}$, passes the same region of the QGP and remains rather unchanged
whereas that due to flavor creation is strongly weakened. There are two
reasons. One is the collective radial velocity of the QGP, which is
communicated by collisions to the heavy quarks. If heavy quarks are emitted
back to back, the one with a large momentum component anti-parallel to the
radial collective velocity gets decelerated due to collisions and can inverse
the direction of its momentum, what is not the case for the other quark which
has a large momentum component parallel to the collective velocity. The other
reason is the path length difference in the QGP of back to back emitted
particles . Both processes weaken the azimuthal correlations between the heavy
quarks and depend on the size of the collective expansion velocity (which
increases with the system size) as well as on the size of the QGP. As a
consequence, the correlation function for heavy systems at $\Delta\phi=0$ is
enhanced and the peak at $\Delta\phi=\pi$ disappears. A similar evolution of
the correlation function, especially for flavour creation, has been observed
in PHSD calculations for Au-Au collisions at $\sqrt{s_{\rm NN}}$ = 200 GeV,
using initial $c\bar{c}$ correlations provided by a tuned Pythia [30]
calculation, as well as in other transport approaches [31, 28].

Figure 15: Correlations of $D\bar{D}$ in p-p, O-O Ar-Ar, and Pb-Pb collisions
at $\sqrt{s_{\rm NN}}=5.02$ TeV and central rapidity $|y|<1$.

Because it is in discussion to extend the rapidity range of the ALICE detector
at LHC, in which heavy quarks can be detected, we display in Fig. 14 the
results if both HQ can be observed in a rapidity interval of $|y|<4$. We see
that the correlations are only slightly modified. The flavor excitation
process remains rather flat and the enhancement around $\Delta\phi\approx 0$
is reduced for the gluon splitting process. In the right column we see that
for this larger rapidity interval the forward-backward enhancement, seen in
p-p, is recovered for Pb-Pb collisions if the $p_{T}$ of both HQ is larger
than 6 GeV.

$\Delta\phi>5\pi/6$ | Flavor creation | Flavor excitation | Gluon splitting
---|---|---|---
$p_{T,D},p_{T,{\bar{D}}}$ all | 13.6% | 81.3% | 5.1%
$p_{T,D}>4,2<p_{T,{\bar{D}}}<4$ GeV | 38.3% | 54.6% | 7.1%
$p_{T,D},p_{T,{\bar{D}}}>6$ GeV | 64.3% | 30.6% | 5.1%
$\Delta\phi<\pi/6$ | Flavor creation | Flavor excitation | Gluon splitting
$p_{T,D},p_{T,{\bar{D}}}$ all | 1.4% | 70.8% | 27.8%
$p_{T,D}>4,2<p_{T,{\bar{D}}}<4$ GeV | 1.9% | 35.9% | 62.2%
$p_{T,D},p_{T,{\bar{D}}}>6$ GeV | 1.9% | 23.7% | 74.4%
Table 1: The fractions of $D\bar{D}$ are from different processes when did the
selection, $\Delta\phi>5\pi/6$ and $\Delta\phi<\pi/6$, in p-p collisions with
$\sqrt{s_{\rm NN}}=5.02$ TeV.

Figure 16: Correlations of $D\bar{D}$ (top), $DD$ (middle), $D\bar{D}-DD$
(bottom) in p-p, O-O, Ar-Ar, Kr-Kr, and Pb-Pb collisions at $\sqrt{s_{\rm
NN}}=5.02$ TeV and central rapidity $|y|<4$. Here $D\bar{D}-DD$ correlation is
counted event-by-event. Figure 17: Correlations of $D\bar{D}$ in Pb-Pb
collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV and central rapidity $|y|<4$. From
left to right are with different centrality bins, [0-10%], [10-30%], [30-50%],
and [60-80%].

Up to now we have investigated the ideal situation that one knows at which
vertex the $Q\bar{Q}$ is created. In experiments the vertex cannot be
identified and therefore the experimentally observed correlation function is
that created by all possible $D\bar{D}$ pairs. Their number increases with the
system size and the $D\bar{D}$ pairs from different vertices are uncorrelated.
The correlation function, obtained by including pairs from different vertices,
means the combinatorial background, is shown in Fig. 15 for different system
sizes. We see that pairs coming from different vertices, whose number
increases strongly as compared to those coming from the same vertex, wash out
rather fast the correlation function. For O-O we see for $p_{T}>$ 3 GeV for
both heavy quarks still a structure of the correlation function but very
reduced as compared to p-p. For Ar-Ar the structure of the correlation
function has disappeared. This is also the case if one widens the rapidity
interval as can be seen from the top row of Fig. 16, where we display the same
correlation function. The combinatorial background destroys also here the
structure of the correlation function of pairs coming from the same vertex.

The middle row show the correlation function of two $D$-mesons. In heavy-ion
collisions the two $D$-mesons come in their large majority from different
pomerons and therefore from different production vertices and should not show
a structure in the correlation function. However, if both $D$ mesons contain
heavy quarks from the same pomeron they are correlatedin EPOS4 [19]. This is
the case for p-p collisions.

To quantify these correlations we display in Table 1 the fraction of
$D\bar{D}$ mesons which come from the different processes in the intervals
$\Delta\phi>5\pi/6$ and $\Delta\phi<\pi/6$, for p-p collisions with
$\sqrt{s_{\rm NN}}=5.02$ TeV. We display this percentage for different cuts in
the transverse momentum of the $D$ and $\bar{D}$. We see that without any cut
flavour excitation is dominant in both intervals but cuts in the transverse
momentum allows to select the pQCD process. For $\Delta\phi>5\pi/6$ this is
the flavour creation and for $\Delta\phi<\pi/6$ gluon splitting. So with with
a combination of a high $p_{T}$ cut and a cut in the opening angle one can
exclusively study these pQCD processes individually.

Can this combinatorial background be removed? for experiments event mixing is
a possibility, which is however too costly for transport approaches. Having
the complete information about all $Q$ and $\bar{Q}$ quarks available a
simpler method is to subtract from the $D\bar{D}$ correlation the $DD$
correlation, for which the number of possible pairs is almost identical. In
the bottom row of Fig. 16 we display the difference between the $D\bar{D}$ and
$DD$ correlation and indeed we recover a structure, which resembles the
structure of the original $D\bar{D}$ correlation function of $D\bar{D}$ pairs,
which are coming from the same vertex, shown in Fig. 14. It remains to be seen
how this recovery potential depends on the acceptance of the detectors.

In heavy ion collisions the structure of the correlation function depends only
slightly on the centrality of the collision. This is shown in Fig. 17, where
we display for three centrality intervals [0-10%], [30-50%] and [60-80%] the
$D\bar{D}$ correlation function for three different $p_{T}$ cuts. In the
[30-50%] centrality bin a slight enhancement at $\Delta\phi\approx\pi$
appears, which becomes even more visible in peripheral collisions, [60-80%]
centrality, where also a small enhancement at $\Delta\phi\approx 0$ appears.
The structure remains, however less pronounced than in the combinatorial
background corrected central collisions (see Fig. 16.)

##  VII Conclusion

We presented a study about the system size dependence of several heavy hadron
observables, which have been identified as being sensitive to the interaction
between heavy quarks and the QGP, created in heavy ion reactions. These
observables have the potential to elucidate further the properties of the
expanding QGP, the key objective of the experimental study of
ultrarelativistic heavy ion collisions. For this investigation we used the
recently advanced EPOS4HQ approach, which has been able to describe the heavy
meson single particle spectra in p-p as well as in heavy ion collisions and
reproduced in addition quite nicely the presently available data on heavy
meson correlations.

We observe that indeed the interaction of the heavy quarks with the QGP
hadrons changes the heavy quark momentum distribution or, the other way round,
that the heavy quark spectra carry information about the constituents of the
QGP and the interaction strength. The momentum change depends strongly on the
system size. The other source of the momentum difference between the initial
$c$-quark and the finally observed hadron is hadronization, which is a process
which depends only mildly on the system size. A systematic study of the
momentum change (or the change of the form of the $p_{T}$ spectrum) may
separate the two sources.

The heavy baryon enhancement at low $p_{T}$ is a smooth function of the system
size and continuous from p-p to Pb-Pb. If verified experimentally, this
strengthen the conjecture that this enhancement is a consequence of the
production of a QGP, even in a system as small as p-p.

Azimuthal correlations between the $Q$ and the $\bar{Q}$ at creation are only
weakly modified by shadowing effect. The form of the correlation function is
determined by the three different creation mechanisms of $Q\bar{Q}$ pairs.
Depending on the creation mechanism the momentum difference between the pair
particles is rather differently modified when the heavy quarks pass the QGP.
Also hadronization modifies the different $Q\bar{Q}$ correlation functions. If
one sums over all different production mechanisms the correlation function is
much smoother than that for the individual processes.

In p-p, a cut in $p_{T}$ and in the opening angle allows to separate
experimentally gluon splitting and flavour creation from the other $Q\bar{Q}$
production mechanisms. This is a promising perspective to study these pQCD
processes separately.

The heavy hadrons coming from the same vertex are still correlated, even in
Pb-Pb collisions. Experimentally the vertex, at which the $Q\bar{Q}$ pair is
produced, cannot be identified. With increasing system size this leads to a
the strong increase of the multiplicity of heavy quark pairs from different
vertices (as compared to those from the same vertex). This combinatorial
background makes the correlation function featureless and insensitive to the
interaction of the heavy quark with QGP partons.

The correlation function can be recovered if one subtracts from the
correlation function for $Q\bar{Q}$ pairs that for $QQ$ pairs. If this is
experimentally possible the system size dependence of the correlation function
is another observable which sheds light on the interaction of heavy quarks
with the QGP.

Acknowledgement: This work is supported by the European Union’s Horizon 2020
research and innovation program under grant agreement No 824093 (STRONG-2020).

## References

  * Agakishiev  _et al._ [2012] G. Agakishiev _et al._ (STAR), Phys. Rev. Lett. 108, 072301 (2012), arXiv:1107.2955 [nucl-ex] .
  * Adare  _et al._ [2007] A. Adare _et al._ (PHENIX), Phys. Rev. Lett. 98, 162301 (2007), arXiv:nucl-ex/0608033 .
  * Khachatryan  _et al._ [2017] V. Khachatryan _et al._ (CMS), JHEP 04, 039 (2017), arXiv:1611.01664 [nucl-ex] .
  * Matsui and Satz [1986] T. Matsui and H. Satz, Phys. Lett. B 178, 416 (1986).
  * Song  _et al._ [2015] T. Song, H. Berrehrah, D. Cabrera, J. M. Torres-Rincon, L. Tolos, W. Cassing, and E. Bratkovskaya, Phys. Rev. C 92, 014910 (2015), arXiv:1503.03039 [nucl-th] .
  * Song  _et al._ [2016] T. Song, H. Berrehrah, D. Cabrera, W. Cassing, and E. Bratkovskaya, Phys. Rev. C 93, 034906 (2016), arXiv:1512.00891 [nucl-th] .
  * Plumari  _et al._ [2018] S. Plumari, V. Minissale, S. K. Das, G. Coci, and V. Greco, Eur. Phys. J. C 78, 348 (2018), arXiv:1712.00730 [hep-ph] .
  * Scardina  _et al._ [2017] F. Scardina, S. K. Das, V. Minissale, S. Plumari, and V. Greco, Phys. Rev. C 96, 044905 (2017), arXiv:1707.05452 [nucl-th] .
  * Minissale  _et al._ [2021] V. Minissale, S. Plumari, and V. Greco, Phys. Lett. B 821, 136622 (2021), arXiv:2012.12001 [hep-ph] .
  * Cao  _et al._ [2016] S. Cao, T. Luo, G.-Y. Qin, and X.-N. Wang, Phys. Rev. C 94, 014909 (2016), arXiv:1605.06447 [nucl-th] .
  * He and Rapp [2020] M. He and R. Rapp, Phys. Rev. Lett. 124, 042301 (2020), arXiv:1905.09216 [nucl-th] .
  * He and Rapp [2019] M. He and R. Rapp, Phys. Lett. B 795, 117 (2019), arXiv:1902.08889 [nucl-th] .
  * Cao  _et al._ [2015] S. Cao, G.-Y. Qin, and S. A. Bass, Phys. Rev. C 92, 024907 (2015), arXiv:1505.01413 [nucl-th] .
  * Beraudo  _et al._ [2022] A. Beraudo, A. De Pace, M. Monteno, M. Nardi, and F. Prino, Eur. Phys. J. C 82, 607 (2022), arXiv:2202.08732 [hep-ph] .
  * Beraudo  _et al._ [2023] A. Beraudo, A. De Pace, D. Pablos, F. Prino, M. Monteno, and M. Nardi, (2023), arXiv:2306.02152 [hep-ph] .
  * Werner and Guiot [2023] K. Werner and B. Guiot, Phys. Rev. C 108, 034904 (2023), arXiv:2306.02396 [hep-ph] .
  * Werner [2024a] K. Werner, Phys. Rev. C 109, 014910 (2024a), arXiv:2306.10277 [hep-ph] .
  * Werner [2024b] K. Werner, Phys. Rev. C 109, 034918 (2024b), arXiv:2310.09380 [hep-ph] .
  * Werner [2023] K. Werner, Phys. Rev. C 108, 064903 (2023), arXiv:2301.12517 [hep-ph] .
  * Zhao  _et al._ [2023a] J. Zhao _et al._ , (2023a), arXiv:2311.10621 [hep-ph] .
  * Zhao  _et al._ [2023b] J. Zhao, J. Aichelin, P. B. Gossiaux, and K. Werner, (2023b), arXiv:2310.08684 [hep-ph] .
  * Gossiaux and Aichelin [2008] P. B. Gossiaux and J. Aichelin, Phys. Rev. C 78, 014904 (2008), arXiv:0802.2525 [hep-ph] .
  * Aichelin  _et al._ [2014] J. Aichelin, P. B. Gossiaux, and T. Gousset, Phys. Rev. D 89, 074018 (2014), arXiv:1307.5270 [hep-ph] .
  * Braaten  _et al._ [1995] E. Braaten, K.-m. Cheung, S. Fleming, and T. C. Yuan, Phys. Rev. D 51, 4819 (1995), arXiv:hep-ph/9409316 .
  * Cacciari  _et al._ [2005] M. Cacciari, P. Nason, and R. Vogt, Phys. Rev. Lett. 95, 122001 (2005), arXiv:hep-ph/0502203 .
  * Xu  _et al._ [2019] Y. Xu _et al._ , Phys. Rev. C 99, 014902 (2019), arXiv:1809.10734 [nucl-th] .
  * Gossiaux  _et al._ [2009] P. B. Gossiaux, R. Bierkandt, and J. Aichelin, Phys. Rev. C 79, 044906 (2009), arXiv:0901.0946 [hep-ph] .
  * Nahrgang  _et al._ [2014] M. Nahrgang, J. Aichelin, P. B. Gossiaux, and K. Werner, Phys. Rev. C 90, 024907 (2014), arXiv:1305.3823 [hep-ph] .
  * Aaij  _et al._ [2012] R. Aaij _et al._ (LHCb), JHEP 06, 141 (2012), [Addendum: JHEP 03, 108 (2014)], arXiv:1205.0975 [hep-ex] .
  * Song  _et al._ [2017] T. Song, H. Berrehrah, J. M. Torres-Rincon, L. Tolos, D. Cabrera, W. Cassing, and E. Bratkovskaya, Phys. Rev. C 96, 014905 (2017), arXiv:1605.07887 [nucl-th] .
  * Zhu  _et al._ [2007] X. Zhu, M. Bleicher, S. L. Huang, K. Schweda, H. Stoecker, N. Xu, and P. Zhuang, Phys. Lett. B 647, 366 (2007), arXiv:hep-ph/0604178 .


# Direct Prediction of Steady-State Flow Fields in Meshed Domain with Graph
Networks

Lukas Harsch, Stefan Riedelbauch
Institute of Fluid Mechanics and Hydraulic Machinery
University of Stuttgart
Stuttgart, Germany 70569
<EMAIL_ADDRESS>
###### Abstract

We propose a model to directly predict the steady-state flow field for a given
geometry setup. The setup is an Eulerian representation of the fluid flow as a
meshed domain. We introduce a graph network architecture to process the mesh-
space simulation as a graph. The benefit of our model is a strong
understanding of the global physical system, while being able to explore the
local structure. This is essential to perform direct prediction and is thus
superior to other existing methods.

##  1 Introduction

Using deep learning to model physical problems has become an area of great
research interest. Learning the physical equations of a systems can be used as
model to control robots (Hoppe et al., 2019) or to solve numerical problems
(Raissi, 2018; Sirignano & Spiliopoulos, 2018). In case of computational fluid
dynamics (CFD), learned simulations can reduce the computational cost and thus
improve the engineering process of, in our case, hydraulic components. During
the development process of hydraulic machinery, a large amount of designs has
to be simulated and evaluated with respect to their fluid dynamical
properties. To speed up this process, stationary simulations are often
preferred in the early stage of the design process. For this reason, we
propose a method, which enables fast predictions of the fluid dynamics,
focusing on predicting the steady-state flow field.

In recent years, deep learning approaches (Tompson et al., 2017; Xie et al.,
2018) have shown promising results inferring fluid flow simulations. Data-
driven approaches are used to learn the underlying physical behavior of the
fluid dynamics from CFD simulations. A generalized solution of the fluid
dynamics should be learned to infer the flow field for new problem setups.

A major difference of the algorithms is the way the fluid domain is
represented. Modeling complex physical systems in general requires a
discretized simulation domain. Mesh-based representations are a common choice
to solve the underlying partial differential equations (PDEs) of the system,
e.g. the Navier-Stokes equations for CFD. The mesh can be structured as a
Cartesian grid to store the flow field. So image processing algorithms can be
applied (Guo et al., 2016; Thuerey et al., 2018), which has its limitations
for complex structured geometries. In common numerical solvers the flow field
is represented as a mesh of arbitrary polygons, which is widely used in
mechanical engineering. With the success of graph neural networks (Battaglia
et al., 2018), it is possible to use deep learning for modeling such physical
problems. Further, the fluid can be treated as particle-based representation,
which is common used in smoothed particle hydrodynamics (SPH) and position
based fluids (PBF) for simulating free-surface liquids. Particles can be
processed as a point cloud. In Schenck & Fox (2018), Ummenhofer et al. (2020)
and Sanchez-Gonzalez et al. (2020), this representation is used to model the
highly dynamic movement of fluids and the interaction with rigid objects.

In our work, we focus on mesh-based representation on an irregular grid. The
combination of a differentiable solver with graph convolutions (GCN) (Kipf &
Welling, 2017) is used in De Avila Belbute-Peres et al. (2020) for high-
resolution predictions on the meshed domain around airfoils. Our method comes
without an additional solver to be more independent from domain knowledge
which is necessary for the solver. Pfaff et al. (2021) provides MeshGraphNets,
which is especially designed for predicting the dynamical states of the system
over time. In section 4 we show, that our method is better suited for steady-
state predictions.

##  2 Model

$N^{v}\times F^{v}$InputEdgeConv$N^{v}\times 128$EdgeConv$N^{v}\times
128$EdgeConv$N^{v}\times 128$$N^{v}\times 384$ptGlobal
Features$1024$$N^{v}\times 1024$$N^{v}\times 1408$Decoder$N^{v}\times
3$Outputmlp {1024}max poolingrepeat Figure 1: The model architecture to
predict the flow field, given the node features of shape $N^{v}\times F^{v}$.
Top row: EdgeConv layers as local feature extractors. Top row: Pooling, to
form an 1D global descriptor. Concatenation of the local and global descriptor
and transformation by the decoder. $\oplus$: concatenation.

The task is a steady-state prediction, so the model should learn to directly
predict the final steady-state velocity and pressure fields, given only
information about structure of the meshed domain and additional domain
information, i.e. node positions $\textbf{u}_{i}$ and additional quantities
$\textbf{n}_{i}$, respectively. We propose a graph neural network model
related to Dynamic Graph CNN (DGCNN) (Wang et al., 2019). Figure 1 shows the
visual scheme of our model architecture.

###  2.1 Edge Function

The meshed domain can be seen as a bidirectional graph $G=(V,E)$, with nodes
$V$ connected by mesh edges $E$. Each node $i\in V$ contains the mesh-space
coordinate $\textbf{u}_{i}$, as well as quantities $\textbf{n}_{i}$ describing
the domain, e.g. angle of attack, Mach number or the node type to distinguish
between the fluid domain and solid objects like the geometry. In addition,
each node is carrying the fluid dynamical quantities $\textbf{q}_{i}$, which
we want to model. The node feature vector
$\mathrm{\textbf{v}}_{i}=[\mathrm{\textbf{u}}_{i},\mathrm{\textbf{n}}_{i}]$ is
defined as the concatenation of the nodes position $\mathrm{\textbf{u}}_{i}$
and the quantities $\mathrm{\textbf{n}}_{i}$. The total amount of nodes is
defined by $N^{v}$ with a node feature vector size of $F^{v}$. The edge
feature vectors $\mathrm{\textbf{e}}_{ij}\in E$ are defined as
$\mathrm{\textbf{e}}_{ij}=h_{\Theta}(\mathrm{\textbf{v}}_{i},\mathrm{\textbf{v}}_{j})$,
where $h_{\Theta}$ is a nonlinear function with a set of learnable parameters
$\Theta$.

Extraction both local and global features is essential for the task of
directly predicting a global field, since it is necessary to understand the
local properties of the mesh as well as its global structure. Therefore, we
are using an edge function called EdgeConv introduced by DGCNN, which is
defined as
$h_{\Theta}(\mathrm{\textbf{v}}_{i},\mathrm{\textbf{v}}_{j})=\bar{h}_{\Theta}(\mathrm{\textbf{v}}_{i},\mathrm{\textbf{v}}_{j}-\mathrm{\textbf{v}}_{i})$.
It is implemented as

$\mathrm{\textbf{e}}^{\prime}_{ij}=\mathrm{ReLU}(\Theta\cdot(\mathrm{\textbf{v}}_{j}-\mathrm{\textbf{v}}_{i})+\phi\cdot\mathrm{\textbf{v}}_{i})\;,\qquad\mathrm{\textbf{v}}^{\prime}_{i}=\underset{j:(i,j)\in
E}{\mathrm{max}}\mathrm{\textbf{e}}^{\prime}_{ij}\;,$ (1)

with $\Theta$ and $\phi$ implemented as a shared MLP and $\mathrm{max}$ as
aggregation operation.

The local neighborhood information are captured by
$(\mathrm{\textbf{v}}_{j}-\mathrm{\textbf{v}}_{i})$, which includes the
relative displacement vector in mesh space
$\mathrm{\textbf{u}}_{ij}=\mathrm{\textbf{u}}_{i}-\mathrm{\textbf{u}}_{j}$ and
the difference of the local quantities $\mathrm{\textbf{n}}_{i}$. The global
shape structures are captured by the patch center $\mathrm{\textbf{v}}_{i}$
including the absolute node position $\mathrm{\textbf{u}}_{i}$ and its
quantities $\mathrm{\textbf{n}}_{i}$. The MLP is embedding the feature vector
$\mathrm{\textbf{v}}_{i}$ to a latent vector of size 128. The EdgeConv
formulation has properties lying between translation-invariance and non-
locality. More details about this model type can be found in Wang et al.
(2019).

###  2.2 Architecture

The DGCNN architecture was originally designed for point clouds, using a
k-nearest neighbor algorithm to construct the graph. Since the simulation
domain is already present as a mesh we rely on the graph structure of the mesh
instead of recomputing the graph. We use the structure of the DGCNN
architecture to increase the ability of extracting local and global properties
as displayed in Figure 1. The model has $L$ layers of EdgeConv. The outputs of
the EdgeConv layers can be seen as local feature descriptors. Shortcut
connections are included to extract multi-scale features. The features from
all EdgeConv layers are concatenated and aggregated by a fully-connected
layer. Then, a global max pooling is used to get the global feature vector of
the mesh. Afterwards, the global feature vector is repeated and concatenated
with the local features. At last, the decoder MLP is transforming the latent
features to the desired output $\mathrm{\textbf{p}}_{i}$. Further information
about the network and hyperparameter settings can be found in Section A.1.

(a) Airfoil

(b) AirfoilRot

(c) Channel

Figure 2: Examples of the predicted velocity field in $x$-direction for all
three data sets, each with the ground truth simulation to the right. The
simulation mesh is displayed in the background. For Channel, the geometry and
walls are highlighted in black.

##  3 Experimental setup

For our experiments, three data sets with fluid dynamical simulations on two-
dimensional meshes are used. The flow field is represented by the simulation
mesh. Each node is holding the fluid dynamical quantities
$\mathrm{\textbf{q}}_{i}=[\mathrm{\textbf{w}}_{i},p_{i}]$, with the velocity
in $x$\- and $y$-direction $\mathrm{\textbf{w}}_{i}$ and the pressure $p_{i}$.
All samples in the data set Airfoil are steady-state, compressible, inviscid
CFD simulations, performed by solving the Euler equations. The same simulation
conditions apply to data set AirfoilRot, which is a modified version of
Airfoil. The CFD simulations of the data set Channel are performed by solving
the steady-state, incompressible, inviscid case of the Navier-Stokes
equations.

The first data set Airfoil is taken from Belbute-Peres et al. (2020). The
simulations are describing the flow around the cross-section of a NACA0012
airfoil for varying angle of attacks $\alpha$ and different Mach numbers $m$.
It includes a single geometry setup, so the mesh is identical for every sample
in the data set. Thus, the simulated quantities $\mathrm{\textbf{q}}_{i}$ vary
only due to the parameter setting of $\alpha$ and $m$. The node feature
vectors $\mathrm{\textbf{v}}_{i}$ include the node positions
$\mathrm{\textbf{u}}_{i}$, the global $m$ and $\alpha$ and the node type,
indicating, whether the node belongs to the airfoil geometry, the fluid domain
or the boundary of the flow field. Since the node positions are the same for
every simulation, the model predictions mainly depend on the parameter setup
of $\alpha$ and $m$. Consequently, it is difficult to show the advantage of
our method to capture local and global structures. This data set is used as a
baseline to compare our method with the results from literature.

A second data AirfoilRot set is constructed, by rotating the whole simulation
domain of Airfoil by the value of $\alpha$. This can be interpreted as
$\alpha=0$, but with a rotated placement of the simulation mesh and its flow
field. Here, $\alpha$ will be discarded from the feature vector, since it is
equal for all samples. Now, the model has to rely more on the geometric
structure for predicting the flow field. All samples in the data set still
have the same geometry profile. Nevertheless, the orientation of the airfoil
results in strong variations of the flow field, so the model has to be
sensitive in understanding the geometric structure.

To increase the the requirement of understanding the geometric structure, we
generate a third data set Channel, which is a two-dimensional channel flow.
Generic objects like circles, squares, triangles and rectangle are randomly
placed inside the channel. The number of objects varies between one and three,
which results in a large variety of fluid dynamics i.e., various flows, back
flows, partial blockage and small gaps. The feature vector for Channel only
consists of the node positions $\mathrm{\textbf{u}}_{i}$ and the node type,
i.e. fluid domain, wall or internal object. All samples in Channel are
simulated with the same boundary conditions, so there is no additional
information encoded in the feature vector. Thus, the focus lies on the
investigation for the understanding of the geometric structure. Examples for
all three data sets are shown in Figure 2. More details about the data sets
are explained in Section A.2.

##  4 Results

Table 1: Test root mean square error (RMSE) for all tasks (RMSE$\times 10^{-2}$). Method | Airfoil | AirfoilRot | Channel
---|---|---|---
GCN | 1.40 | - | -
MeshGraphNets | 0.95 | 7.34 | 13.30
MeshGraphNets (Abs) | 0.73 | 1.28 | 13.04
MeshGraphNets (Abs+Pool) | 1.25 | 0.99 | 6.58
Ours | 1.40 | 0.87 | 5.81

We compare our approach with two baseline methods, a base GCN and
MeshGraphNets (Pfaff et al., 2021). The results for the base GCN are taken
from De Avila Belbute-Peres et al. (2020). They are only available for the
Airfoil data set. Then, we use a re-implementation of MeshGraphNets for
comparison with the data sets studied in this work.

MeshGraphNets was designd for modeling the system dynamics over time. The
models spatial equivariance bias is counteracting the direct prediction of a
final flow field. However, the model can achieve good results for Airfoil, but
fails to predict the flow field for the other data sets. We investigate two
extensions of MeshGraphNets, indicated by the suffix (Abs) and (Abs+Pool).
First, the absolute positions $\mathrm{\textbf{u}}_{i}$ is concatenated to the
node features (Abs). Second, we include the pooling operation for global
feature extraction (Abs+Pool). Both extensions are adopted from the DGCNN
architecture and aim to increase the understanding of the global structure.
The adapted model is described more detailed in Section A.3.

The visual results of our method are displayed in Figure 2, showing the
predicted and the ground truth field of the velocity component in
$x$-direction. There is almost no visually noticeable difference between the
predicted and the ground truth flow field. This shows, that our model is able
to produce high-quality predictions. More visual results of the velocity
fields are shown in Figure 4.

Table 1 shows the prediction error in terms of RMSE for all five models on the
three data sets. In the simplest data set Airfoil, with respect to the
complexity of the geometric structure, all models successfully learn to
predict the flow field. As mentioned in Pfaff et al. (2021) MeshGraphNets
outperform the base GCN, which is also the case with our re-implementation.
Surprisingly, (Abs) outperforms all the other methods. In terms of RMSE, our
method falls behind, but visually the fluid dynamics remain nearly identical
to the ground truth.

On the richer AirfoilRot task, MeshGraphNets was unable to obtain a correct
prediction of the flow field. This results in a high RMSE, which is also
confirmed by visual investigations. In contrast, the adapted MeshGraphNets
capture enough geometric information to successfully predict the flow field.
This experiment shows the benefit of using both the local displacement and the
absolute node position as features. In this case, our model outperforms the
other methods in terms of RMSE.

Beside the velocities and pressure, the flow field of a compressible flow
should also be represented by the density, but Airfoil and AirfoilRot do not
provide values for the density. Nevertheless, including the density as an
additional output quantity will increase the complexity of the problem setup
and thus, increase the prediction error slightly. Further, it has to be
mentioned, that the test set of Airfoil and AirfoilRot is within the Mach
number range of the training set. So the test case is an interpolation task.
Since we are not using any additional physical equations or numerical solvers,
our model fails to solve difficult extrapolation tasks as presented in De
Avila Belbute-Peres et al. (2020).

For the most challenging data set Channel, also (Abs) fails. Only (Abs+Pool)
and our method achieve a sufficient low RMSE, but with an improvement of the
RMSE with our method. The same findings are indicated by visual investigations
of the predicted flow fields. This experiment shows, that our method
generalizes well to new unseen geometric setups of the channel flow. In this
experiment, both models with a global feature extractor are successful. This
shows, that extracting the global features with the pooling operation is
necessary to capture a high level of understanding of the global geometric
structure. Nevertheless, our method is better in terms of model complexity and
processing time, which makes it superior to MeshGraphNets (Abs+Pool).

##  5 conclusion

We propose a mesh-based method to model the fluid dynamics for an accurate and
efficient prediction of steady-state flow fields. The experiments demonstrate
the models strong understanding of the geometric structure. The method is
independent of the structure of the simulation domain, being able to capture
highly irregular meshes. Further, our method works well for different kind of
systems, i.e. Eulerian and Navier-Stokes systems. Beside this, we show, that
the model does not require any a-priori domain information, e.g. inflow
velocity or material parameters. Thus, the model can be used for other
systems, e.g. viscous fluids or other domains which are represented as field
data like rigid body simulations.

## References

  * Ba et al. (2016) Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E. Hinton.  Layer normalization, 2016.
  * Battaglia et al. (2018) Peter Battaglia, Jessica Blake Chandler Hamrick, Victor Bapst, Alvaro Sanchez, Vinicius Zambaldi, Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, Caglar Gulcehre, Francis Song, Andy Ballard, Justin Gilmer, George E. Dahl, Ashish Vaswani, Kelsey Allen, Charles Nash, Victoria Jayne Langston, Chris Dyer, Nicolas Heess, Daan Wierstra, Pushmeet Kohli, Matt Botvinick, Oriol Vinyals, Yujia Li, and Razvan Pascanu.  Relational inductive biases, deep learning, and graph networks.  _arXiv_ , 2018.  URL https://arxiv.org/pdf/1806.01261.pdf.
  * Belbute-Peres et al. (2020) Filipe de Avila Belbute-Peres, Thomas D. Economon, and J. Zico Kolter.  Combining Differentiable PDE Solvers and Graph Neural Networks for Fluid Flow Prediction.  In _International Conference on Machine Learning (ICML)_ , 2020.
  * De Avila Belbute-Peres et al. (2020) Filipe De Avila Belbute-Peres, Thomas Economon, and Zico Kolter.  Combining differentiable PDE solvers and graph neural networks for fluid flow prediction.  In Hal Daumé III and Aarti Singh (eds.), _Proceedings of the 37th International Conference on Machine Learning_ , volume 119 of _Proceedings of Machine Learning Research_ , pp. 2402–2411. PMLR, 13–18 Jul 2020.  URL http://proceedings.mlr.press/v119/de-avila-belbute-peres20a.html.
  * Guo et al. (2016) Xiaoxiao Guo, Wei Li, and Francesco Iorio.  Convolutional Neural Networks for Steady Flow Approximation.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, pp. 481–490, New York, NY, USA, 2016. Association for Computing Machinery.  ISBN 9781450342322.  doi: 10.1145/2939672.2939738.
  * Hoppe et al. (2019) Sabrina Hoppe, Zhongyu Lou, Daniel Hennes, and Marc Toussaint.  Planning Approximate Exploration Trajectories for Model-Free Reinforcement Learning in Contact-Rich Manipulation.  _IEEE Robotics and Automation Letters_ , PP:1–1, 2019.  doi: 10.1109/LRA.2019.2928212.
  * Ioffe & Szegedy (2015) Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In Francis Bach and David Blei (eds.), _Proceedings of the 32nd International Conference on Machine Learning_ , volume 37 of _Proceedings of Machine Learning Research_ , pp. 448–456, Lille, France, 07–09 Jul 2015\. PMLR.  URL http://proceedings.mlr.press/v37/ioffe15.html.
  * Kipf & Welling (2017) Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  In _International Conference on Learning Representations (ICLR)_ , 2017.
  * Pfaff et al. (2021) Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzalez, and Peter Battaglia.  Learning mesh-based simulation with graph networks.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=roNqYL0_XP.
  * Raissi (2018) Maziar Raissi.  Deep Hidden Physics Models: Deep Learning of Nonlinear Partial Differential Equations.  _J. Mach. Learn. Res._ , 19:25:1–25:24, 2018.
  * Sanchez-Gonzalez et al. (2020) Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter W. Battaglia.  Learning to simulate complex physics with graph networks.  In _International Conference on Machine Learning_ , 2020.
  * Schenck & Fox (2018) Connor Schenck and Dieter Fox.  SPNets: Differentiable Fluid Dynamics for Deep Neural Networks.  _CoRR_ , abs/1806.06094, 2018.
  * Sirignano & Spiliopoulos (2018) Justin Sirignano and Konstantinos Spiliopoulos.  DGM: A deep learning algorithm for solving partial differential equations.  _Journal of Computational Physics_ , 375:1339–1364, 2018\.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2018.08.029.
  * Thuerey et al. (2018) Nils Thuerey, Konstantin Weissenow, Harshit Mehrotra, Nischal Mainali, Lukas Prantl, and Xiangyu Hu.  Well, how accurate is it? A Study of Deep Learning Methods for Reynolds-Averaged Navier-Stokes Simulations.  _CoRR_ , abs/1810.08217, 2018.
  * Tompson et al. (2017) Jonathan Tompson, Kristofer Schlachter, Pablo Sprechmann, and Ken Perlin.  Accelerating Eulerian Fluid Simulation with Convolutional Networks.  In _Proceedings of the 34th International Conference on Machine Learning - Volume 70_ , ICML’17, pp. 3424–3433. JMLR.org, 2017.
  * Ummenhofer et al. (2020) Benjamin Ummenhofer, Lukas Prantl, Nils Thuerey, and Vladlen Koltun.  Lagrangian Fluid Simulation with Continuous Convolutions.  In _International Conference on Learning Representations_ , 2020.
  * Wang et al. (2019) Yue Wang, Yongbin Sun, Ziwei Liu, Sanjay E. Sarma, Michael M. Bronstein, and Justin M. Solomon.  Dynamic graph cnn for learning on point clouds.  _ACM Trans. Graph._ , 38(5), October 2019.  ISSN 0730-0301.  doi: 10.1145/3326362.  URL https://doi.org/10.1145/3326362.
  * Xie et al. (2018) You Xie, Erik Franz, Mengyu Chu, and Nils Thürey.  tempoGAN: a temporally coherent, volumetric GAN for super-resolution fluid flow.  _ArXiv_ , abs/1801.09710, 2018.

##  Appendix A Appendix

###  A.1 Additional Model Details

All MLPs in the EdgeConv layers are using three shared fully-connected layers
of size 128. The two hidden layers come with ReLU activation, followed by a
non-activated output layer. The MLPs in the EdgeConv layers are followed by a
BatchNorm (Ioffe & Szegedy, 2015) layer. To aggregate the features of all
EdgeConv layers, a fully-connected layer of size 1024 is used. The decoder MLP
has four fully-connected layers of size (1024, 512, 256, 3) to transform the
latent feature vector to the output shape. We are using $L=5$ EdgeConv layers,
which is a good trade-off between model capacity and computational cost. The
whole model consists of 20 fully-connected layers with a total number of
approximately 3.3 million parameters.

We train our model in a supervised fashion using the per-node output features
$\mathrm{\textbf{p}}_{i}$ from the decoder with a $L_{2}$ loss between
$\mathrm{\textbf{p}}_{i}$ and the corresponding ground truth values
$\mathrm{\textbf{p}}_{i}$. A batch size of 2 is chosen for Airfoil and
AirfoilRot. In Channel the number of nodes is smaller compared to the former
data sets, so a batch size of 6 is used. The models are trained on a single
RTX 3090 GPU with the Adam optimizer, using a learning rate of $10^{-4}$.

###  A.2 Dataset

The simulation domain of Airfoil and AirfoilRot is a two-dimensional, mixed
triangular and quadrilateral mesh with 6648 nodes. The mesh is highly
irregular with varying scale of the edge lengths in different regions of the
mesh. The training set is defined by

$\alpha_{\mathrm{train}}=\\{-10,-9,\ldots,9,10\\}\;,\qquad
m_{\mathrm{train}}=\\{0.2,0.3,0.35,0.4,0.5,0.55,0.6,0.7\\}\;.$

The test set is defined similar by

$\alpha_{\mathrm{test}}=\\{-10,-9,\ldots,9,10\\}\;,\qquad
m_{\mathrm{test}}=\\{0.25,0.45,0.65\\}\;.$

The training and test pairs are sampled uniformly from $\alpha\times m$. This
results in 168 samples for the training set and 63 samples for the test set,
which is split to 32 samples for validation and 31 samples for the test. The
samples from Airfoil and AirfoilRot come with a min-max normalization to the
value range of $[-1,1]$ separately for each velocity component and the
pressure.

The simulation domain of Channel is a two-dimensional, unstructured triangular
mesh with an equal scale of the edge lengths. Depending on shape and number of
objects inside the channel, the number of nodes is variable, but is
approximately 1750. This data set has 1000 samples for training and 150
samples each for validation and test. In Channel, we normalize the samples by
the maximum velocity and the maximum pressure.

The node features used as inputs as well as the outputs are summarized in the
following table.

Method | inputs $\mathrm{\textbf{v}}_{i}$ | outputs $\mathrm{\textbf{p}}_{i}$
---|---|---
Airfoil |  $\mathrm{\textbf{u}}_{i}$, $\mathrm{n}_{i}$, $m$, $\alpha$ | $\mathrm{\textbf{w}}_{i},p_{i}$
AirfoilRot |  $\mathrm{\textbf{u}}_{i}$, $\mathrm{n}_{i}$, $m$ | $\mathrm{\textbf{w}}_{i},p_{i}$
Channel |  $\mathrm{\textbf{u}}_{i}$, $\mathrm{n}_{i}$ | $\mathrm{\textbf{w}}_{i},p_{i}$

The inputs are defined by the node position $\mathrm{\textbf{u}}_{i}$, the
node type $\mathrm{n}_{i}$, the angle of attacks $\alpha$ and the Mach numbers
$m$, while the outputs are the velocity in $x$\- and $y$-direction
$\mathrm{\textbf{w}}_{i}$ and the pressure $p_{i}$.

###  A.3 MeshGraphNets Adaption

$N^{v}\times F^{v}$, $N^{e}\times F^{e}$InputEncoder$N^{v}\times 128$,
$N^{e}\times 128$Message Passing$N^{v}\times 128$, $N^{e}\times 128$Message
Passing$N^{v}\times 128$, $N^{e}\times 128$Message Passing$N^{v}\times 128$,
$N^{e}\times 128$$N^{v}\times 640$ptGlobal Feature$1024$$N^{v}\times
1024$$N^{v}\times 1664$Decoder$N^{v}\times 3$Outputmlp {1024}max poolingrepeat
Figure 3: Adapted model architecture using message passing block with local
and global descriptor. $\oplus$: concatenation.

For the adaption of MeshGraphNets we are following their Encoder-Processor-
Decoder architecture, but with a processor adapted by the idea of DGCNN.
Instead of EdgeConv layers, message passing blocks (Battaglia et al., 2018)
are used to update node and edge features. Figure 3 shows the visual scheme of
the adapted architecture.

Encoder The edge feature vectors $\mathrm{\textbf{e}}_{ij}\in E$ are a
concatenation of the relative displacement vector in mesh space
$\mathrm{\textbf{u}}_{ij}=\mathrm{\textbf{u}}_{i}-\mathrm{\textbf{u}}_{j}$ and
its norm $|\mathrm{\textbf{u}}_{ij}|$ with an edge feature vector size of
$F^{e}$ and a total amount of edges $N^{e}$. The node feature vectors
$\mathrm{\textbf{v}}_{i}\in V$ is defined by the quantities
$\mathrm{\textbf{n}}_{i}$ and can be extended by other features with a node
feature vector size of $F^{v}$ a total number of nodes $N^{v}$. In the
encoder, MLPs $\epsilon^{M}$ and $\epsilon^{V}$ are used to embed the feature
vectors $\mathrm{\textbf{e}}_{ij}$ and $\mathrm{\textbf{v}}_{i}$, respectively
to latent vectors of the size 128.

Processor The processor uses $L$ identical message passing blocks to update
the edges e and nodes v into embeddings $\mathrm{\textbf{e}}^{\prime}$ and
$\mathrm{\textbf{v}}^{\prime}$ respectively by

$\mathrm{\textbf{e}}^{\prime}_{ij}\leftarrow
f^{M}(\mathrm{\textbf{e}}_{ij},\mathrm{\textbf{v}}_{i},\mathrm{\textbf{v}}_{j})\;,\qquad\mathrm{\textbf{v}}^{\prime}_{i}\leftarrow
f^{V}(\mathrm{\textbf{v}}_{i},\sum_{j}^{\phantom{x}}\mathrm{\textbf{e}}^{\prime}_{ij})\;,$
(2)

where $f^{M}$ and $f^{V}$ are implemented as MLPs with residual connections.
The MLPs are using the concatenated node and edge feature vectors as input.

Decoder Finally, the decoder MLP is transforming the latent node features from
the last processor block into the output features. More details about this
model type can be found in Sanchez-Gonzalez et al. (2020) and Pfaff et al.
(2021).

Local and global Features EdgeConv captures both the local neighborhood
information $(\mathrm{\textbf{v}}_{j}-\mathrm{\textbf{v}}_{i})$ and the global
shape structures by $\mathrm{\textbf{v}}_{i}$. Therefore, we are adapting the
node features by using a concatenation of the quantities
$\mathrm{\textbf{n}}_{i}$ and the absolute node position
$\mathrm{\textbf{u}}_{i}$. This node features help to capture the global shape
structure, while the edge features capture the local neighborhood information.
This adapted version is called MeshGraphNets (Abs). Further, we adopt the
structure of the DGCNN architecture to increase the ability of extracting
local and global properties. Instead of stacking $L$ message passing block,
all the message passing blocks are concatenated to calculate the global
feature vector. Then, again the local and global features are concatenated and
transformed to the desired output by the decoder MLP. This adapted version is
called MeshGraphNets (Abs+Pool).

All MLPS in the encoder, message passing blocks and decoder are using three
fully-connected layers of size 128. The two hidden layers come with ReLU
activation, followed by a non-activated output layer. The MLPs except the
decoder are followed by a LayerNorm (Ba et al., 2016) layer. Similar to our
method, we are using $L=5$ message passing block.

###  A.4 Additional Results

(a) Airfoil

(b) AirfoilRot

(c) Airfoil

(d) AirfoilRot

(e) Channel

(f) Channel

(g) Channel

Figure 4: Additional examples of the predicted velocity field in $x$-direction
for all three data sets, each with the ground truth simulation to the right.