#  Production locality and spatial diffusion of heavy flavour quarks
at high energy densities

M. Gazdzicki  Jan Kochanowski University, Kielce, Poland  D. Kikoła  Warsaw
University of Technology, Warsaw, Poland  I. Pidhurskyi  Jan Kochanowski
University, Kielce, Poland  European Organization for Nuclear Research, CERN,
Geneva, Switzerland  L. Tinti  Jan Kochanowski University, Kielce, Poland

###### Abstract

Heavy-ion collisions are a unique tool for testing the behaviour of matter
under extreme conditions. The momentum correlations of charm and bottom
hadrons have been considered for testing heavy quarks’ thermalization in the
hot, dense medium produced by the collisions. In this respect, two effects
have been considered: the decrease of the initial back-to-back correlations
and the increase of correlations due to heavy-quark interactions with
collectively flowing medium.

Here, we show that, in the case of a single charm and anti-charm hadron pair
production, the collective flow allows for testing heavy-quark production
locality and spatial diffusion. Using an example of central Pb+Pb collisions
at the CERN SPS energies, we demonstrate that the azimuthal correlations of
charm and anti-charm hadrons are particularly sensitive to their spatial
correlations. We argue that the existing experimental technology and beam
intensities at the CERN SPS should allow for the corresponding measurements
soon. The correlation measurements in collisions with a single heavy-quark
pair produced will provide a unique input constraining the diffusion of charm
quarks and verifying assumptions concerning production locality of a charm and
anti-charm quark pair.

##  I Introduction

Collisions of heavy ions at relativistic energies provide insights into
fascinating features of nuclear matter at high energy densities. This includes
the creation and properties of the Quark-Gluon Plasma (QGP) [1] \- a state of
matter with quark and gluon degrees of freedom expected to exist in the
Universe’s first moments. Moreover, there is a possibility of discovering the
critical point of strongly interacting matter; for example, see Refs. [2, 3]
and references therein. Impressive progress has been made in experimental and
theoretical studies in the last decades. Still, many properties of high-
density matter and particle creation in the medium remain to be uncovered. In
particular, we propose testing the popular assumption of the local production
of particles and antiparticles. This should be experimentally possible by
studying heavy-quark (charm and bottom) correlations at properly selected
collision energies.

Measurements of correlations between the charm meson and its anti-particle
have been proposed to test the equilibration of charm [4, 5] in momentum
space. In a semi-classical picture, the initial back-to-back momentum
correlations between the $c$ and $\bar{c}$ quarks are reduced by the
interactions with the medium and hadronization of the quarks (see, for
instance, Ref. [6] and references therein). Thus, the charm hadron
correlations provide means for quantifying transport properties of the
strongly interacting medium, complementary to measurements of collective
effects (via elliptic flow $\mathrm{v_{2}}$) [7] and modification of momentum
spectra via nuclear modification factor $\mathrm{R_{AA}}$ [8].

In this paper, we present and discuss a different effect related to the
creation mechanism of heavy quarks in a deconfined, partonic system. Using
azimuthal correlations of charmed hadrons at low transverse momentum
$p_{\textrm{T}}$, one can study the spatial distribution of charm quarks.
Specifically, we show that the observed correlations at low $p_{\textrm{T}}$
provide direct insights into whether heavy quarks ($Q$ and $\bar{Q}$) are
produced close to each other in the coordinate space (in other words, if there
are produced locally), or if the production is non-local, i.e. production
points of $Q$ and $\bar{Q}$ are distant. Clearly, the information on the
creation locality is of fundamental importance. Moreover, the creation
mechanism is an input assumption in models of charm and bottom quark
interaction with the quark-gluon plasma, and all the modern experiments at
SPS, RHIC and the LHC conduct programs that aim to quantify the QGP parameters
using heavy quarks. The interpretation of the experimental observations (for
instance, long-range correlations measured with a large difference in rapidity
$\Delta y$) could drastically differ depending on whether one assumes the
heavy quarks are initially created close or distant in the coordinate space.
Thus, experimental verification of the usually assumed local production of a
$Q\bar{Q}$ pair is crucial for understanding particle production in strong
interactions. The study requires selecting heavy-ion collisions having at most
a single $Q\bar{Q}$ pair created. This condition is fulfilled for charm quarks
at the SPS energy and bottom quarks at RHIC at sufficiently low collision
energies.

The idea utilizes the collective flow of charm hadrons measured in heavy-ion
collisions at high energies [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. We
assume that final-state momenta of charm hadrons are given by the
superposition of the flow and a random (thermal) contribution due to
statistical hadronization. Hadronic rescattering and final-state interactions
are neglected, supported by the recent measurements of interaction parameters
of $D$-mesons with hadrons [20, 21]. The flow contribution depends only on the
emission point in the freeze-out hypersurface, whereas the thermal
contribution is a random effect, uncorrelated for different hadrons. This
paper is motivated by the fact that, even for a locally thermalized and
expanding medium, the momenta of charm and anti-charm hadrons originating from
the same $c$ and $\bar{c}$ pair are correlated. Depending on the charm’s
creation points and spatial diffusion properties in the medium, the charm
hadron and the anti-hadron emission points can be either close or distant. In
a locally thermalized and expanding system, the charm hadrons have an average
momentum dependent on the fluid cell’s drift speed (flow). If the emission
points of the hadrons are close, they will have a similar drift, and thus
their momenta will be correlated. In the case of distant emission points the
correlation should vanish.

Here, we stress that the arguments are generally valid, but to directly
measure the wanted correlations, one should have no more than one $c\bar{c}$
pair produced per collision. Otherwise, the measured two-particle correlation
function includes pairs of $c$\- and $\bar{c}$-hadrons coming from different
and likely independent charm production processes. The magnitude of this
unwanted contribution to the momentum-correlation results strongly depends on
the multiplicity distribution of heavy-flavour pairs. This effect is
especially important in the heavy-ion collisions at RHIC and the LHC. On
average, one expects $\simeq 3$ $c\bar{c}$ pairs in the 10% most central Au+Au
collisions at $\mathrm{\sqrt{s_{NN}}=200}$ GeV at RHIC [22, 23], and a few
tens at the LHC (for example, $\simeq 30$ $c\bar{c}$ pairs in the 10% most
central Pb+Pb reactions at $\mathrm{\sqrt{s_{NN}}=5.02}$ TeV [24, 25]). The
multiplicity distribution of heavy flavour is difficult to access
experimentally, and the wanted correlations at very high energies cannot be
extracted in a model-independent way. Thus, the measurements should be
performed at sufficiently low collision energies to minimize the bias due to
unwanted correlations, where the mean multiplicity of $c\bar{c}$ pairs is
below one. For this reason, we consider an example of central Pb+Pb collisions
at the CERN SPS energies. This example can be straightforwardly extended to
bottom hadron production at RHIC or at the LHC.

The heavy-flavour production and azimuthal correlations in heavy-ion
collisions at very high energies were addressed theoretically in the past; for
review, see Ref. [26]. In particular, they were considered as a tool for
uncovering a mechanism behind the jet suppression [27, 28] and the study of
charm energy-loss mechanism [29, 30, 31, 32]. The heavy-quark spatial
diffusion in QCD matter was discussed recently in Refs. [33, 34, 35], see also
references therein. The ATLAS experiment measured the azimuthal-angle
correlations of muon pairs originating from heavy-flavour decays in 5.02 TeV
Pb+Pb collisions [36]. One notes that the measured muon pairs come from jet-
like correlations of high transverse-momentum heavy-flavour hadrons. However,
a possibility to distinguish between local and non-local $Q\bar{Q}$ production
is addressed for the first time in this work.

The paper is organized as follows. First, we briefly discuss theoretical
challenges in predicting ab initio hadron correlations in heavy-ion
collisions, Sec. II. Then, the qualitative idea of testing production locality
is quantified using simple modelling presented in Sec. III. The physics
meaning of different assumptions on charm-hadron spatial correlations in the
emission volume is also discussed in this section. The feasibility of the
corresponding measurements is estimated in Sec. IV, and the results are
summarized in Sec. V.

##  II Theoretical challenges

Quantum Chromodynamics (QCD) is the commonly accepted theory of strong
interactions. However, attempts to derive precise quantitative predictions for
multi-particle production in high-energy collisions from the QCD have been
unsuccessful. Predictions of QCD–inspired models suffer from large
uncertainties. Here, we discuss them in aspects relevant to this work.

The most popular QCD-inspired approaches to predict hadron production in
heavy-ion collisions are based on classical approximations. For instance,
these are the relativistic kinetic theory and hydrodynamic models [37]. In the
context of heavy-ion collisions, the quantum effects are expected to be large
or even dominant with respect to the classical predictions, at least regarding
the flow of energy and momentum [38]. This is because the typical action scale
of the system, a few hundred (at most) of MeV’s of temperature and spatial
changes in a fraction of a femtometer, is smaller than $\hbar c\simeq 200$
MeV$\cdot$fm. The surprising success of hydrodynamical models in describing
nuclear reactions can probably be traced back either to the attractor dynamics
[39] or the generalized off-shell hydrodynamic expansion [40]. In any case,
these arguments hold only for the hydrodynamic variables and not, for
instance, for the two-particle correlations. Operators’ expectation values,
such as the energy density in hydrodynamics, are generally considered. Still,
their fluctuations (e.g., variance and higher-order moments) and related
correlations are more difficult to deal with.

It is important to stress that, despite the success in predicting some
observables, the approximations used in quantitative models are known to be
inadequate for the full description of the data. For instance, hydrodynamics
cannot properly address the high $p_{T}$ part of the spectra; relativistic
kinetic theory assumes molecular chaos and removes two and three-particle
correlations already at the classical level. Most physically proven quantum
effects (diffraction, entanglement, etc.) are neglected and cannot be
addressed by the current models. It is not simple to estimate the size of the
quantum effects, lacking quantitative models including them. Thus, it is
unclear whether the correlations produced by a classical treatment like
relativistic kinematics are enough to describe the experimental results. For
further discussion see Appendix A.

Considering the above, we follow a model-agnostic approach in this work. We do
not consider questionable ab initio calculations of $Q$ and $\bar{Q}$ spatial
and momentum correlations. Instead, we test whether extreme assumptions on
heavy-quark emission locality lead to experimentally distinguishable
predictions.

##  III Quantitative predictions and discussion

The following assumptions are made to quantify the intuitive expectations for
the considered momentum correlation:

  1. (i)

The production of charm hadrons in head-on Pb+Pb collisions is considered. The
collision energy is assumed to be adjusted to have a mean charm multiplicity
below one, allowing to neglect production of more than one $c$\- and $\bar{c}$
hadron pair in a single collision. This likely corresponds to the top CERN SPS
energy ($\sqrt{s_{NN}}\approx 17$ GeV) [41].

  2. (ii)

The charm hadrons are emitted from the freeze-out hyper-surface of a spherical
fireball undergoing a Hubble-like expansion. That is, the three velocity reads
$\vec{v}=\vec{r}/t$, with $\vec{r}={x,y,z}$ being the distance from the centre
of the fireball, and four-velocity
$u^{\mu}=x^{\mu}/\tau=x^{\mu}/\sqrt{t^{2}-r^{2}}$. It was recently
demonstrated that Hubble-like expansion is an appropriate approximation of
velocity fields in heavy-ion collisions in the energy range of our interest
[42].

  3. (iii)

The freeze-out hyper-surface is set by the freeze-out time $\tau=\tau_{fo}$
and the maximal radius $r\leq R_{max}$. They are assumed to be of the order of
the Pb nucleus radius: $\tau_{fo}=9$ fm/$c$ and $R_{max}=6$ fm.

  4. (iv)

Emission probability of charm hadrons is independent of the fluid cell on the
freeze-out hyper-surface. This is consistent with the method of predicting the
spectra within the relativistic hydrodynamics approach. Note that the
considered correlations are given by the conditional probability of the charm
hadrons to be emitted from the same cell or another one with respect to the
anti-charm one.

  5. (v)

In the rest frame of the flow, the charm hadron momentum $p$ distribution at
the freeze-out hyper-surface is assumed to be the statistical one:

$\frac{d^{3}N}{dp\leavevmode\nobreak\ d^{2}\Omega}\leavevmode\nobreak\
\propto\leavevmode\nobreak\ p^{2}\exp{\left(\frac{-\sqrt{m^{2}+p^{2}}}{T_{\rm
fo}}\right)}\leavevmode\nobreak\ ,$ (1)

where $m=1.869$ $\mbox{Ge\kern-1.00006ptV}\\!/\\!c$ is the charm hadron mass
assumed to be equal to the $D^{0}$ meson mass, and the temperature parameter
is $T_{\rm fo}$ = 150 MeV. The statistical momenta of charm hadrons are drawn
independently.

  6. (vi)

To calculate the hadron momentum in the collision rest frame, the obtained
statistical four-momentum is boosted with the flow velocity,
$\vec{v}=\vec{u}\leavevmode\nobreak\ /\leavevmode\nobreak\ \sqrt{1+u^{2}}$.

Note that for simplicity, we do not consider correlations between momenta $c$
and $\bar{c}$-hadrons resulting from the energy-momentum conservation and
dynamics of the pair creation process. The change of these correlations during
the system evolution was discussed in Refs. [4, 5] for heavy-ion collision at
the top RHIC energy ($\sqrt{s_{NN}}=200$ GeV) and at the LHC. Given that we
focus on producing charmed mesons with low $p_{\textrm{T}}$ in low-energy
collisions, we expect the back-to-back correlation will not play a significant
role in the measurement we consider in this work.

Then, the results are calculated for three different space correlations of the
$c$\- and $\bar{c}$-hadrons. These are

  1. (a)

The $c$\- and $\bar{c}$-hadrons are emitted from the same fluid cell. Thus,
the average of their momenta is set by the drift velocity of the cell. Their
actual momenta are different because of the independence of their momenta in
the fluid rest frame. This ansatz is labelled the local emission.

  2. (b)

The emission points of charm hadrons are independent of each other. Hence,
they don’t have a common drift velocity, This ansatz is labelled the
independent emission.

  3. (c)

The intermediate case is modelled, assuming the correlation function of the
emission points to be the 3D Gaussian with
$\sigma=\sigma_{x}=\sigma_{y}=\sigma_{z}=2$ fm. Note that the points are
required to be within the fireball volume. The flow components of $c$\- and
$\bar{c}$-hadrons are different but correlated, leading to the correlation of
their hadron momenta. Clearly in the limits of $\sigma\rightarrow 0$ and
$\sigma\rightarrow\infty$ one recovers the local and independent emissions,
respectively. This ansatz is labelled the correlated emission.

Figure 1 shows the distribution of $c$-$\bar{c}$ hadron pairs in the
difference of azimuthal angles $\Delta\phi$ (left) and transverse momenta
$\Delta p_{T}$ (right) for local, independent and correlated emission. The
results are obtained using the Monte Carlo technique with $10^{7}$ events
generated. The distributions of the pairs in $\Delta\phi$ significantly differ
for local, independent and correlated emissions. The differences are smaller
in the case of the transverse momentum difference. The uniform distribution in
$\Delta\phi$ for the independent emission is independent of the flow and
random momentum contributions modelling. The distributions in $\Delta\phi$
decrease monotonically for the local and correlated emission from
$\Delta\phi=0$ to $\Delta\phi=\pi$, but the quantitative properties of this
qualitative behaviour depend on model details. Nonetheless, the effect of
correlation at $\Delta\phi\approx 0$ is remarkably different compared to the
back-to-back correlations expected for charm pair production in hard parton
scatterings [4, 5]. Thus, we expect experimental data will allow
discrimination between these two different types of correlations.

Figure 1:  The distribution of $c$-$\bar{c}$ pairs in the difference of
azimuthal angles $\Delta\phi$ (left) and transverse momenta $\Delta p_{T}$
(right) for local, independent and correlated ($\sigma=2$ fm emission).

It is clear that in the case of the $\Delta\phi$ distribution rather limited
data statistics (see the next section) should allow us to distinguish between
predictions obtained assuming different space correlations between the emitted
charm hadrons and different production mechanisms of charm quarks. Encouraged
by this conclusion, we turn to the standard approach to heavy-ion collisions
[37] and, within it, discuss the implications of different possible outcomes
of the experimental measurements.

The approach pictures heavy-ion collisions at high energies as a time sequence
of the following stages:

  1. (1)

Initial stage \- a high-density quark-gluon plasma is created. QCD is assumed
to be a valid theory. Charm-anti-charm quark pairs are produced locally and in
a limited number because of the high energy threshold.

  2. (2)

Expansion stage \- the plasma expands [43], reaching the hadronisation
temperature $T_{H}\approx 150$ MeV. The (anti-)charm quarks thermalize by
interactions with the medium and acquire the medium flow.

  3. (3)

Hadronization stage \- the plasma, including the $c$ and $\bar{c}$ quarks, is
converted to hadrons and resonances following the statistical rules [44, 45]
applied in the rest frame of a plasma fluid element. Thus, the flow and
hadronization (local statistical process) contributions give the momenta of
charm hadrons.

  4. (4)

Free-streaming stage \- resonances decay, and non-interacting hadrons freely
stream in the vacuum to a detector.

Many additional details, conceptual and quantitative, can be added [37], about
the hydrodynamic evolution or the rescattering after hadronization. The stages
listed above are the most relevant to this paper.

Within the standard heavy-ion approach,

  1. (A)

The experimental data consistent with the local emission would imply a small
space separation of $c$\- and $\bar{c}$-quarks during the expansion stage.
This should be confronted with the charm-quark spatial diffusion calculated
using the QCD-inspired approaches; for recent examples, see Refs. [33, 34,
35].

  2. (B)

The experimental data consistent with the independent emission would imply a
large spatial diffusion of the charm quarks in the plasma. Ultimately, for
more accurate models, the experimental results may be inconsistent with a
local creation of heavy quarks and anti-quarks because of the limited
expansion time and maximum velocity (speed of light) in a semi-classical
transport. This might imply non-local creation or charm transport with
velocities larger than the speed of light in a high-energy-density medium. The
idea is qualitatively discussed in Appendix A.

For the historical record, this paper was motivated by the non-local,
indeterministic toy model [46] requiring the teleportation transitions in its
most symmetric version.

  3. (C)

The data consistent with the correlated emission would give a sensitive input
for restricting the charm-quark spatial diffusion in the plasma.

  4. (D)

It is always wise to leave a door open to the unexpected. The approximations
used to compute the spectra might fail to describe the correlations, and the
experimental results could qualitatively disagree with the expectations.

##  IV Feasibility of experimental measurements

This section discusses the requirements for the experimental measurements of
correlations between charm and anti-charm hadrons produced in head-on heavy-
ion collisions. The important physics condition is a mean multiplicity of
charm being small enough to neglect the production of two or more pairs of
charm hadrons. This requirement implies the measurements at relatively small
collision energies, probably close to the top SPS energy of
$\sqrt{s_{NN}}\approx 20$ GeV. It also suggests collecting data in the fixed
target mode, which, due to the Lorentz boost of the centre-of-mass allows for
high detection acceptance and efficiency. For now, we only consider
measurements of the most abundant open-charm hadrons, $D^{0}$ and
$\bar{D^{0}}$ mesons. The required statistics of recorded Pb+Pb central
collisions can be derived from the average number of reconstructed
$D^{0}\bar{D^{0}}$-pairs, $\langle
D^{0}\bar{D^{0}}\rangle_{\textnormal{rec}}$. In Appendix B, we estimate that
modern experiments at the CERN SPS should be able to record sufficient data to
measure 1000 or more $D^{0}\bar{D^{0}}$-pairs. Figure 2 demonstrates the
statistical precision of a signal from 1000 $D^{0}\bar{D^{0}}$-pairs, assuming
that the statistical fluctuations of background pairs can be neglected. We
conclude that this is sufficient to distinguish between the independent and
local emission.

Figure 2: The projection for statistical precision of measurement of the
azimuthal correlation $\Delta\phi$ assuming the experiment registered
$N=1000\leavevmode\nobreak\ D^{0}\bar{D^{0}}$ pairs. The local, independent,
and correlated emission is assumed.

The discussed measurement of correlations between $c$ and $\bar{c}$ is only
meaningful for the quarks produced as a pair. However, if multiple pairs of
$c\bar{c}$ quarks were produced within the same event, we will observe an
unavoidable background due to combining $c$\- and $\bar{c}$-hadrons
originating from different pairs. Quantifying this background suffices taking
a ratio between a multiplicity of produced $c\bar{c}$ pairs in the event to a
number of all possible combinations of $c$\- and $\bar{c}$-hadrons that could
be observed in the event. To compute this ratio for different values of
$\langle c\bar{c}\rangle$, it was assumed that $c\bar{c}$-multiplicity follows
a Poisson distribution, parameterized by the given $\langle c\bar{c}\rangle$.
This yields a probability of having more than one $c\bar{c}$ pair within the
same event, in which case one has unwanted combinations of $c$\- and
$\bar{c}$-hadrons. Obtained values of the ratio pairs to combinations are
given in Appendix B, Table 1. Values around $50\%$ indicate that about half of
the combinations of $c$\- and $\bar{c}$-hadrons are unwanted. For the real-
world analysis, background due to the misidentification of open-charm hadrons
would likely be significant. Moreover, lower values of $\langle
c\bar{c}\rangle$ also imply higher requirements for the event statistics. It
suggests that $\langle c\bar{c}\rangle$ should be selected to maximize the
results’ significance. Realistic estimates of $\langle c\bar{c}\rangle$ for
head-on Pb+Pb collisions at the top CERN SPS energies range between $0.1$ and
$1$ [41]. This further supports the conclusion that the measurements at the
CERN SPS should considered.

##  V Summary

In this work, we propose to test the locality of heavy-quark creation by
studying the momentum correlation of $c$\- and $\bar{c}$-hadrons produced in
heavy-ion collisions at collision energies with the mean multiplicity of
$c\bar{c}$ pairs below one.

We show that the azimuthal correlations of charm hadrons observed in an
experiment are sensitive to whether $c$ and $\bar{c}$ quarks are emitted
locally, close to each other in the coordinate space, or they have distant,
ultimately uncorrelated, emission points. Beyond its fundamental significance,
the locality test will provide important input for modelling charm quark
interaction with a nuclear medium.

Since the emission of multiple, uncorrelated, pairs of $c$\- and
$\bar{c}$-hadrons from a single collision would spoil the wanted correlations,
it is mandatory to perform the measurements at sufficiently low collision
energies, granting a low production probability of multiple-charm pairs. The
proposed method can also be used for bottom hadrons.

As a quantitative example, we consider charm hadron measurements in head-on
Pb+Pb collisions at the CERN SPS. Assuming typical values of data-taking
parameters for the NA61/SHINE experiment at SPS, we show that the required
measurements would need a data-taking rate of 10k Hz or more. These rates are
easily allowed by the current detector technologies. Thus, the corresponding
measurements should be possible by the upgraded NA61/SHINE and the new NA60++
experiments after the CERN LS3 upgrade period.

###### Acknowledgements.

This work is partially supported by the Polish National Science Centre grants
2018/30/A/ST2/00226, 2018/30/E/ST2/00089 and 2020/39/D/ST2/02054.

Appendices

A. Classical vs quantum-mechanical approach to heavy-quark production. Here,
we discuss why the classical approximation for charm production in the limit
of a single pair is incorrect. The simplest model to address is hydrodynamics.
The main equation is the local four-momentum conservation or, more precisely,
the expectation value of it

$0=\partial_{\mu}T^{\mu\nu}=\partial_{\mu}{\rm
tr}\left(\hat{\rho}\,\widehat{T}^{\mu\nu}\right),$ (2)

with respect to the density matrix $\hat{\rho}$ of the system. Additional
equations, the equation of state and the treatment of the non-ideal part
(transport coefficients) allow us to solve the system for the expectation
values $T^{\mu\nu}$. Statistical hadronization is then used to calculate
predictions for particle production. Scattering after the hadronization is
usually considered with a separate transport phase. In some cases, the baryon
number conservation equation

$0=\partial_{\mu}J^{\mu}_{B}=\partial_{\mu}{\rm
tr}\left(\hat{\rho}\,\widehat{J}^{\mu}_{B}\right),$ (3)

is added to the hydrodynamics equations. In principle, the electric current
and the other conserved charges should also be considered when calculating the
charge densities. Moreover, the link between tensors in space-time and
particles in phase space, necessary at the hadronization stage, is through the
relativistic Wigner distribution $W(x,p)$ [47]. The latter is the
generalization of the classical distribution function $W(x,p)\xrightarrow{\rm
classical\;limit}\propto\delta(p^{2}-m^{2})f(x,{\bf p})$. It depends on the
bilinearity of the fields and the one-particle reduced density matrix. It does
not depend on the two-particle ones and higher orders. By construction,
regardless of the ansatz (local equilibrium, viscous corrections), the
hadronization formula is for the one-particle observables only. All of the
content about particle correlations must come from somewhere else.
Despite being a very different model, similar considerations hold for the
relativistic Boltzmann equation because it stems from $W(x,p)$ too.
Relativistic kinetic theory is a limit of the evolution of the Wigner
distribution. As explained in [47], the approximations needed to use the
Relativistic Boltzmann equation instead of the more general equations for the
evolution of $W(x,p)$ include both arbitrarily small gradients and arbitrarily
weak interaction. Then, one can neglect the coupling with the two-particle
reduced density matrix, and the only “quantum leftover” is the cross-section,
which must be evaluated in the framework of axiomatic field theory. These two
conditions are enough to question whether relativistic kinetic theory can be
used for the QGP. Strong interactions and large gradients are needed to fit
the experimental data on top of a realistic (non-ideal) state equation that
already requires phenomenological modifications to the simple relativistic
Boltzmann equation. All the phenomenological modifications used in the state-
of-the-art models (temperature-dependent masses, off-shell cross sections,
etc.) do not insert any contribution from the $n$-particles reduced density
matrix. This sector of the microscopic theory is systematically neglected. The
spectra can be deduced from $W(x,p)$ alone, a one-particle object. If the
evolution of $W(x,p)$ couples mostly to itself, one can argue that these
extensions of the relativistic kinetic theory have a good chance to reproduce
the spectra (and $v_{2}$ and other one particle objects). Still, the same
cannot be said about correlations.
If one prefers a more intuitive approach to quantum fields, some
considerations must be made from first principles. Because of the Heisenberg
uncertainty, one cannot have an arbitrarily sharp wave function in both
position and momentum at the same time. The more $c\bar{c}$ pair is well-
defined in momentum, the more it must be delocalized. If the quarks are
assumed to be produced as close as possible to momentum eigenstates, to forget
about the details of the wavefunction in momentum space, one has to consider
them substantially delocalized in space. They cannot be considered in one
cell, and the wave function in the configuration space gives a weight
regarding which part of the medium is ‘’felt more” by the heavy flavours. In
any case, neither hydrodynamics nor kinetic theory are equipped to treat such
wave functions dynamically.

B. Example estimate of event statistics and data-taking time. Here, we present
a simple estimate of the event statistics and data-taking time assuming
detector setup and performance similar to the NA61/SHINE experiment at CERN
[48] recording Pb+Pb collisions at $\sqrt{s_{NN}}=17.3$ GeV. Assuming that
processes that impact the reconstruction of a $D^{0}$ and a $\bar{D^{0}}$
mesons within an event are approximately uncorrelated, we estimate the average
number of reconstructed pairs as

$\langle D^{0}\bar{D^{0}}\rangle_{rec}\approx\langle
c\bar{c}\rangle\cdot\left(P(c\rightarrow
D^{0})\cdot\textnormal{BR}(D^{0}\rightarrow K\pi)\cdot
P(\textnormal{acc})\cdot P(\textnormal{sel})\cdot
P(\textnormal{rec})\right)^{2},$ (4)

where $\langle c\bar{c}\rangle$ is the average number of $c\bar{c}$-pairs per
event. The $P(c\rightarrow D^{0})$ = 0.31 is a probability for $c$-quark to
hadronize into the $D^{0}$ meson evaluated within the PHSD model [49],
BR$(D^{0}\rightarrow K^{+}\pi^{-})$ = 3.98% is a branching ratio of decay
channel used in the measurements [50], $P(\textnormal{acc})$ = 0.5 is a
probability for $D^{0}$ to be within an acceptance region of the detector,
$P(\textnormal{sel})$ = 0.2 is a probability for $D^{0}$ to pass background-
suppressing selection of charm meson candidates, and $P(\textnormal{rec})$ =
0.9 is a probability of reconstructing the meson. The value of
$P(\textnormal{acc})$ was evaluated using the Geant4 simulation with the
detector setup for November 2022, $P(\textnormal{sel})$ is taken from the
pilot analysis of $D^{0}$ and $\bar{D^{0}}$ production [51], and
$P(\textnormal{rec})$ was obtained from a Geant4 simulation with the setup for
November 2022 and reconstruction software used for previous open charm
analysis using 2017 and 2018 data [52, 51].

Finally, given $\langle D^{0}\bar{D^{0}}\rangle_{rec}$, an estimate of the
required event statistics can be obtained via

$\\{\textnormal{number of head-on events to
collect}\\}\approx\frac{\\{\textnormal{number of {\hbox{D^0\bar{D^0}}} pairs
to reconstruct}\\}}{\langle D^{0}\bar{D^{0}}\rangle_{rec}}.$ (5)

The $\langle c\bar{c}\rangle$ value is neither reliably predicted by models
nor measured by experiments. However, considering available estimates [41], we
expect that the value of $\langle c\bar{c}\rangle$ for head-on Pb+Pb at
$\sqrt{s_{NN}}\approx 17$ GeV should range from 0.1 up to 1.

Putting all together, estimates on the run time needed to collect 1000
$D^{0}\bar{D^{0}}$-pairs for different event rates of the upgraded NA61/SHINE
experiment and for different values of $\langle c\bar{c}\rangle$ are given in
Table 1.

| $\langle c\bar{c}\rangle=0.1$ | $\langle c\bar{c}\rangle=0.2$ | $\langle c\bar{c}\rangle=0.5$ | $\langle c\bar{c}\rangle=1$
---|---|---|---|---
1 kHz |  $1000$ days |  $500$ days |  $200$ days |  $100$ days
10 kHz |  $100$ days |  $50$ days |  $20$ days |  $10$ days
100 kHz |  $10$ days |  $5$ days |  $2$ days |  $1$ day
$N_{pair}/N_{comb}$ | 91% | 83% | 66% | 50%
Table 1:  Estimate of the duration of a data-taking period needed to collect
1000 $D^{0}\bar{D^{0}}$-pairs (first three rows). In the calculations, the
duty cycle of 30% was assumed. The last row shows the ratio of the produced
pairs of $c\bar{c}$ quarks to all combinations of them.

A typical ion beam period at CERN is about four weeks. Entries in Table 1 with
a data-taking time of 100 days or more correspond to scenarios where the
measurement may take longer than a period between the CERN accelerators’ long
shutdowns. Moreover, at the moment, the event rate of 100 kHz would require a
significant upgrade of the NA61/SHINE detector and its beamline. However, a
setup corresponding to 10 kHz should be achievable within the nearest years.
Thus, we find that for $\langle c\bar{c}\rangle>0.2$, it should be possible to
perform the measurements of $c\bar{c}$-correlations by NA61/SHINE in the CERN
Run 4 period (2028-2032). The additional possibility for the experimental
study would be constructing a new experiment optimized for charm measurements.
The corresponding letter of intent was recently submitted to the CERN SPSC
[53].

## References

  * [1] E. V. Shuryak, “Quantum Chromodynamics and the Theory of Superdense Matter,” Phys. Rept. 61 (1980) 71–158.
  * [2] M. Gazdzicki and P. Seyboth, “Search for Critical Behaviour of Strongly Interacting Matter at the CERN Super Proton Synchrotron,” Acta Phys. Polon. B 47 (2016) 1201, arXiv:1506.08141 [nucl-ex].
  * [3] A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanov, and N. Xu, “Mapping the Phases of Quantum Chromodynamics with Beam Energy Scan,” Phys. Rept. 853 (2020) 1–87, arXiv:1906.00936 [nucl-th].
  * [4] X. Zhu, M. Bleicher, S. L. Huang, K. Schweda, H. Stoecker, N. Xu, and P. Zhuang, “D anti-D correlations as a sensitive probe for thermalization in high-energy nuclear collisions,” Phys. Lett. B 647 (2007) 366–370, arXiv:hep-ph/0604178.
  * [5] S. Cao, G.-Y. Qin, and S. A. Bass, “Modeling of heavy-flavor pair correlations in Au-Au collisions at 200A GeV at the BNL Relativistic Heavy Ion Collider,” Phys. Rev. C 92 no. 5, (2015) 054909, arXiv:1505.01869 [nucl-th].
  * [6] M. He and R. Rapp, “Hadronization and charm-hadron ratios in heavy-ion collisions,” Physical Review Letters 124 no. 4, (Jan, 2020) . https://doi.org/10.1103%2Fphysrevlett.124.042301.
  * [7] H. Song and U. W. Heinz, “Suppression of elliptic flow in a minimally viscous quark-gluon plasma,” Phys. Lett. B 658 (2008) 279–283, arXiv:0709.0742 [nucl-th].
  * [8] D. Kharzeev, Y. V. Kovchegov, and K. Tuchin, “Nuclear modification factor in d+Au collisions: Onset of suppression in the color glass condensate,” Phys. Lett. B 599 (2004) 23–31, arXiv:hep-ph/0405045.
  * [9] STAR Collaboration, J. Adam et al., “First Observation of the Directed Flow of $D^{0}$ and $\overline{D^{0}}$ in Au+Au Collisions at $\sqrt{s_{\rm NN}}$ = 200 GeV,” Phys. Rev. Lett. 123 no. 16, (2019) 162301, arXiv:1905.02052 [nucl-ex].
  * [10] STAR Collaboration, L. Adamczyk et al., “Measurement of $D^{0}$ Azimuthal Anisotropy at Midrapidity in Au+Au Collisions at $\sqrt{s_{NN}}$=200 GeV,” Phys. Rev. Lett. 118 no. 21, (2017) 212301, arXiv:1701.06060 [nucl-ex].
  * [11] STAR Collaboration, L. Adamczyk et al., “Elliptic flow of electrons from heavy-flavor hadron decays in Au + Au collisions at $\sqrt{s_{\rm NN}}=$ 200, 62.4, and 39 GeV,” Phys. Rev. C 95 no. 3, (2017) 034907, arXiv:1405.6348 [hep-ex].
  * [12] ALICE Collaboration, B. B. Abelev et al., “Azimuthal anisotropy of D meson production in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV,” Phys. Rev. C 90 no. 3, (2014) 034904, arXiv:1405.2001 [nucl-ex].
  * [13] ALICE Collaboration, B. Abelev et al., “D meson elliptic flow in non-central Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 2.76TeV” Phys. Rev. Lett. 111 (2013) 102301, arXiv:1305.2707 [nucl-ex].
  * [14] ALICE Collaboration, S. Acharya et al., “Elliptic Flow of Electrons from Beauty-Hadron Decays in Pb-Pb Collisions at $\sqrt{s_{NN}}$ = 5.02 TeV,” Phys. Rev. Lett. 126 no. 16, (2021) 162001, arXiv:2005.11130 [nucl-ex].
  * [15] ALICE Collaboration, S. Acharya et al., “Transverse-momentum and event-shape dependence of D-meson flow harmonics in Pb–Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV,” Phys. Lett. B 813 (2021) 136054, arXiv:2005.11131 [nucl-ex].
  * [16] ATLAS Collaboration, M. Aaboud et al., “Measurement of the suppression and azimuthal anisotropy of muons from heavy-flavor decays in Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=2.76$ TeV with the ATLAS detector,” Phys. Rev. C 98 no. 4, (2018) 044905, arXiv:1805.05220 [nucl-ex].
  * [17] ALICE Collaboration, S. Acharya et al., “$D$-meson azimuthal anisotropy in midcentral Pb-Pb collisions at $\mathbf{\sqrt{s_{\rm NN}}=5.02}$ TeV,” Phys. Rev. Lett. 120 no. 10, (2018) 102301, arXiv:1707.01005 [nucl-ex].
  * [18] ALICE Collaboration, S. Acharya et al., “Elliptic Flow of Electrons from Beauty-Hadron Decays in Pb-Pb Collisions at $\sqrt{s_{NN}}$ = 5.02 TeV,” Phys. Rev. Lett. 126 no. 16, (2021) 162001, arXiv:2005.11130 [nucl-ex].
  * [19] PHENIX Collaboration, A. Adare et al., “Heavy Quark Production in $p+p$ and Energy Loss and Flow of Heavy Quarks in Au+Au Collisions at $\sqrt{s_{NN}}=200$ GeV,” Phys. Rev. C 84 (2011) 044905, arXiv:1005.1627 [nucl-ex].
  * [20] ALICE Collaboration, S. Acharya et al., “First study of the two-body scattering involving charm hadrons,” Phys. Rev. D 106 no. 5, (2022) 052010, arXiv:2201.05352 [nucl-ex].
  * [21] F. Grosa, “ALICE determines the scattering parameters of D mesons with light-flavor hadrons,” 2022.  Quark Matter 2022.
  * [22] STAR Collaboration, J. Adam et al., “Centrality and transverse momentum dependence of $D^{0}$-meson production at mid-rapidity in Au+Au collisions at ${\sqrt{s_{\rm NN}}=\rm{200\,GeV}}$,” Phys. Rev. C 99 no. 3, (2019) 034908, arXiv:1812.10224 [nucl-ex].
  * [23] STAR Collaboration, L. Adamczyk et al., “Measurements of $D^{0}$ and $D^{*}$ Production in $p+p$ Collisions at $\sqrt{s}=200$ GeV,” Phys. Rev. D 86 (2012) 072013, arXiv:1204.4244 [nucl-ex].
  * [24] ALICE Collaboration, S. Acharya et al., “Charm-quark fragmentation fractions and production cross section at midrapidity in pp collisions at the LHC,” Phys. Rev. D 105 no. 1, (2022) L011103, arXiv:2105.06335 [nucl-ex].
  * [25] ALICE Collaboration, S. Acharya et al., “Prompt D0, D+, and D∗+ production in Pb–Pb collisions at $\sqrt{s_{\mathrm{NN}}}$ = 5.02 TeV,” JHEP 01 (2022) 174, arXiv:2110.09420 [nucl-ex].
  * [26] A. Andronic et al., “Heavy-flavour and quarkonium production in the LHC era: from proton–proton to heavy-ion collisions,” Eur. Phys. J. C 76 no. 3, (2016) 107, arXiv:1506.03981 [nucl-ex].
  * [27] M. Attems, J. Brewer, G. M. Innocenti, A. Mazeliauskas, S. Park, W. van der Schee, and U. A. Wiedemann, “The medium-modified $g\to c\overline{c}$ splitting function in the BDMPS-Z formalism,” JHEP 01 (2023) 080, arXiv:2203.11241 [hep-ph].
  * [28] M. Attems, J. Brewer, G. M. Innocenti, A. Mazeliauskas, S. Park, W. van der Schee, and U. Wiedemann, “Medium-enhanced $c\bar{c}$ radiation,” arXiv:2209.13600 [hep-ph].
  * [29] M. Rohrmoser, Study of correlations of heavy quarks in heavy ion collisions and their role in understanding the mechanisms of energy loss in the quark gluon plasma.  PhD thesis, SUBATECH, Nantes, 2017.
  * [30] M. Nahrgang, J. Aichelin, P. B. Gossiaux, and K. Werner, “Azimuthal correlations of heavy quarks in pb + pb collisions at $\sqrt{s}=2.76$ tev at the cern large hadron collider,” Phys. Rev. C 90 (Aug, 2014) 024907. https://link.aps.org/doi/10.1103/PhysRevC.90.024907.
  * [31] H. Wang, J.-H. Chen, Y.-G. Ma, and S. Zhang, “Charm hadron azimuthal angular correlations in Au + Au collisions at $\sqrt{s_{\mathrm{NN}}}=200$ GeV from parton scatterings.” 2019.
  * [32] H. Wang and J.-H. Chen, “Study on open charm hadron production and angular correlation in high-energy nuclear collisions,” Nucl. Sci. Tech. 32 no. 1, (2021) 2.
  * [33] M. L. Sambataro, S. Plumari, and V. Greco, “Impact of off-shell dynamics on the transport properties and the dynamical evolution of Charm Quarks at RHIC and LHC temperatures,” Eur. Phys. J. C 80 no. 12, (2020) 1140, arXiv:2005.14470 [hep-ph].
  * [34] F. Capellino, A. Beraudo, A. Dubla, S. Floerchinger, S. Masciocchi, J. Pawlowski, and I. Selyuzhenkov, “Fluid-dynamic approach to heavy-quark diffusion in the quark-gluon plasma,” Phys. Rev. D 106 no. 3, (2022) 034021, arXiv:2205.07692 [nucl-th].
  * [35] S. Satapathy, S. De, J. Dey, and S. Ghosh, “Spatial diffusion of heavy quarks in background magnetic field,” arXiv:2212.08933 [hep-ph].
  * [36] ATLAS Collaboration, “Measurements of the azimuthal correlations of muon pairs from heavy-flavor decays in 5.02 TeV Pb+Pb collisions with the ATLAS detector,” tech. rep., CERN, Geneva, 2022.  ATLAS-CONF-2022-022.
  * [37] W. Florkowski, Phenomenology of Ultra-relativistic Heavy-ion Collisions.  World Scientific, 2010.
  * [38] L. Tinti, “Quantum free-streaming: out of equilibrium expansion for the free scalar fields,” arXiv:2304.00109 [hep-ph].
  * [39] J. Jankowski and M. Spaliński, “Hydrodynamic attractors in ultrarelativistic nuclear collisions,” Prog. Part. Nucl. Phys. 132 (2023) 104048, arXiv:2303.09414 [nucl-th].
  * [40] L. Tinti, “Off-shell hydrodynamic expansion,” Phys. Rev. D 108 no. 3, (2023) 036015.
  * [41] NA61/SHINE Collaboration, A. Snoch, “Charm Program of NA61/SHINE: Motivation and Measurements,” arXiv:1803.01692 [nucl-ex].
  * [42] N. S. Tsegelnik, E. E. Kolomeitsev, and V. Voronyuk, “Helicity and vorticity in heavy-ion collisions at energies available at the JINR Nuclotron-based Ion Collider facility,” Phys. Rev. C 107 no. 3, (2023) 034906, arXiv:2211.09219 [nucl-th].
  * [43] L. D. Landau, “On the multiparticle production in high-energy collisions,” Izv. Akad. Nauk Ser. Fiz. 17 (1953) 51–64.
  * [44] R. Hagedorn and J. Rafelski, “Hot Hadronic Matter and Nuclear Collisions,” Phys. Lett. B 97 (1980) 136.
  * [45] F. Becattini, J. Manninen, and M. Gazdzicki, “Energy and system size dependence of chemical freeze-out in relativistic nuclear collisions,” Phys. Rev. C 73 (2006) 044905, arXiv:hep-ph/0511092.
  * [46] M. Gazdzicki, M. Gorenstein, I. Pidhurskyi, O. Savchuk, and L. Tinti, “Equilibration and Locality,” Acta Phys. Polon. B 53 no. 8, (2022) 2, arXiv:2206.01151 [nucl-th].
  * [47] S. R. De Groot, Relativistic Kinetic Theory. Principles and Applications.  1980\.
  * [48] NA61 Collaboration, N. Abgrall et al., “NA61/SHINE facility at the CERN SPS: beams and detector system,” JINST 9 (2014) P06005, arXiv:1401.4699 [physics.ins-det].
  * [49] W. Cassing and E. Bratkovskaya, “Parton–hadron–string dynamics: An off-shell transport approach for relativistic energies,” Nuclear Physics A 831 no. 3-4, (Dec, 2009) 215–242. https://doi.org/10.1016%2Fj.nuclphysa.2009.09.007.
  * [50] Particle Data Group Collaboration, R. L. Workman and Others, “Review of Particle Physics,” PTEP 2022 (2022) 083C01.
  * [51] A. Merzlaya, Open charm measurements at the NA61/SHINE experiment at CERN SPS with the new Vertex Detector.  PhD thesis, Jagiellonian University, 2021.  https://cds.cern.ch/record/2771816.  Presented 07 Jun 2021.
  * [52] M. L. Bajda, “Development of track reconstruction for the upgraded NA61/SHINE vertex detector,” Master’s thesis, Jagiellonian University, 2023\.
  * [53] NA60+ Collaboration, C. Ahdida, G. Alocco, et al., “Letter of Intent: the NA60+ experiment,” tech. rep., CERN, Geneva, 2022.  https://cds.cern.ch/record/2845241.


# Dynamic mode decomposition as an analysis tool for time-dependent partial
differential equations

Miha Rot12, Martin Horvat3, Gregor Kosec2 The authors would like to
acknowledge the financial support of the Slovenian Research Agency (ARRS)
research core funding No. P2-0095 and the Young Researcher program PR-10468.
1Jozef Stefan International Postgraduate School, Jamova cesta 39, 1000
Ljubljana, Slovenia  2“Jožef Stefan” Institute, Parallel and Distributed
Systems Laboratory, Jamova cesta 39, 1000 Ljubljana, Slovenia  3University of
Ljubljana, Faculty of Mathematics and Physics, Department of Physics,
Jadranska cesta 19, SI-1000 Ljubljana.
<EMAIL_ADDRESS>
###### Abstract

The time-dependent fields obtained by solving partial differential equations
in two and more dimensions quickly overwhelm the analytical capabilities of
the human brain. A meaningful insight into the temporal behaviour can be
obtained by using scalar reductions, which, however, come with a loss of
spatial detail. Dynamic Mode Decomposition is a data-driven analysis method
that solves this problem by identifying oscillating spatial structures and
their corresponding frequencies. This paper presents the algorithm and
provides a physical interpretation of the results by applying the
decomposition method to a series of increasingly complex examples.

###### Index Terms:

dynamic mode decomposition; data-driven analysis; wave equation; non-Newtonian
natural convection

##  I Introduction

Humans posses a remarkable ability for pattern recognition that is however
still quickly overwhelmed when analysing complex systems. When dealing with
time series of two or three dimensional field snapshots, as is common in e.g.
hydrodynamics, we are often forced to use reductions that decrease complexity
and provide an approachable insight into the system’s behaviour. An example of
such reduction when dealing with natural convection is tracking the ratio
between heat convection and conduction, also known as the Nusselt number, on a
significant boundary. The Nusselt number provides us with a scalar value that
reflects the behaviour in our system. Analysing time series of such values can
provide insight on whether the system is stable, oscillatory or maybe even
chaotic, with Fourier analysis available as a potent tool to extract the
characteristic frequencies reflected in the chosen reduction. No matter how
good the reduction is we still lose the majority of the spatial detail which
is often important for understanding the dynamics.

Spatial structure of the analysed system is required for good understanding of
the dynamics, which is especially important when dealing with engineering
challenges like the airflow around aircraft structures [1]. A common approach
to identifying important spatial structures in analysed fields is modal
decomposition, with Proper Orthogonal Decomposition (POD)111Also known as
principal component analysis (PCA) and various other names in different
fields., proposed for hydrodynamics by Lumley [2] in 1967, and it’s variations
[3, 4] as one of the most widely used techniques. POD is a decomposition
technique that identifies orthogonal modes that best222Optimally in $L_{2}$
sense that corresponds to energy when decomposing flow field [5]. represent
the dataset and provides an ordering based on their importance.

Modes obtained with POD provide information about energetically important
parts of the system but are completely oblivious to the temporal behaviour
meaning that the results of the decomposition would stay the same even if
analysis was performed on reordered snapshots. Dynamic mode decomposition
(DMD) was proposed by Schmid in a 2008 talk and the subsequent paper in 2010
[6] as a technique that joins the spatial aspects of POD and temporal aspects
of Fourier transform [4]. DMD identifies characteristic frequencies,
corresponding spatial structures and whether they amplify, decay or remain
constant through the sampled timespan. The method was initially proposed for
hydrodynamics but has since found other applications including infectious
disease spread [7] and computer vision [8] as it is completely data-driven and
independent of the underlying dynamics.

This paper introduces the DMD algorithm and interpretation in Sec. II,
demonstrates the results on vibrating membrane examples in Sec. III and
concludes with the decomposition applied to a more interesting case of
oscillatory non-Newtonian natural convection in Sec. IV.

##  II Dynamic Mode Decomposition

###  II-A Interpretation

The DMD algorithm is based on finding eigenvalues $\lambda$ and eigenvectors
$\boldsymbol{\varphi}$ of the linear mapping represented by the matrix
$\boldsymbol{A}$

$\boldsymbol{v}_{i+1}=\boldsymbol{A}\boldsymbol{v}_{i},$ (1)

that connects the subsequent states $\boldsymbol{v}_{i}$ and
$\boldsymbol{v}_{i+1}$ of the analysed system. As such the algorithm is
completely independent of the governing equations that drive the underlying
dynamics and can be applied to any experimental or simulated data.

The linear nature of eigendecomposition does not preclude us from analysing
non-linear systems as the dynamics of any such system can be expressed with an
infinite-dimensional linear Koopman operator [9] with the DMD providing a good
approximation of it’s eigendecomposition as long as the quality and quantity
of input data is sufficient [10].

DMD eigenvalues $\lambda$ provide information about the temporal behaviour of
the spatial structure described in the eigenvector. Eigenvalues are complex
numbers with the mode’s angular frequency information contained in the
argument

$\nu_{i}=\frac{\arg(\lambda_{i})}{2\pi\Delta t},$ (2)

with sampling interval $\Delta t$ providing scaling to the output, that is
otherwise unaffected by the frequency of data snapshot sampling. The magnitude
of eigenvalues provides information about changes in mode’s strength, with
$|\lambda|>1$ for amplifying and $|\lambda|<1$ for decaying modes.

The linear combination of DMD modes calculated from real valued input data has
to be real valued, which can be satisfied either by real
eigenvalues/eigenvectors or complex-conjugate pairs of complex
eigenvalues/eigenvectors. Real eigenvalues represent modes with frequency 0.
At least one such background mode with frequency 0 and eigenvalue magnitude of
1 is present whenever the decomposed data has non-zero average but multiple
such real-valued growing/decaying non-oscillatory modes can appear. Complex
conjugate pairs represent oscillating modes and can be treated as a single
mode and will be in most of the later visualisations.

###  II-B Algorithm

The algorithm for exact DMD proposed by Tu [11] starts by taking subsequent
snapshots $\boldsymbol{v}_{i}$ of analysed fields separated by $\Delta t$.
Uniformly sampled snapshots simplify calculation and frequency interpretation
but are not explicitly required for DMD. Snapshots are flattened from their
arbitrary shape into vectors with length $M$ and arranged as columns in a
matrix

$\boldsymbol{V}_{i}^{N}=\matrixquantity[\boldsymbol{v}_{i}&\boldsymbol{v}_{i+1}&\cdots&\boldsymbol{v}_{N}],$
(3)

with $N$ denoting the number of snapshots. The number of data points in a
snapshot $M$ is usually much larger than the number of snapshots $N$, leading
to tall and narrow matrices $\boldsymbol{V}_{i}^{N}$ that can be exploited to
efficiently calculate the eigen-decomposition of the large $M\times M$ matrix
$\boldsymbol{A}$.

The snapshot matrix $\boldsymbol{V}_{i}^{N}$ needs to have a
sufficient333Generally at least two times the number of relevant oscillatory
modes. rank to describe the dynamics. Rank deficiency issues usually appear
when dealing with standing waves [10, 12] that cause a linear dependency
between snapshots in columns of $\boldsymbol{V}_{i}^{N}$ and can be solved by
stacking subsequent snapshots

$\boldsymbol{V}_{i}^{N}=\matrixquantity[\boldsymbol{v}_{i}&\boldsymbol{v}_{i+1}&\cdots&\boldsymbol{v}_{N-m}\\\
\boldsymbol{v}_{i+1}&\boldsymbol{v}_{i+2}&\cdots&\boldsymbol{v}_{N-m+1}\\\
\vdots&\vdots&\ddots&\vdots\\\
\boldsymbol{v}_{i+m}&\boldsymbol{v}_{i+m+1}&\cdots&\boldsymbol{v}_{N}],$ (4)

ensuring $\rank(\boldsymbol{V}_{i}^{N})>=m$ at the cost of effectively larger
system, and only using the first $M$ values from the now larger DMD
eigenvectors. Time shift augmentation can also help with noisy data or in
cases where snapshots are small [12].

For simplicity we continue our discussion with the notation from Eq. (3),
apply the relationship from Eq. (1) to the matrix

$\boldsymbol{V}_{1}^{N}=\matrixquantity[\boldsymbol{A}\boldsymbol{v}_{0}&\boldsymbol{A}\boldsymbol{v}_{i}&\cdots&\boldsymbol{A}\boldsymbol{v}_{N-1}]=\boldsymbol{A}\boldsymbol{V}_{0}^{N-1},$
(5)

and perform the compact singular value decomposition (SVD) on the data matrix

$\boldsymbol{V}_{0}^{N-1}=\boldsymbol{U}\boldsymbol{S}\boldsymbol{W}^{\dagger}.$
(6)

The results of the SVD can already be assigned physical meaning as it is one
of the two main algorithms used to compute POD. Singular values denoting the
importance of POD modes are contained in the diagonal $\boldsymbol{S}$ with
the corresponding eigenvectors in columns of $\boldsymbol{U}$. The complexity
of the underlying dynamics can be estimated from the spectrum of singular
values. Cases with exponentially decreasing singular value magnitudes signify
a low-dimensional behaviour that can be described with a handful of modes.
Slow decay is not as simple to explain as it can be a product of either very
complex dynamics or noisy data.

Output of the SVD algorithm can be truncated based on singular values or prior
knowledge about the system to retain only $r$ singular values and
corresponding singular vectors in $\boldsymbol{U}$ and $\boldsymbol{W}$.
Truncation is beneficial as it reduces the number of numerical operations, as
evident from the following paragraph, and removes noise, but care has to be
taken in order to avoid truncating relevant parts of the dynamics. More
details about truncation with examples and advanced strategies are available
in [12].

The optionally truncated result can be reorganized and rewritten as

$\boldsymbol{U}^{\dagger}\boldsymbol{V}_{1}^{N}\boldsymbol{W}\boldsymbol{S}^{-1}=\boldsymbol{U}^{\dagger}\boldsymbol{A}\boldsymbol{U}:=\widetilde{\boldsymbol{A}},$
(7)

with $\widetilde{\boldsymbol{A}}$ representing the projection of
$\boldsymbol{A}$ onto the POD eigenvectors. We then perform the
eigendecomposition of $\widetilde{\boldsymbol{A}}$ which is a much smaller
$r\times r$ matrix than the full, $M\times M$ sized, matrix $\boldsymbol{A}$.
The decomposition

$\widetilde{\boldsymbol{A}}\boldsymbol{w}=\lambda\boldsymbol{w}.$ (8)

gives us eigenvalues $\lambda_{i}$ and eigenvectors $\boldsymbol{w}_{i}$. The
eigenvalues are the same for both $\widetilde{\boldsymbol{A}}$ and
$\boldsymbol{A}$ while the eigenvectors $\boldsymbol{\varphi}$ for the latter
still need to be transformed from their projected representation

$\boldsymbol{\varphi}=\frac{1}{\lambda}\boldsymbol{V}_{1}^{N}\boldsymbol{W}\boldsymbol{S}^{-1}\boldsymbol{w}.$
(9)

###  II-C Mode ordering

The main disadvantage of DMD is that the hierarchy of the identified modes is
unclear. Their relative importance can not be deduced from either the unitary
eigenvectors or the eigenvalues that carry temporal data. This deficiency is
bypassed by calculating the projection coefficients for expressing the
measured datasets in terms of DMD eigenvectors

$\boldsymbol{v}(t)=\sum_{k}\boldsymbol{\varphi}_{k}b_{k}(t).$ (10)

As the relative importance of eigenvectors in an oscillating dataset changes
with time it is best to calculate their importance for all available snapshots
expressed in matrix form as

$\boldsymbol{V}_{0}^{N-1}=\boldsymbol{\Phi}\boldsymbol{B},$ (11)

with eigenvectors collected in columns of matrix
$\boldsymbol{\Phi}=\matrixquantity[\boldsymbol{\varphi}_{0}&\boldsymbol{\varphi}_{1}&\cdots&\boldsymbol{\varphi}_{r}]$
and coefficients corresponding to datasets in $\boldsymbol{V}_{0}^{N-1}$
collected as columns in matrix
$\boldsymbol{B}=\matrixquantity[\boldsymbol{b}_{0}&\boldsymbol{b}_{1}&\cdots,\boldsymbol{b}_{N-1}]$.
The complex number $\boldsymbol{B}_{ij}$ corresponds to the projection
coefficient of the $i$-th eigenvector for the $j$-th data sample. Columns of
$\boldsymbol{B}$ can be thus interpreted as the time evolution of mode
importance that can be condensed into a single value by calculating the
average value. We use the mode power ($P$)

$P_{i}=\frac{1}{N}\sum_{j=0}^{N-1}\absolutevalue{B_{ij}}$ (12)

to rank the calculated modes.

An efficient algorithm for calculating $\boldsymbol{B}$ has also been proposed
by Tu et al. [10]. The algorithm avoids additional computation beyond
calculating the left eigenvectors $\boldsymbol{Z}$ during the
eigendecomposition of $\widetilde{\boldsymbol{A}}$ by expressing the matrix as

$\boldsymbol{B}=\boldsymbol{Z}^{\dagger}\boldsymbol{S}\boldsymbol{W}^{\dagger}.$
(13)

##  III Vibrating membrane

###  III-A Numerical model

We demonstrate the DMD algorithm with a damped wave equation on a 2D square
membrane. This case is convenient as it still allows for comparison against
known analytic results but is closer to possible real world applications. We
are solving the wave equation

$\partialderivative[2]{h}{t}=c^{2}\laplacian
h-\gamma\partialderivative{h}{t},$ (14)

with $c$ denoting the velocity of wave propagation and $\gamma$ the damping
factor, on a square domain $\Omega=\\{(x,y):x,y\in[0,1]\\}$ with Dirichlet
boundary condition $h_{\partial\Omega}=0$. A constant velocity $c=1$ is used
for all of the presented examples.

We start the dynamics with a Gaussian initial condition with $\sigma=0.1$ and
amplitude of $1$ placed at an arbitrarily chosen position $(0.3,0.4)$ to
simulate an impact against the membrane. We would like to show that DMD
identifies eigenmodes for the membrane and analyse how damping might affect
the accuracy.

###  III-B Analytic solution

We can obtain an analytic solution for the partial differential Eq. (14) by
using the separation of variables. For now we ignore damping as it is not
relevant for the eigenstate spatial configuration. The solution is a linear
combination of modes defined with two integer quantization indices $m,n>=1$.
Analytical eigenfrequencies and eigenvectors are expressed in terms of $m$ and
$n$

$\nu=\frac{c}{2}\sqrt{m^{2}+n^{2}},$ (15) $h(x,y,m,n)=\sin(m\pi x)\sin(n\pi
y),$ (16)

with eigenvectors shown in Figure 1. Only one eigenvector is shown for the
degenerate frequencies.

Figure 1: Analytic eigenvectors for the oscillation of a square membrane,
ordered by increasing frequency.

###  III-C Results

The partial differential equation is solved numerically for different damping
rates using the meshless RBF-FD [13] method implemented with the Medusa [14]
C++ library. A subset of system state snapshots for the 0 damping case that
that we use in DMD is shown in the right graphs of Figure 2. The dynamics are
disordered due to the non-symmetric position of the initial jolt and do not
offer much insight at the first glance. We can use the sum of membrane heights
as an observable that reduces the system’s behaviour into a single scalar
value to visualise the oscillation and the effect of different damping rates
as shown in the left graph of Figure 2.

Figure 2: left: The sum of height values provides a simple insight into the
system’s behaviour. The time evolution of thus derived observable is shown for
different damping rates $\gamma$. right: A selection of snapshots used in the
decomposition providing an insight into the system’s behaviour during the
observed time. A heat map display is used to visualise the scalar values with
blues for negative, white for zero and reds for positive values.

We use DMD on 1000 uniformly sampled snapshots from the interval $t\in[0,10]$
that are composed into the data matrix $\boldsymbol{V}$ with stacking $m=10$.
Each snapshot is composed of 3985 height values in computational nodes for the
mesh-free partial differential equation solving, which are uniformly
positioned within $\Omega$.

The resulting frequencies, calculated from DMD eigenvalues as described in Eq.
(2), are shown in the top graph of Figure 3. The results provide a good
estimation for the analytic results with sub percent relative error for the
first $\sim 30$ modes as shown in the central graph. Estimated frequencies
remain good even for the relatively strongly damped $\gamma=1$ case showing
promise for eigenfrequency estimation on experimental data obtained from
complex objects.

There is a slight inconsistency in the results with a sharp increase in the
relative error for the $\gamma=0$ frequencies. We would expect that the least
damped case would give the best results but there seems to be an additional
non-physical mode at index 11 that has caused a shift in the following
frequencies. Fortunately we can use the mode power described in Sec. II-C and
shown in the bottom graph of Figure 3 to detect and remove the highlighted
spurious mode with a significantly lower strength. It is surprising how little
nonsensical modes we identified even though we initially disturbed the
membrane in a non-symmetric location and how chaotic the actual system shown
in the heat maps of Figure 2 looks to a human observer.

Figure 3: top: Eigenfrequencies calculated with DMD for different damping
rates compared to analytic values. centre: The relative error between DMD
eigenfrequencies and analytic values for the square membrane. bottom: Mode
powers for the calculated DMD modes.

We can also compare the analytic eigenvectors shown in Figure 1 and DMD
eigenvectors in Figure 4. The DMD vector positions in the grid are shifted by
one compared to the analytic due to the non-oscillatory background mode that
is always present in the decomposition of data with non-zero average. Some
eigenvectors match almost perfectly while others wildly differ. This might
seem inconsistent at the first glance but can be explained with degeneracy in
the analytical spectrum. Linear combinations of degenerate eigenvectors still
constitute a valid eigenvector. The non-symmetric eigenvector $\varphi_{11}$
can quickly be identified as the one corresponding to the previously mentioned
spurious mode.

Figure 4: DMD mode eigenvectors ordered by frequency for the $\gamma=0$ case.
The positions of modes are shifted by one relative to Figure 3 because of the
constant, zero frequency, background mode. The grayscale image displays
absolute value with contours ranging from blue to red for the imaginary part.

###  III-D Oddly shaped membrane

Now that we have shown that DMD can be used to identify inherent oscillatory
modes we apply it to an irregularly shaped membrane. The duck-shaped membrane
is again disturbed with a Gaussian initial state, resulting in the oscillation
shown in Figure 5. This combination of meshless computational nodes, the time
evolution of the system, the DMD spectrum, and the DMD eigenvectors ordered by
frequency provides a condensed overview of the dynamics and shows how DMD
could be used in engineering analysis. The DMD spectrum is a convenient and
commonly used display of DMD mode frequency and power in the same graph, which
can, when sorted by frequency, be used similarly to the Fourier transform.

Figure 5: The DMD decomposition applied to a duck-shaped oscillating membrane.
top left: Positioning of computational nodes with yellow for interior nodes
where we solve the wave equation and purple for boundary nodes with enforced
Dirichlet condition. top right: A subset of membrane state snapshots used for
DMD. centre: The DMD spectrum showing the frequencies and powers of modes
ordered by the former. bottom: Eigenvectors of DMD modes ordered by frequency.

##  IV Non-Newtonian fluid

The last example is a practical use case for DMD applied to hydrodynamics. We
solve, with details available in [15], the system of partial differential
equations

$\displaystyle\divergence\boldsymbol{v}$ $\displaystyle=0,$ (17)
$\displaystyle\rho(\partialderivative{\boldsymbol{v}}{t}+\boldsymbol{v}\cdot\gradient{\boldsymbol{v}})$
$\displaystyle=-\gradient
p+\divergence(\eta\gradient\boldsymbol{v})-\boldsymbol{g}\rho\beta
T_{\Delta},$ (18) $\displaystyle\rho
c_{p}(\partialderivative{T}{t}+\boldsymbol{v}\cdot\gradient{T})$
$\displaystyle=\divergence(\lambda\gradient T),$ (19) $\displaystyle\eta$
$\displaystyle=\eta_{0}\left(\frac{1}{2}\norm{\gradient{\boldsymbol{v}}+(\gradient{\boldsymbol{v}})^{T}}\right)^{\frac{n-1}{2}}.$
(20)

with $\boldsymbol{v}$, $T$, $p$, $\rho$, $\boldsymbol{g}$, $\beta$,
$T_{\Delta}$, $c_{p}$, $\eta_{0}$, $n$ representing the flow velocity field,
temperature field, pressure field, density, gravity, thermal expansion
coefficient, temperature offset, heat capacity, viscosity constant and non-
Newtonian index respectively.

The system is used to describe the natural convection in an incompressible
non-Newtonian fluid. Non-Newtonian fluids have a variable viscosity meaning
that the relationship between the shear-strain and the shear-stress is non-
linear. We use a simple power-law model to describe the relationship with non-
Newtonian power index $n$ as the main parameter of behaviour. Smaller values
of $n$ indicate stronger shear-thinning non-Newtonian effects while $n=1$
equals to a normal, Newtonian, fluid.

We will be using the dimensionless Rayleigh number (Ra) that has a similar
meaning for natural convection as Reynolds number, meaning stronger dynamics
for higher values, as a case parameter, with details unimportant for the
demonstration.

The model is solved in a square cavity with differentially heated walls
schematically shown in the upper left graph of Figure 6. The left and right
walls are kept at a temperature differential causing the fluid to form a
vortex as it heats and rises at the right boundary and cools and descends on
the left. This circulation is stationary for low Ra with a transition into an
oscillatory and later chaotic regime as Ra increases. The upper right graphs
in Figure 6 show a selection of system temperature and velocity field
snapshots for the oscillatory behaviour at Ra$=10^{6}$. A set of parameters
where the transition into the oscillatory behaviour occurs has been identified
by Kosec et al. [16] and we use DMD to analyse the differences in dynamics as
we push past that point with increasing Ra and decreasing $n$.

Figure 6: The setup and behaviour for the differentially heated cavity non-
Newtonian natural convection case. upper left:  The case setup with boundary
conditions. upper right: A selection of snapshots sampled from the system with
Ra$=10^{6}$ and $n=1$. The red-blue background colour shows the local fluid
temperature with added streamlines to illustrate the flow direction. bottom:
Time evolution of Nusselt number for different parameters.

As mentioned in the introduction we can use the Nusselt number, shown in the
bottom graph of Figure 6, as an observable into the systems behaviour. The
system is stationary for Ra$=5\cdot 10^{3}$ and $n=1$ but oscillations occur
both when we increase Ra and when we decrease $n$.

We use DMD to decompose the 5 cases utilising 500 snapshots of the system
between $t=10$ and $t=11$ with stacking $m=5$. The results are shown in Figure
7a for increasing Ra and Figure 7b for decreasing $n$. The mode power defined
in Eq. (12) is used to identify the strongest modes and display the
corresponding eigenvectors. The first observation that is common to both modes
of spurring the dynamics pertains to the DMD spectrum. The initial stationary
case has the vast majority of its power in the first constant mode with others
most likely only containing noise. As the intensity of the dynamics increases
we get more distinct modes but the power of the background modes also
increases as expected with ever wilder dynamics heading towards chaos, where
all modes would be present and indistinguishable by power.

(a) Different Ra, Newtonian $n=1$.

(b) Different non-Newtonian $n$, Ra$=5\cdot 10^{3}$.

Figure 7: The DMD spectra and eigenvectors for the strongest modes.

We can be reasonably certain that the oscillatory dynamics are in fact
different for the Ra and $n$ spurred dynamics by using the eigenvectors of DMD
modes sorted by power to identify the important parts. The strongest areas are
clearly different with eigenvectors for the Ra induced dynamics in Figure 7a
present mainly as perturbations to the main central vortex while the modes for
the $n$ induced dynamics in Figure 7b point to activity in corners where
counter vortices occur. This is consistent with the typical shear-thinning
non-Newtonian behaviour where the larger velocity gradients are penalised less
providing better conditions for smaller vortices.

##  V Conclusion

We have presented the algorithm for dynamic mode decomposition and applied it
to various test and example cases. The examples progressed from the most
basic, where we were able to verify the results agains known closed-form
values, towards a concrete hydrodynamics example where DMD provided valuable
insight into the system that would be difficult to obtain otherwise. The DMD
has proved to be a very practical tool for dynamical system analysis and it is
recommended that anyone dealing with hydrodynamics or similarly complex
systems should be at least familiar with its fundamentals. The algorithm is
only a decade old and has an active community working on extensions for the
concept of Koopman operator decomposition and applications for the existing
versions of DMD.

## References

  * [1] T. R. Ricciardi, W. Wolf, R. L. Speth, and P. Bent, “Analysis of noise sources in realistic landing gear configurations through high fidelity simulations,” in _AIAA Scitech 2019 Forum_ , 2019. [Online]. Available: https://arc.aiaa.org/doi/abs/10.2514/6.2019-0003
  * [2] J. L. Lumley, “The structure of inhomogeneous turbulent flows,” in _Proceedings of the International Colloquium on the Fine Scale Structure of the Atmosphere and its Influence on Radio Wave Propagation_ , A. M. Yaglam and V. I.Tatarsky, Eds., Moscow, 1967, pp. 166–178.
  * [3] C. W. Rowley, “Model reduction for fluids, using balanced proper orthogonal decomposition,” _International Journal of Bifurcation and Chaos_ , vol. 15, no. 03, pp. 997–1013, 2005.
  * [4] K. Taira, S. L. Brunton, S. T. M. Dawson, C. W. Rowley, T. Colonius, B. J. McKeon, O. T. Schmidt, S. Gordeyev, V. Theofilis, and L. S. Ukeiley, “Modal analysis of fluid flows: An overview,” _AIAA Journal_ , vol. 55, no. 12, pp. 4013–4041, 2017. [Online]. Available: https://doi.org/10.2514/1.J056060
  * [5] P. Holmes, J. L. Lumley, G. Berkooz, and C. W. Rowley, _Turbulence, Coherent Structures, Dynamical Systems and Symmetry_ , 2nd ed., ser. Cambridge Monographs on Mechanics. Cambridge University Press, 2012.
  * [6] P. J. Schmid, “Dynamic mode decomposition of numerical and experimental data,” _Journal of Fluid Mechanics_ , vol. 656, p. 5–28, 2010.
  * [7] J. Proctor and P. Welkhoff, “Discovering dynamic patterns from infectious disease data using dynamic mode decomposition,” _International health_ , vol. 7, pp. 139–45, 03 2015.
  * [8] J. Grosek and J. N. Kutz, “Dynamic mode decomposition for real-time background/foreground separation in video,” _arXiv e-prints_ , p. arXiv:1404.7592, Apr. 2014.
  * [9] C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. S. Henningson, “Spectral analysis of nonlinear flows,” _Journal of Fluid Mechanics_ , vol. 641, p. 115–127, 2009.
  * [10] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. N. Kutz, “On dynamic mode decomposition: Theory and applications,” _Journal of Computational Dynamics_ , vol. 1, no. 2, pp. 391–421, 2014.
  * [11] J. Tu, “Dynamic mode decomposition: Theory and applications,” Ph.D. dissertation, Princeton University, sep 2013.
  * [12] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor, _Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems_ , ser. Other Titles in Applied Mathematics. Society for Industrial and Applied Mathematics, 2016. [Online]. Available: https://books.google.si/books?id=wbedDAEACAAJ
  * [13] A. I. Tolstykh and D. A. Shirobokov, “On using radial basis functions in a “finite difference mode” with applications to elasticity problems,” _Computational Mechanics_ , vol. 33, no. 1, pp. 68–79, Dec 2003. [Online]. Available: https://doi.org/10.1007/s00466-003-0501-9
  * [14] J. Slak and G. Kosec, “Medusa: A c++ library for solving pdes using strong form mesh-free methods,” _ACM Trans. Math. Softw._ , vol. 47, no. 3, Jun. 2021. [Online]. Available: https://doi.org/10.1145/3450966
  * [15] M. Rot and G. Kosec, “Refined rbf-fd analysis of non-newtonian natural convection,” 2022. [Online]. Available: https://arxiv.org/abs/2202.08095
  * [16] G. Kosec and B. Šarler, “Solution of a low prandtl number natural convection benchmark by a local meshless method,” _International Journal of Numerical Methods for Heat & Fluid Flow_, vol. 23, p. 22, 01 2013.