Data efficient surrogate modeling for engineering  design: Ensemble-free batch mode deep active learning for regression

add1,add2]Harsh Vardhancorrauth
[corrauth]Corresponding author
add2]Umesh Timalsina
add1,add2]Peter Volgyesi
add1,add2,add3]Janos Sztipanovits
[add1]Department of Computer Science, Vanderbilt University, 1400 18th Ave. S., Nashville, TN 37212-2846, United States
[add2]Institute for Software Integrated Systems, Vanderbilt University, 1025 16th Ave. S., Nashville, TN 37212-2328, United States
[add3]Department of Electrical and Computer Engineering, Vanderbilt University, 400 24th Ave. S., Nashville, TN 37212, United States

In a computer-aided engineering design optimization problem that involves notoriously complex and time-consuming simulator, the prevalent approach is to replace these simulations with a data-driven surrogate that approximates the simulator's behavior at a much cheaper cost. The main challenge in creating an inexpensive data-driven surrogate is the generation of a sheer number of data using these computationally expensive numerical simulations.
In such cases, Active Learning (AL) methods have been used that attempt to learn an input–output behavior while labeling the fewest samples possible. The current trend in AL for a regression problem is dominated by the Bayesian framework that needs training an ensemble of learning models that makes surrogate training computationally tedious if the underlying learning model is Deep Neural Networks (DNNs). However, DNNs have an excellent capability to learn highly nonlinear and complex relationships even for a very high dimensional problem.  To leverage the excellent learning capability of deep networks along with avoiding the computational complexity of the Bayesian paradigm, in this work we propose a simple and scalable approach for active learning that works in a student-teacher manner to train a surrogate model. By using this proposed approach, we are able to achieve the same level of surrogate accuracy as the other baselines like DBAL and Monte Carlo sampling with up to 40% fewer samples. We empirically evaluated this method on multiple use cases including three different engineering design domains:finite element analysis, computational fluid dynamics, and propeller design.

§ ACKNOWLEDGMENTS

This work is supported by DARPA’s Symbiotic Design for CPS project (FA8750-20-C-0537).

[Forrester et al., 2008]
authorA. Forrester, authorA. Sobester,
authorA. Keane, titleEngineering design via surrogate
modelling: a practical guide, publisherJohn Wiley & Sons,

[Köppen, 2000]
authorM. Köppen,
titleThe curse of dimensionality,
in: booktitle5th online world conference on soft computing
in industrial applications (WSC5), volume volume1,
year2000, pp. pages4–8.

[Rasmussen, 2003]
authorC. E. Rasmussen,
titleGaussian processes in machine learning,
in: booktitleSummer school on machine learning,
organizationSpringer, year2003, pp.

[Vardhan and Sztipanovits, 2022]
authorH. Vardhan, authorJ. Sztipanovits,
titleDeep learning-based fea surrogate for sub-sea
pressure vessel,
journalarXiv preprint arXiv:2206.03322

[Raina et al., 2019]
authorA. Raina, authorC. McComb,
authorJ. Cagan,
titleLearning to design from humans: Imitating human
designers through deep learning,
journalJournal of Mechanical Design volume141

[Vardhan et al., 2021]
authorH. Vardhan, authorP. Volgyesi,
authorJ. Sztipanovits,
titleMachine learning assisted propeller design,
in: booktitleProceedings of the ACM/IEEE 12th
International Conference on Cyber-Physical Systems, year2021,
pp. pages227–228.

[Ren et al., 2021]
authorP. Ren, authorY. Xiao,
authorX. Chang, authorP.-Y. Huang,
authorZ. Li, authorB. B. Gupta,
authorX. Chen, authorX. Wang,
titleA survey of deep active learning,
journalACM computing surveys (CSUR) volume54
(year2021) pages1–40.

[Shannon, 2001]
authorC. E. Shannon,
titleA mathematical theory of communication,
journalACM SIGMOBILE mobile computing and communications
review volume5 (year2001) pages3–55.

[Bishop, 1997]
authorC. M. Bishop,
titleBayesian neural networks,
journalJournal of the Brazilian Computer Society
volume4 (year1997) pages61–68.

[MacKay, 1992]
authorD. J. MacKay,
titleA practical bayesian framework for backpropagation
journalNeural computation volume4
(year1992) pages448–472.

[Konda and Tsitsiklis, 1999]
authorV. Konda, authorJ. Tsitsiklis,
titleActor-critic algorithms,
journalAdvances in neural information processing systems
volume12 (year1999).

[Sutton and Barto, 2018]
authorR. S. Sutton, authorA. G. Barto,
titleReinforcement learning: An introduction,
publisherMIT press, year2018.

[Simpson et al., 2001]
authorT. W. Simpson, authorJ. Poplinski,
authorP. N. Koch, authorJ. K. Allen,
titleMetamodels for computer-based engineering design:
survey and recommendations,
journalEngineering with computers volume17
(year2001) pages129–150.

[Pisner and Schnyer, 2020]
authorD. A. Pisner, authorD. M. Schnyer,
titleSupport vector machine,
in: booktitleMachine learning,
publisherElsevier, year2020, pp.

[Breiman, 2001]
authorL. Breiman,
titleRandom forests,
journalMachine learning volume45
(year2001) pages5–32.

[Palar et al., 2019]
authorP. S. Palar, authorR. P. Liem,
authorL. R. Zuhal, authorK. Shimoyama,
titleOn the use of surrogate models in engineering design
optimization and exploration: The key issues,
in: booktitleProceedings of the Genetic and Evolutionary
Computation Conference Companion, year2019, pp.

[LeCun et al., 2015]
authorY. LeCun, authorY. Bengio,
authorG. Hinton,
titleDeep learning,
journalnature volume521
(year2015) pages436–444.

[Vesting and Bensow, 2014]
authorF. Vesting, authorR. E. Bensow,
titleOn surrogate methods in propeller optimisation,
journalOcean engineering volume88
(year2014) pages214–227.

[Liang et al., 2018]
authorL. Liang, authorM. Liu,
authorC. Martin, authorW. Sun,
titleA deep learning approach to estimate stress
distribution: a fast and accurate surrogate of finite-element analysis,
journalJournal of The Royal Society Interface
volume15 (year2018) pages20170844.

[del Rio-Chanona et al., 2019]
authorE. A. del Rio-Chanona, authorJ. L. Wagner,
authorH. Ali, authorF. Fiorelli,
authorD. Zhang, authorK. Hellgardt,
titleDeep learning-based surrogate modeling and
optimization for microalgal biofuel production and photobioreactor design,
journalAIChE Journal volume65
(year2019) pages915–923.

[Settles, 2009]
authorB. Settles,
titleActive learning literature survey

[Fedorov, 2013]
authorV. V. Fedorov, titleTheory of optimal
experiments, publisherElsevier, year2013.

[MacKay, 1992]
authorD. J. MacKay,
titleInformation-based objective functions for active data
journalNeural computation volume4
(year1992) pages590–604.

[Matsushita et al., 2018]
authorK. Matsushita, authorK. Matsushita,
authorHasebe, titleDeep active learning,
publisherSpringer, year2018.

[Sener and Savarese, 2017]
authorO. Sener, authorS. Savarese,
titleActive learning for convolutional neural networks: A
core-set approach,
journalarXiv preprint arXiv:1708.00489

[Gal and Ghahramani, 2015]
authorY. Gal, authorZ. Ghahramani,
titleBayesian convolutional neural networks with bernoulli
approximate variational inference,
journalarXiv preprint arXiv:1506.02158

[Kirsch et al., 2019]
authorA. Kirsch, authorJ. Van Amersfoort,
authorY. Gal,
titleBatchbald: Efficient and diverse batch acquisition
for deep bayesian active learning,
journalAdvances in neural information processing systems
volume32 (year2019).

[Ash et al., 2019]
authorJ. T. Ash, authorC. Zhang,
authorA. Krishnamurthy, authorJ. Langford,
authorA. Agarwal,
titleDeep batch active learning by diverse, uncertain
gradient lower bounds,
journalarXiv preprint arXiv:1906.03671

[Tsymbalov et al., 2018]
authorE. Tsymbalov, authorM. Panov,
authorA. Shapeev,
titleDropout-based active learning for regression,
in: booktitleInternational conference on analysis of
images, social networks and texts, organizationSpringer,
year2018, pp. pages247–258.

[Uesato et al., 2018]
authorJ. Uesato, authorA. Kumar,
authorC. Szepesvari, authorT. Erez,
authorA. Ruderman, authorK. Anderson,
authorN. Heess, authorP. Kohli, et al.,
titleRigorous agent evaluation: An adversarial approach to
uncover catastrophic failures,
journalarXiv preprint arXiv:1812.01647

[Vardhan and Sztipanovits, 2021]
authorH. Vardhan, authorJ. Sztipanovits,
titleRare event failure test case generation in
in: booktitle2021 6th International Conference on Machine
Learning Technologies, year2021, pp. pages34–40.

[Wei et al., 2015]
authorK. Wei, authorR. Iyer,
authorJ. Bilmes,
titleSubmodularity in data subset selection and active
in: booktitleInternational conference on machine
learning, organizationPMLR, year2015, pp.

[Wolf, 2011]
authorG. W. Wolf, titleFacility location: concepts,
models, algorithms and case studies. series: Contributions to management
science: edited by zanjirani farahani, reza and hekmatfar, masoud,
heidelberg, germany, physica-verlag, 2009, 549 pp.,€ 171.15,
$219.00, 144.00, isbn 978-3-7908-2150-5 (hardprint),
978-3-7908-2151-2 (electronic), year2011.

[Zhdanov, 2019]
authorF. Zhdanov,
titleDiverse mini-batch active learning,
journalarXiv preprint arXiv:1901.05954

[Vardhan and Sztipanovits, 2022]
authorH. Vardhan, authorJ. Sztipanovits,
titleDeepal for regression using epsilon-weighted hybrid
query strategy,
journalarXiv preprint arXiv:2206.13298

[Auer et al., 2002]
authorP. Auer, authorN. Cesa-Bianchi,
authorP. Fischer,
titleFinite-time analysis of the multiarmed bandit
journalMachine learning volume47
(year2002) pages235–256.

[Riegel et al., 2016]
authorJ. Riegel, authorW. Mayer,
authorY. van Havre, titleFreecad,

[Epps et al., 2009]
authorB. Epps, authorJ. Chalfant,
authorR. Kimball, authorA. Techet,
authorK. Flood, authorC. Chryssostomidis,
titleOpenprop: An open-source parametric design and
analysis tool for propellers,
in: booktitleProceedings of the 2009 grand challenges in
modeling & simulation conference, year2009, pp.

[Jasak et al., 2007]
authorH. Jasak, authorA. Jemcov,
authorZ. Tukovic, et al.,
titleOpenfoam: A c++ library for complex physics
in: booktitleInternational workshop on coupled methods in
numerical dynamics, volume volume1000,
organizationIUC Dubrovnik Croatia, year2007, pp.

[Hyder and Asif, 2008]
authorM. J. Hyder, authorM. Asif,
titleOptimization of location and size of opening in a
pressure vessel cylinder using ansys,
journalEngineering Failure Analysis volume15
(year2008) pages1–19.

[Shigley, 1972]
authorJ. E. Shigley, titleMechanical engineering
design, publisherMcGraw-Hill Companies, year1972.

[Coney, 1989]
authorW. B. Coney, titleA method for the design of a
class of optimum marine propulsors, Ph.D. thesis, Massachusetts Institute of
Technology, year1989.

[Myring, 1976]
authorD. Myring,
titleA theoretical study of body drag in subcritical
axisymmetric flow,
journalAeronautical quarterly volume27
(year1976) pages186–194.

[Sagaut et al., 2013]
authorP. Sagaut, authorM. Terracol,
authorS. Deck, titleMultiscale and multiresolution
approaches in turbulence-LES, DES and Hybrid RANS/LES Methods: Applications
and Guidelines, publisherWorld Scientific, year2013.

[Menter, 1992]
authorF. R. Menter, titleImproved two-equation k-omega
turbulence models for aerodynamic flows, typeTechnical Report,


aainstitutetext: Institut für Kernphysik, Westfälische Wilhelms-Universität
Münster,
Wilhelm-Klemm-Str. 9, Münster, Germanybbinstitutetext: Research Division and
ExtreMe Matter Institute EMMI, GSI Helmholtzzentrum für Schwerionenforschung
GmbH, Planckstr. 1, Darmstadt, Germanyccinstitutetext: Physikalisches
Institut, Ruprecht-Karls-Universität Heidelberg,
Im Neuenheimer Feld 226, Heidelberg, Germanyddinstitutetext: Lund University,
Professorsgatan 1, Lund, Sweden

# Transverse dynamics of charmed hadrons in ultra-relativistic nuclear
collisions

Anton Andronic b,c Peter Braun-Munzinger c Hjalmar Brunßen c Jana Crkovská c
Johanna Stachel d Vytautas Vislavicius c Martin Völkl  andronic@uni-
muenster.de<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Transverse momentum $p_{\text{T}}$ spectra and anisotropic flow distributions
are studied for charmonia and charmed hadrons produced in Pb–Pb collisions and
measured with the ALICE detector at the CERN Large Hadron Collider (LHC). The
investigations are performed within the framework of the Statistical
Hadronization Model with the transverse dynamics evaluated using predictions
from relativistic viscous hydrodynamics as implemented in the computer codes
MUSIC and Fluid$u$M. With this essentially parameter-free approach good
agreement is obtained for $p_{\text{T}}$ spectra in the range $p_{\text{T}}$
$<10$ GeV/c. The observed wide distribution in $p_{\text{T}}$ of anisotropic
flow coefficients v2 and v3 for charmonia is also well reproduced, while their
magnitude is generally somewhat over predicted. This finding may be connected
to a difference in spatial distribution between light and charmed hadrons due
to a different diffusion of light and heavy quarks in the hot fireball.

##  1 Introduction

The production and hadronization of charm quarks is of central interest in
current quark-gluon plasma (QGP) research ALICE:2022wpn . In particular,
hadrons with one or more charm quarks such as charmonia or (multi-)charm
hadrons play an increasing role in our quest to quantitatively understand the
expansion and hadronization of the QGP fireballs formed in ultra-relativistic
nuclear collisions with heavy beams. In addition, the production yields of
such hadrons contain unique information about the degree of deconfinement of
charm quarks reached in such collisions Andronic:2003zv ; Andronic:2017pug ;
Andronic:2021erx ; Gross:2022hyw .

The production probabilities of charmed hadrons in central Pb–Pb collisions
are very well described within the framework of the statistical hadronization
model (SHMc), see Braun-Munzinger:2000csl ; Andronic:2017pug  and the
description in the following section 2. Within this framework alone there is
no information about the transverse dynamics and, in particular, the
transverse momentum spectra of the produced charm hadrons. Recently, we have
shown Andronic:2019wva ; Andronic:2021erx  that transverse momentum spectra
can be well described if one couples hydrodynamical Song:2007ux  or hydro-
inspired models such as developed within the ’blast wave’ approach
Siemens:1978pb ; Schnedermann:1993ws  to the SHMc. Using these models it is
assumed that the particles under consideration, i.e. in our case the charm
quarks, participate in the collective expansion and flow developed in the QGP
fireball as a consequence of the pressure build-up and and pressure gradients
in the hot medium. This success implies that a much more direct route to
transverse dynamics in the SHMc is to implement the charmed hadrons, e.g.
charmonia, directly into one of the now available relativistic hydrodynamics
simulation frameworks. In particular, we focus on hydrodynamic modeling of
open charm hadrons and charmonia in the framework of MUSIC Schenke:2010rr  and
of Fluid$u$M Floerchinger:2018pje . How this is done in detail is described in
section 3 below. A key ingredient of the current approach is that the
integrated yield for each charm hadron is taken directly from the SHMc
calculation, and that, in order to describe the transverse spectral shape, the
hydrodynamic evolution is stopped at the QCD phase boundary. Here our explicit
assumption is that, for charmed hadrons, not only the yields but also the
spectra freeze out at the phase boundary. For light valence flavor hadrons
there is, on the other hand, indication that, while hadron yields are frozen,
the expansion continues for a couple of fm/c in the now non-equilibrium
hadronic phase and that transverse momentum spectra are determined at kinetic
freeze out at a somewhat lower effective temperature. At LHC energies this
implies that for charmed hadrons the evolution is stopped when the QGP
temperature reaches T${}_{pc}=156.5$ MeV. All charmed hadrons are then
produced at the same thermal conditions as hadrons with light (u,d,s) valence
quarks Andronic:2005yp ; Andronic:2017pug , with the one major difference of
the introduction of a charm quark fugacity $g_{c}$ Braun-Munzinger:2000csl ;
Andronic:2021erx . The quantity $g_{c}$ is introduced since charm quarks are
dominantly produced in initial hard collisions, i.e. never reach chemical
equilibrium but thermalize in terms of their momentum distributions
subsequently in the hot QGP medium. Thus, $g_{c}$ is not a free parameter but
determined from the measured charm quark content of the fireball. Of course,
as in all statistical hadronization models, the approach needs information
about the full mass spectrum of charmed hadrons and their strong decays, see
Andronic:2021erx .

The main idea pursued in the present manuscript is to study the implications
of the detailed hydrodynamic evolution on the transverse dynamics. In
particular, we compare the predictions of the above cited two, in their
technical implementation very different, hydrodynamic model codes with the
measured transverse momentum spectra of charmonia and charm hadrons. In
addition, we investigate how well the recently measured azimuthal anisotropies
of charmonia ALICE:2013xna ; ALICE:2020pvw  and charm hadrons ALICE:2013olq
are reproduced. This is a particularly interesting question as first attempts
in this direction missed the experimental data in particular at higher values
of $p_{\text{T}}$ $>$ 4 GeV/c Zhou:2014kka ; Du:2015wha . In a recent
publication He:2021zej  good agreement was obtained by employing Langevin
simulations to more accurately account for spatial anisotropies and
introducing space-momentum correlations of the charm quarks during the
hydrodynamic expansion of the hot fireball. Below we show that both transverse
momentum spectra and anisotropic flow distributions are in rather good
agreement with SHMc predictions under the assumption that charm quarks fully
participate in the hydrodynamic expansion. This is realized by directly
inserting the charmonia when the hydrodynamic evolution is stopped at the
pseudo critical temperature of the chiral phase transition such that charmonia
inherit the full flow of the medium at this instance. The direct use of the
hydrodynamic models rather than the simplification of a blast wave model (as
in Andronic:2021erx ) retains complexities of the medium evolution. For
example, while in blast wave models freeze out is realized as a simple
integration along a freeze out contour in a plane of proper time vs radius, in
a hydrodynamic model individual fluid cells are followed in time and freeze
out occurs when a given energy density or temperature is reached. This implies
complex geometries of the freeze out in space and time, with velocity vectors
of a given fluid cell and freeze out surface normals pointing in all possible
directions. Inspecting the radial profile of velocities of fluid cells, one
finds indeed a strong correlation (as employed in the blast wave model) but,
in particular for larger radii of 5-10 fm, a very wide distribution of
velocities Andronic:2021erx . We also note that freeze out at early times
actually gives a substantial contribution to particle spectra from fluid cells
at the surface of the fireball and with a large velocities. In addition, the
use of anisotropic modeling (already in 2+1 dimensional hydrodynamics) allows
for calculation of the flow coefficients.

In section 6 we show how a blast wave parametrization can be optimally matched
to reproduce the full hydrodynamics calculation of the ${\rm J}/\psi$
spectrum. Finally in section 7 we discuss and estimate the consequence of the
charm spatial distribution possibly being more compact than the energy
distribution at freeze out.

##  2 Brief summary of the statistical hadronization model for heavy quarks,
SHMc

Here we summarize the main ideas behind the SHMc. For much more detail on the
original development see Andronic:2021erx ; Andronic:2017pug ; Andronic:2003zv
; Braun-Munzinger:2000csl . Our main emphasis will be on the production of
charmonia and in particular on the connection between SHMc for yields, and
hydrodynamic models for the description of transverse dynamics. The key idea
is based on the recognition that, contrary to what happens in the (u,d,s)
sector, the heavy (mass $\sim$ 1.2 GeV) charm quarks are not thermally
produced. Rather, production takes place in initial hard collisions. The
produced charm quarks then thermalize in the hot fireball, but the total
number of charm quarks is conserved during the evolution of the fireball
Andronic:2006ky  since charm quark annihilation is very small. In essence,
this implies that charm quarks can be treated like impurities. Their thermal
description then requires the introduction of a charm fugacity $g_{c}$ Braun-
Munzinger:2000csl ; Andronic:2021erx . The value of $g_{c}$ is not a free
parameter but experimentally determined by measurement of the total charm
cross section, relying on the precisely measured cross section for D0
production in Pb–Pb collisions ALICE:2021rxa  and fragmentation as in the
SHMc, see the section on quark matter in Gross:2022hyw . For 0-10% central
Pb–Pb collisions at 5.02 TeV, the number of $\mathrm{c\bar{c}}$ pairs per unit
of rapidity is 13.8$\pm$2.1 (16.3$\pm$2.4 with an enhanced charm-baryon
spectrum, see below) resulting in $g_{c}=31.5\pm 4.7$ and in a D0 to ${\rm
J}/\psi$ ratio of 53, in very good agreement with recent precision data
ALICE:2023gco . Note that, with this approach, unlike in the case in which one
starts with the charm production cross section in pp collisions, the
uncertainties related to gluon shadowing in the Pb nucleus are avoided,
leading to reduced uncertainties in $g_{c}$ (currently $\pm$15%).

We note here that the charm balance equation should contain canonical
corrections for more peripheral collisions or for lighter collision systems
(or lower collision energies), i.e., whenever the number of charm pairs is not
large compared to 1 Gorenstein:2000ck ; BraunMunzinger:2003zd ;
Andronic:2021erx . For central Pb–Pb collisions at 5.02 TeV the canonical
correction factors are in fact very close to 1 Andronic:2021erx . For more
peripheral collisions or smaller systems the canonical corrections are
important. We follow everywhere the canonical treatment described in detail in
Andronic:2021erx .

The large value of $g_{c}\approx 30$ for charm production in central Pb–Pb
collisions at mid rapidity implies very large enhancements for charmed hadrons
compared to what is obtained in the purely thermal case. In the absence of
canonical corrections the enhancement factor is about 900 for charmonia and
doubly charmed, and $2.6\cdot 10^{4}$ for triply charmed hadrons.

The statistical hadronization approach is applied to the core of the nuclear
overlap region. As the density drop off in nuclei is gradual, we define a
corona region. We assume this to start where the density of one of the two
colliding nuclei has dropped to 10 % of the central density. The rationale is
that in the overlap of this tail of the distribution with the core or tail of
the density distribution of the other nucleus, the average number of
collisions is less than 1. The two regions are derived from a Glauber model
calculation as a function of centrality. As an example, for the 0-10 %
centrality bin, the total overlap region contains 361.3 participants, of which
we define 340.6 as the core. In this core there are 1645 binary collisions.

##  3 Extracting momentum information from hydrodynamic calculations

As described in the previous section, the core part of the fireball is assumed
to hadronize statistically while preserving the fixed number of charm and anti
charm quarks produced in hard initial collisions. This implies that the
relative distribution of the charm quarks among hadrons is driven by the
statistical weights at a given temperature. We use the temperature determined
for chemical freeze out obtained when fitting the yields of hadrons with u,d,s
valence quarks, Tchem = 156.5 MeV Andronic:2017pug , which happens to coincide
with the pseudo critical temperature for the chiral phase transition Tpc.
While the statistical hadronization model as such makes no prediction about
the momentum distributions of the produced hadrons, the underlying assumption
is nevertheless that the hadronizing charm quarks are in (or close to) local
thermal equilibrium kinetically. Measurements of the various D-meson spectra
and elliptic flow combined with model comparisons confirm this assumption
ALICE:2021rxa . Based on this statistical approach of charmed hadron formation
a consistent prediction of the transverse dynamics is obtained by using the
hydrodynamic modelling optimized for observables of light valence quark
hadrons, albeit at Tpc. Local thermal equilibrium gives the momentum dependent
yields in the rest frame of the fluid, i.e. in each fluid cell. Employing the
boosts of all fluid cells and their geometric arrangement, transverse momentum
spectra including azimuthal asymmetries can be obtained as described in this
section.

In hydrodynamic simulations, the transition from the hydrodynamic phase to a
hadron gas with little interaction is usually not modeled as part of the
system evolution. Instead, the system evolves in time fully hydrodynamically.
The end of the evolution is reached by defining a transition hypersurface,
delineating the points in space and time where a certain critical temperature
or energy density are reached. At this hypersurface particles are produced
according to the Cooper-Frye formalism Cooper:1974mv :

$E\frac{{\rm d}N}{{\rm d}^{3}p}=g_{i}\int_{\Sigma}f(u^{\mu}p_{\mu})p^{\mu}{\rm
d}^{3}\Sigma_{\mu}~{}.$ (1)

This formalism assumes a thermalized momentum distribution $f$ at the freeze
out hypersurface and then counts the particles crossing the boundary. We
assume that the charmed hadrons are produced at the same rate in all unit
volumes at the transition temperature. This means that information about the
shape of the overall momentum distribution can be obtained by calculating the
freeze out momentum distribution from the hydrodynamic simulations according
to 1. The integral of the distribution is then scaled according to the SHMc as
described above. Compared to a description using a blast wave approach this
gives a much more realistic description of the freeze out process. Here, the
freeze out has significant contributions from hypersurfaces that are volume
like (with time like surface normal vectors ${\rm d}^{3}\Sigma_{\mu}$) and
ones that are surface like (with space like normal vectors). In the blast wave
approach only volume like contributions are considered. Also, the full
complexity of the radius-time-velocity association is reflected in the
spectrum, while in the blast wave approach a simplified tight connection
between these variables is employed. Further and new information on the
physics behind the Cooper-Frye formalism can be found in Kirchner:2023fsj .

For our purposes we use two hydrodynamic models: MUSIC Schenke:2010rr  and
Fluid$u$M Floerchinger:2018pje . MUSIC was used with the settings of
Schenke:2020mbo  and using a 2+1D, boost-invariant setup for the measurements
at mid rapidity. This includes an initial state based on the IP-Glasma model
Schenke:2012hg ; Schenke:2012wb . In Fluid$u$M, the fluid dynamic fields are
decomposed into an azimuthally and Bjorken boost symmetric background
contribution and perturbations. This allows a numerically very efficient
calculation of the hydrodynamic evolution. For the calculation of transverse
momentum spectra, terms of linear order in perturbation do not contribute so
that in this case the azimuthally symmetric background evolution is
sufficient. Fluid$u$M was used with the parameters as in Devetak:2019lsk  and
a normalization of the energy density in the initial state model TRENTo
Moreland_2015  that reproduced the integrated pion, kaon and proton yields
measured in central Pb–Pb collisions. In a very recent paper, the diffusion of
charm quarks and its effect on transverse momentum spectra of charm hadrons is
discussed in the context of the Fluid$u$M model Capellino:2023cxe . The
settings of both models are summarized in Table 1. In both cases, a constant
ratio of shear viscosity to entropy density $\eta/s$ is assumed. The ratio of
bulk viscosity to entropy density $\zeta/s$ is assumed to peak somewhat above
the pseudo critical temperature and the values and positions of the maximum
are given in the table. The precise functions can be found in Schenke:2020mbo
and Devetak:2019lsk , respectively.

Table 1: MUSIC and Fluid$u$M parameters model | initial conditions | $\tau_{0}$ | $\eta/s$ | $\left(\zeta/s\right)_{\text{max}}$ | $T_{\mathrm{peak}}$ | $T_{\text{fo}}$
---|---|---|---|---|---|---
Fluid$u$M | TRENTo | $0.18\,\mathrm{fm/c}$ | 0.16 | 0.06 | $175\,\mathrm{MeV}$ | $156\,\mathrm{MeV}$
MUSIC | IP-Glasma | $0.4\,\mathrm{fm/c}$ | 0.12 | 0.13 | $160\,\mathrm{MeV}$ | $156\,\mathrm{MeV}$

##  4 Calculation of ${\rm J}/\psi$ spectrum with MUSIC and Fluid$u$M

Figure 1 shows the distribution of the freeze out hypersurfaces in time and
radius. In Fluid$u$M, the calculation is started with a rotationally symmetric
solution as the basic output around which perturbations are calculated. This
means that the resulting freeze out hypersurface is azimuthally isotropic,
which gives a line in the $\tau$-$R$ space.

The $p_{\text{T}}$ distribution of thermal ${\rm J}/\psi$ mesons resulting
from this freeze out represents the contribution of the core in the SHMc
model. To this, a pp-like contribution of the corona is added (see above). The
third component to the ${\rm J}/\psi$ spectrum is obtained from feed down from
weak decays of beauty hadrons. This dominantly contributes at high momenta.
Beauty quarks interact with the medium, but are likely not fully thermalized
Andronic:2022ucg . Thus, the transverse momentum distribution of the beauty
hadrons is expected to lie somewhere between the extremes of a pp-like
distribution and a fully thermalized distribution with strong effects of
collective flow as represented by hydrodynamic simulations. We chose the
arithmetic mean of the two cases. For this, we used a pp-like $p_{\text{T}}$
distribution obtained from FONLL Cacciari:1998it  and a thermalized
distribution represented by blast wave parameters described below. We included
B mesons by using PYTHIA8 Bierlich:2022pfr  to generate their weak decays to
${\rm J}/\psi$. For the sum of ${\rm B}^{+}$, ${\rm B}^{-}$, ${\rm B}^{0}$,
and $\bar{\rm B}^{0}$ we used a total yield after all strong decays of ${\rm
d}N/{\rm d}y=0.86$ and 0.116 for 0-10% and 30-50% centrality ranges as given
by SHM calculations for the bottom sector at mid rapidity Andronic:2022ucg .

Figure 1: Left: $R$ and $\tau$ distribution of the freeze out hypersurfaces
for MUSIC (histogram) and Fluid$u$M (black line). The MUSIC hypersurface
elements are weighted by their volume (but not the energy streaming outwards
through this element). Right: Contributions to the ${\rm J}/\psi$
$p_{\text{T}}$-distribution in the SHMc model. The extraction of the blast
wave parameters is discussed in section 6.

Figure 2: Transverse momentum distributions of ${\rm J}/\psi$ production in
Pb–Pb collisions at mid rapidity. The result from the SHMc model using the two
hydrodynamical generators is compared to the measurement by the ALICE
Collaboration ALICE:2023gco . The bottom plot show ratios of model results to
the measurements. There the red data points are displayed at unity to indicate
the experimental uncertainties.

Figure 2 shows the resulting $p_{\text{T}}$ distributions for ${\rm J}/\psi$
mesons produced in central and semi-central collisions at mid rapidity
together with the ALICE measurement at 5 $\rm TeV$. The centrality selection
for both models is based on final state particle multiplicities as used in the
ALICE experiment. The vertical size of the boxes shown reflects the
uncertainty in the charm cross section, i.e. of $g_{c}$, the uncertainty of
the corona contribution due to the experimental uncertainty of the spectrum
measured for pp collisions, as well as an uncertainty from the modeling of the
feed down from beauty hadrons. The overall agreement is rather good
considering that there is no free parameter adjusted to data. The general
shape of the distributions is very similar. For the most central bin, the SHMc
model gives a distribution that is somewhat harder than the experimental data.
The difference is in the intermediate $p_{\text{T}}$ region dominated by the
thermalized core contribution from the hydrodynamic simulations, with very
similar results from MUSIC and Fluid$u$M. In the semi-central case, where the
transverse flow effect is reduced, the agreement with the data is rather
close.

For light hadrons ALICE:2019hno , central as compared to semi-central events
result in a larger average $p_{\text{T}}$ for each particle species. This is
well reproduced by an increased mean transverse velocity
$\langle\beta_{\perp}\rangle$ obtained from hydrodynamic models and also from
simple blast wave fits to the data. We obtain a corresponding shift in the
core contribution using both hydrodynamic generators as seen in Fig. 2. In the
experimental data, however, the maximum of the $p_{\text{T}}$ distribution is,
if anything, slightly lower for the more central events. This may point to
some additional effect beyond the assumption that the charm quarks are carried
along the hydrodynamical flow of the medium.

In Figure 3 we show the corresponding model prediction for the $\psi(2S)$
spectrum. The model uncertainties were obtained as for the ${\rm J}/\psi$ case
above. In comparison we show experimental data from ALICE for Pb–Pb collisions
in the 0-90 % centrality range ALICE:2022jeh . While the overall spectral
shape is in reasonable agreement, the data exceed the model predictions at low
$p_{\text{T}}$ by 1-2 standard deviations. High precision data expected from
ALICE Run 3 for central collisions are needed to establish whether there are
any significant differences between model predictions and data.

Comparing to our first modelling Andronic:2019wva  of charmonium spectra using
a blast wave model, we note that the current spectra based on full
hydrodynamics are somewhat harder. The most probable $p_{\text{T}}$ increases
from about 1.5 GeV/$c$ in the blast wave calculations to 2.5 GeV/$c$ with full
hydrodynamic treatment, a trend which was in fact expected. The blast wave
model parameters were chosen to match the MUSIC distribution of the freeze out
hypersurface with a polynomial of the velocity $\beta$ as function of radius.
While the MUSIC distribution could be well matched with an exponent of n =
0.85 and a value of $\beta_{max}$ = 0.62, it was apparent that the
distribution was very broad for fluid elements at radii between 5 and 10 fm
reaching from velocities of about 0.1 to 0.85 and therefore some arbitrariness
in how to define the best $\beta_{max}$ (see section 6 below).

Figure 3: Calculation of the $\psi(2S)$ $p_{\text{T}}$ spectrum in central
Pb–Pb collisions with the SHMc and MUSIC. In comparison are shown experimental
data from ALICE ALICE:2022jeh . The error band is obtained as for the ${\rm
J}/\psi$ spectrum.

##  5 Modeling of ${\rm J}/\psi$ anisotropic flow coefficients

The flow coefficients quantify the azimuthal anisotropy of particle
distributions, driven in hydrodynamic models by the azimuthal anisotropy and
resulting pressure gradients in the initial state and reduced by dissipation
during the expansion phase. In the MUSIC calculations we have full access to
the IP-Glasma initial state. We use the definition of the minor axis of the
initial state ellipse and triangularity

$\psi_{n}=\frac{{\rm
atan2}(r^{2}\langle\sin(n\phi)\rangle,r^{2}\langle\cos(n\phi)\rangle)+\pi}{n}~{},$
(2)

as introduced in Alver:2010gr . Here, the average was taken over the
distribution of the energy density in the initial state rather than that of
the participants used in the original approach. The flow coefficients $v_{2}$
and $v_{3}$ were then calculated as

$v_{n}=\langle\cos(n(\phi-\psi_{n}))\rangle~{},$ (3)

with the averaging done over the expectation values of the number of ${\rm
J}/\psi$ in each $\phi$ bin. The resulting flow corresponds to fully
thermalized charm quarks flowing with the medium and producing ${\rm J}/\psi$
mesons during hadronization. This was taken as the flow contribution of the
core part of the model. The contribution of the corona was assumed to be
independent of the angle $\phi$, consistent with initial production without
further interaction. The feed down from beauty hadrons was modeled as having a
$v_{2}$ and $v_{3}$ constant in $p_{\text{T}}$ over the range shown, with
values of $v_{2}=0.01\pm 0.01$ and $v_{3}=0.0\pm 0.02$ for the 0-10% and
$v_{2}=0.05\pm 0.02$ and $v_{3}=0.0\pm 0.02$ for 30-50% centrality bins. These
values were estimated from the measurements in ATLAS:2020yxw .

Figure 4: The $v_{2}$ coefficients of the SHMc model prediction combined with
the 2+1D MUSIC calculations. For the error band see text. The results are
compared with data from the ALICE Collaboration at central and forward
rapidity for central (left) and semi-central (right) Pb–Pb collisions
ALICE:2020pvw .

Figure 4 shows the resulting $p_{\text{T}}$-dependent elliptic flow
coefficients for ${\rm J}/\psi$ in central and semi-central events. The error
band contains the uncertainty of the overall charm cross section, the
experimental uncertainty of the pp spectrum and the uncertainty attributed to
the beauty feed down contribution. The elliptic flow at large $p_{\text{T}}$
is suppressed by the increasing contribution from the corona. This leads to a
peak of the elliptic flow around 5 $\mathrm{GeV}/c$. For very high transverse
momenta, the $v_{2}$ of the core contribution becomes negative. As this
happens at transverse momenta where the core contribution is already very
small, this has little effect on the total $v_{2}$ but is nevertheless
interesting. We find that the negative values arise from regions close to the
edge of the medium moving outward very rapidly and freezing out at early
times. The contribution is stronger from fluid cells close to the tip of the
ellipse. It seems plausible that this could be an effect of the curvature of
the surface of the ellipse. The interface with neighboring fluid cells is
smaller at the tip, leading to less slowing down due to dissipation of a given
fluid cell in that region and hence larger momenta in this direction resulting
in a negative $v_{2}$. Further investigations are necessary to demonstrate
this explicitly.

The resulting $v_{2}$ exhibit a substantial flow contribution of thermalized
${\rm J}/\psi$ out to about 8-10 $\mathrm{GeV}/c$ . This suggests that the
flow observed in data even at these large transverse momenta may not be due to
the effect of the path length difference, as commonly assumed for the high
$p_{\text{T}}$ region, but still is of hydrodynamic origin coming from the
anisotropy of the fluid cells and their motion.

Comparison to ALICE data at central and forward rapidities shows good
agreement for central collisions. For semi-central events the range in
$p_{\text{T}}$ of sizable $v_{2}$ values is reasonably well reproduced, while
the model produces about twice the maximum flow amplitude.

Figure 5: The $v_{3}$ coefficients of the SHMc model combined with the 2+1D
MUSIC calculations. The model results are compared with data from the ALICE
Collaboration at forward rapidity for central (left) and semi-central (right)
Pb–Pb collisions ALICE:2020pvw .

The $v_{3}$ coefficients exhibit a similar behaviour for both central and more
peripheral events. Similar to the $v_{2}$, the contribution from the core
gives a substantial contribution up to $8-9~{}\mathrm{GeV}/c$, which is in
line with measurements. The triangular flow coefficients from our model peak
at values around 0.08 for transverse momenta of 5-6 $\mathrm{GeV}/c$. These
peak values are somewhat larger than the experimental observations. For
comparison, the $v_{2}$ flow coefficients measured by ALICE in the forward
direction are larger by roughly a factor of two.

##  6 Matching hydrodynamics to blast wave parametrization and resulting open
charm spectra

To calculate the momentum distributions of open charm hadrons in the SHMc, the
feed down of a large number of charm hadron states needs to be taken into
account. In the SHMc 55 charmed mesons and 74 charmed (anti) baryons with
their decay branching fractions are included. They are not usually included in
the hydrodynamics codes, where we included the charmonia considered in the
previous section explicitly. The inclusion of as complete a mass spectrum as
possible is very important: Strong decays quadruple the D meson yields, for
the $\Lambda_{\rm c}$ the increase is by a factor 5 Andronic:2021erx . To
incorporate this the FastReso FastReso  framework provides the necessary
capability including the decay kinematics of the open charm hadrons. The input
to the calculation is provided by a blast wave model Schnedermann:1993ws ;
Florkowski:2010zz ; Andronic:2019wva .

This model uses a power law dependence of velocity of radius
$\beta(r)=\beta_{\rm max}(r/r_{\rm max})^{n}$ with a maximum transverse medium
velocity $\beta_{\rm max}$ and an exponent $n$ of the radial dependence as
input parameters. While we already employed this ansatz in our previous
publications Andronic:2019wva ; Andronic:2021erx , the matching of the blast
wave parameters to the hydrodynamic calculation was done differently here.
Since we have now at our disposal the ${\rm J}/\psi$ spectrum from full
hydrodynamics as shown in the previous section, we can find an optimum
matching of the parametrization to full hydrodynamics.

In general, the production of particles from a part of the freeze out
hypersurface depends on the mass of the particles, as given in equation 1. The
freeze out hypersurface is in fact quite different in the full hydrodynamic
model when compared to schematic hypersurfaces used in the blast wave
approach. In particular it also has surface like contributions with space like
normal vectors. Thus, the blast wave parameters which best reproduce the
hydrodynamic model depend on the particle mass. However, in the limit of large
masses, this matching becomes easier again. Dividing eq. 1 by the particle
energy, we get:

$\frac{{\rm d}N}{{\rm
d}^{3}p}=g_{i}\int_{\Sigma}f(u^{\mu}p_{\mu})\frac{1}{\gamma}u^{\mu}{\rm
d}^{3}\Sigma_{\mu}~{}.$ (4)

In the limit of large masses, the particle velocity is approximately equal to
the medium velocity as thermal velocity fluctuations are suppressed.
Therefore, we can use $u$ and $\gamma$ corresponding to the medium. Then,
$\frac{1}{\gamma}u^{\mu}{\rm d}^{3}\Sigma_{\mu}$ is independent of the
particle species. If the freeze out hypersurface is taken at constant time
($\Sigma_{\mu}=(\Sigma_{0},\vec{0})$), $\frac{1}{\gamma}u^{\mu}{\rm
d}^{3}\Sigma_{\mu}$ is simply the volume of the fluid cell. This volume can be
used as a weight to take into account how fast the medium flows through the
hypersurface element and thereby how much it contributes to the spectrum of
the (large mass) particle. The velocities of each hypersurface element are
multiplied with it’s volume as a weight to find the best match of a blast wave
parametrization to the full hydrodynamic result for charmed hadrons.

Figure 6: Comparison of velocity distribution of the freeze out hypersurface
elements obtained from MUSIC (green line) to a matched blast wave distribution
with a smoothed edge (short dashed blue line) and a velocity matched blast
wave distribution (red dashed line).

In contrast to the previous method of extraction, this ignores the spatial
position of the individual fluid cells and only uses their velocity
distribution. The effort here is thus to optimally reproduce the resulting
charmed hadron (${\rm J}/\psi$) momentum distribution of the hydrodynamic code
with the blast wave model rather than reproducing the $R-\beta$ distribution.
In general, there is no need to equate the radius parameter $r$ in the blast
wave model with the radius in the hydrodynamic simulation. For a radius-
independent freeze out time, the velocity distribution of the blast wave model
has the form of a power law ${\rm d}N/{\rm d}\beta\sim\beta^{\alpha}$, with an
abrupt drop to $0$ at $\beta_{\rm max}$. Figure 6 shows the volume weighted
distribution of velocities of the freeze out hypersurface elements from MUSIC.
Obviously there is no sharp drop but rather a smooth fall off at large
$\beta$. This is due to the generally broad distribution of hypersurface
elements as shown in Fig. 1 and also event-by-event fluctuations. The matching
to the MUSIC distribution is done in two steps. First the distribution from
MUSIC is fitted with a power law multiplied by an error function at large
$\beta$. As shown in Fig. 6 this distribution (blue short dashed) curve
represents a very close match to the MUSIC output. This fit yields the power
law parameter of the blast wave parametrization, $n=2/(\alpha+1)$. Next, using
the power law exponent $n$, the $\beta_{\rm max}$ of the blast wave
parametrization is fixed such that its $\langle\beta\gamma\rangle$ is equal to
the one from MUSIC. This is represented by the long dashed red line in Fig. 6.
The resulting optimal blast wave parameters are $n$ = 0.69 and $\beta_{max}$ =
0.76. Using these parameters and the blast wave parametrization yields a
momentum distribution for the ${\rm J}/\psi$ very similar to that of the
freeze out from the hydrodynamic simulations themselves (see Fig. 1 right)
verifying that the current procedure produces an optimum match. Comparing
these to the blast wave parameters using in Andronic:2019wva ;
Andronic:2021erx , $n$ = 0.85 and $\beta_{max}$ = 0.62, the former is about 20
% smaller and the latter 20 % larger than the optimum match. While the effect
of the increase in $\beta_{max}$ is more significant for the spectrum, both
changes go into the same direction of making the spectra harder (such as the
distributions shown in Figs. 2, 3) as compared to what is shown in
Andronic:2019wva .

As a further cross check we evaluate the equivalent freeze out volume
represented by the hypersurface elements in MUSIC. The weights
$\frac{1}{\gamma}u^{\mu}{\rm d}^{3}\Sigma_{\mu}$ give the volume relevant for
particle production in the collision center-of-momentum frame. When
integrating over the momenta to get the total yield, the Lorentz contraction
of this volume also needs to be considered, giving another factor of $\gamma$.
Summing up the contributions $u^{\mu}{\rm d}^{3}\Sigma_{\mu}$ thus gives an
estimate for the effective volume for particle production of about 5000 $\rm
fm^{3}$ very close to the result from the SHM fits Andronic:2021erx .

Figure 7: Left: Prompt D meson spectra for central Pb–Pb collisions including
feed down from strong decays of higher mass charmed hadrons using the
optimized blast wave parametrization and FastReso . They are compared to ALICE
data ALICE:2021rxa . Right: The same spectra normalized to those from pp
collisions scaled with the number of binary collisions.

The blast wave parameters were used as the input for a calculation of the
${\rm D}^{0}$ transverse momentum spectrum, considering the core contribution
with feed down from strong decays as in Andronic:2021erx  and corona
contribution based on the measured spectrum from pp collisions. The resulting
$p_{\text{T}}$ distributions are shown in Figure 7 in comparison with data
from the ALICE Collaboration. The uncertainties are determined as above for
the ${\rm J}/\psi$ spectra. Also shown in the right panel are the
corresponding distributions normalized to pp collisions scaled with the number
of binary collisions, i.e. $R_{\text{AA}}$. For $R_{\text{AA}}$ the
uncertainties introduced by the measured spectrum in pp collisions used for
the corona part cancel. It can be seen that the contribution from the
thermalized bulk of the medium is dominant at momenta below a few times the
particle mass, about 5 GeV/$c$. This region is in very good agreement with the
model. At higher momenta there is a deficit in the model. The additional
contribution could be from c-quarks not fully thermalized in the medium that
should show up at high $p_{\text{T}}$.

Figure 8: Left: Spectra of $\Lambda_{\rm c}$ for central Pb–Pb collisions
including feed down from strong decays of higher mass charmed baryons using
the optimized blast wave parametrization and FastReso . For this calculation
additional undiscovered charm baryon states were considered (see
Andronic:2021erx ). Also shown are ALICE data ALICE:2021bib . Right: The same
spectra normalized to those from pp collisions scaled with the number of
binary collisions.

The transverse momentum distributions of all other charmed hadrons can be
computed in an analogue way. We show here in Figure 8 the $p_{\text{T}}$
distribution calculated for $\Lambda_{\rm c}$ together with the corresponding
$R_{\text{AA}}$ distribution. Displayed are results from the SHMc with an
enhanced charmed baryon scenario as outlined in Andronic:2021erx . The
motivation is the observation by the ALICE collaboration of a very large
fragmentation fraction of charm into the $\Lambda_{\rm c}$ baryon in pp
collisions at the LHC ALICE:2021dhb . One explanation, proposed by He and Rapp
He:2019tik , was based on the observation that the charmed baryon spectrum
from the relativistic quark model contains many more charmed baryon states
than established experimentally. This is also the case for the charmed baryon
spectrum from lattice QCD Bazavov:2014yba . The hypothesis was made that these
states exist but were not yet discovered experimentally. This could be due to
a large width and therefore significant overlap. Including these theoretically
predicted states into their statistical model, He and Rapp found that they
could reproduce the enhanced fragmentation into $\Lambda_{\rm c}$. In the same
spirit, we have enhanced the contribution of charmed baryons by the equivalent
amount in our SHMc calculations and find that the $\Lambda_{\rm c}$ yield in
Pb–Pb collisions is about doubled while the D meson sector is practically
unchanged, for a correspondingly larger charm cross section (see
Andronic:2021erx  where both the regular and the enhanced charm baryon results
are shown).

The resulting $\Lambda_{\rm c}$ spectral distribution (Fig. 8) is compared to
recent data from ALICE and a similar picture emerges as for charmonia and D
mesons: The yield at low and moderate $p_{\text{T}}$ and therefore the bulk
yield are well reproduced by the current calculations, in particular
considering the uncertainties. At high $p_{\text{T}}$ there is again a deficit
in the model results. We also show in Figure 9 the ratio of $\Lambda_{\rm c}$
to D meson $p_{\text{T}}$ spectra. The experimental data from ALICE are very
well reproduced. In particular the maximum in this baryon to meson ratio comes
out naturally from the hydrodynamic description combined with statistical
hadronization without invoking an additional coalescence mechanism.

Figure 9: Ratio of $\Lambda_{\rm c}$ to ${\rm D}^{0}$ yields as function of
transverse momentum for central Pb–Pb collisions as shown separately in Figs:
7, 8 compared to ALICE data ALICE:2021bib ; ALICE:2021rxa .

##  7 Estimating systematic effects of charm space time distribution in the
medium

In the SHMc approach to charm hadronization we assume that the charm quarks,
at hadronization, are thermalized in terms of the momentum distributions, i.e.
the momentum distribution is determined by a (common) freeze out temperature
and the fluid expansion at this instance. This is supported by observables of
D mesons compared to transport models where a range is determined for the
spatial diffusion coefficient $1.5<2\pi D_{s}<4.5$ at the pseudo critical
temperature ALICE:2021rxa . From this a kinetic equilibration time
$\tau_{kin}$ = 2.5 - 7.6 fm/$c$ can be estimated. Recently, for the first time
spatial diffusion coefficients were derived in lattice QCD calculations with
light dynamical quarks Altenkort:2023oms . In the region between 1 and 2 times
the pseudo critical temperature the values were found to be significantly
smaller than in previous quenched lQCD calculations, implying very fast
hydrodynamization of charm quarks within 1-2 fm/$c$. This rapid thermalization
does however not imply that the spatial distribution of charm quarks at
hadronization is the same as that of light quarks and gluons.

While the initial energy density is distributed roughly following the areal
density of participants $N_{part}$, the charm quark density follows the areal
density of binary collisions $N_{coll}$ which is more compact. Therefore, at
early times the outermost region of the nuclear overlap region is less densely
populated with charm quarks as compared to energy, or equivalently, light
quarks and gluons.

We have no experimental observable at present that gives any information about
the spatial distribution of charm quarks at hadronization. In fact, since the
momentum equilibration takes some time (several collisions), it is plausible
that the expansion of charm quarks lags somewhat behind that of the overall
medium and even at late times the outermost regions of the freeze out
hypersurface may not be equally populated with charm quarks. In order to give
an indication how charmed hadron momentum distributions would be affected, we
present in the following two very schematic estimates basically removing
contributions from (i) early times or (ii) large radii at late times.

The distribution of hypersurface elements at freeze out is shown in Fig. 1 and
it is visible that the outermost regions of the nuclear overlap freeze out
within the first few fm/$c$. These are regions that, immediately after the
collision, are underpopulated in charm as compared to energy. And at these
early times it is also questionable that charm quarks are equilibrated in
terms of their momenta. To test the effect on D meson spectra, we remove
contributions from the hypersurface elements corresponding to the earliest 5 %
of the effective volume at freeze out. This happens to correspond to radii
where the density of $N_{coll}$ is less than half the density of $N_{part}$.
Since in statistical hadronization ${\rm J}/\psi$ mesons are formed from
uncorrelated charm and anticharm quarks, we cut off the earliest 15 % for this
study (considering the density of $N_{coll}^{2}$ vs density of $N_{part}$).
This corresponds to removing the hypersurface elements located within the
first 4.5 fm/$c$ (see Fig. 1) for D meson production and 6.7 fm/$c$ for ${\rm
J}/\psi$ production. We recalculate blast wave parameters matching the
remaining distribution and obtain for D mesons $n=0.647$ and $\beta_{\rm
max}=0.74$. In fact some of the elements with the largest velocity were
removed by this cut. The resulting spectra for D mesons and ${\rm J}/\psi$ are
shown in Fig. 10 as the short dashed lines in comparison to the full spectra.
Note that for this discussion only the core component is shown. Once can see a
visible but small effect in terms of softening of the $p_{\text{T}}$ spectra.

Figure 10: Comparison of core components of D meson and ${\rm J}/\psi$ spectra
with the default freeze out, to those calculated removing the early freeze-out
(short dashed lines) and to those removing, in the blast wave parametrization,
the outermost radii (long dashed lines). The spectra are normalized to the
same integral for easier comparison. For the ${\rm D}^{0}$ mesons, only the
direct contribution is shown.

To test the scenario that even at late times the charm quark spatial
distribution could be more compact than the overall fireball, we resort to the
blast wave distribution of velocities matched to the MUSIC freeze out
hypersurface and remove the outermost 1 fm. This simply implies $\beta_{max}$
= 0.69 instead 0.76. The corresponding spectra are also shown in Fig. 10 as
the long dashed lines. Now a considerable softening of the distributions is
observed. The most probable $p_{\text{T}}$ shifts downward for D mesons by 0.3
GeV/$c$ and for ${\rm J}/\psi$ by nearly 1 GeV/$c$.

While the two scenarios discussed are very schematic, it is obvious that they
both go in the direction of obtaining an even better match to experimental
data for ${\rm J}/\psi$, where we currently see a deviation, while much less
affecting the D meson spectrum, where we already see a good match between
model and data. It is also obvious that a slightly more compact final state
has a bigger effect than removal of early freeze out.

##  8 Summary and Outlook

In this paper we have explored the possibilities to describe quantitatively
and with essentially no free parameters the transverse momentum spectra and
anisotropic flow distributions for charmonia and charmed hadrons produced in
Pb–Pb collisions and measured with the ALICE detector at the CERN Large Hadron
Collider (LHC). The basis for this work is the realization that the SHM
originally developed to describe yields of hadrons composed of (u,d,s) quarks
can be expanded into the charm sector with no free additional parameters by
treating charm quarks as ’impurities’ in the hot and dense fireball with the
charm quark number determined by the measured total open charm cross section
in Pb–Pb or Au–Au collisions. Furthermore, the transverse dynamics in the
charm sector is more straightforward to implement than in the light quark
sector as we have evidence that all hadrons with heavy quarks are formed at or
very near the QCD phase boundary.

The specific focus of the present paper is indeed the transverse dynamics in
the charm quark sector which is further developed by coupling the SHMc to the
computer codes MUSIC and Fluid$u$M. For charmonia and hadrons with one charm
quark a rather quantitative description is obtained in the low $p_{\text{T}}$
sector for spectra of D-mesons, $\Lambda_{\rm c}$ baryons and charmonia. The
observed wide distribution in $p_{\text{T}}$ of anisotropic flow coefficients
v2 and v3 for charmonia is also well reproduced, but their magnitude is
generally somewhat over predicted.

This finding may be connected to a difference in spatial distribution between
light and charmed hadrons due to a different diffusion of light and heavy
quarks in the hot fireball. To make progress in this area two scenarios are
explored by modifying the distribution of charm quarks at late and early
times. The first results are promising but further development is needed to
make these studies quantitative.

On the theory side, to make progress requires more thorough modelling for
instance in the direction employed in Capellino:2023cxe . On the experimental
side, a measurement of Hanbury Brown-Twiss correlations between D mesons would
give experimental evidence about the charm quark spatial distribution.

Finally, a completely open field awaits us in the area of hadrons with 2 or 3
charm quarks. These will be experimentally explored in LHC Run 3 and Run 4
with the upgraded ALICE apparatus and ultimately with the proposed new ALICE 3
detector ALICE:2022wwr . It will be very exciting to test SHMc predictions in
the multi-charm sector and to explore hadronization and deconfinement
scenarios there.

###### Acknowledgements.

We would like to thank A. Kirchner, C. Shen and B. Schenke for their support
concerning Fluid$u$M and MUSIC. We acknowledge inspiring discussions with S.
Floerchinger, B. Friman, A. Mazeliauskas and K. Redlich. This work was
performed within and supported by the DFG Collaborative Research Centre ”SFB
1225 (ISOQUANT)”. V. V. gratefully acknowledges funding from the Knut and
Alice Wallenberg Foundation (the CLASH project).

## References

  * (1) ALICE collaboration, _The ALICE experiment – A journey through QCD_ , 2211.04384.
  * (2) A. Andronic, P. Braun-Munzinger, K. Redlich and J. Stachel, _Statistical hadronization of charm in heavy ion collisions at SPS, RHIC and LHC_ , _Phys. Lett. B_ 571 (2003) 36 [nucl-th/0303036].
  * (3) A. Andronic, P. Braun-Munzinger, K. Redlich and J. Stachel, _Decoding the phase structure of QCD via particle production at high energy_ , _Nature_ 561 (2018) 321 [1710.09425].
  * (4) A. Andronic, P. Braun-Munzinger, M.K. Köhler, A. Mazeliauskas, K. Redlich, J. Stachel et al., _The multiple-charm hierarchy in the statistical hadronization model_ , _JHEP_ 07 (2021) 035 [2104.12754].
  * (5) F. Gross et al., _50 Years of Quantum Chromodynamics_ , 2212.11107.
  * (6) P. Braun-Munzinger and J. Stachel, _(Non)thermal aspects of charmonium production and a new look at J / psi suppression_ , _Phys. Lett. B_ 490 (2000) 196 [nucl-th/0007059].
  * (7) A. Andronic, P. Braun-Munzinger, M.K. Köhler, K. Redlich and J. Stachel, _Transverse momentum distributions of charmonium states with the statistical hadronization model_ , _Phys. Lett. B_ 797 (2019) 134836 [1901.09200].
  * (8) H. Song and U.W. Heinz, _Causal viscous hydrodynamics in 2+1 dimensions for relativistic heavy-ion collisions_ , _Phys. Rev. C_ 77 (2008) 064901 [0712.3715].
  * (9) P.J. Siemens and J.O. Rasmussen, _Evidence for a blast wave from compress nuclear matter_ , _Phys. Rev. Lett._ 42 (1979) 880.
  * (10) E. Schnedermann, J. Sollfrank and U.W. Heinz, _Thermal phenomenology of hadrons from 200-A/GeV S+S collisions_ , _Phys. Rev. C_ 48 (1993) 2462 [nucl-th/9307020].
  * (11) B. Schenke, S. Jeon and C. Gale, _Elliptic and triangular flow in event-by-event (3+1)D viscous hydrodynamics_ , _Phys. Rev. Lett._ 106 (2011) 042301 [1009.3244].
  * (12) S. Floerchinger, E. Grossi and J. Lion, _Fluid dynamics of heavy ion collisions with mode expansion_ , _Phys. Rev. C_ 100 (2019) 014905 [1811.01870].
  * (13) A. Andronic, P. Braun-Munzinger and J. Stachel, _Hadron production in central nucleus-nucleus collisions at chemical freeze-out_ , _Nucl. Phys. A_ 772 (2006) 167 [nucl-th/0511071].
  * (14) ALICE collaboration, _J/ $\psi$ Elliptic Flow in Pb-Pb Collisions at $\sqrt{s_{\rm NN}}$ = 2.76 TeV_, _Phys. Rev. Lett._ 111 (2013) 162301 [1303.5880].
  * (15) ALICE collaboration, _J/ $\psi$ elliptic and triangular flow in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV_, _JHEP_ 10 (2020) 141 [2005.14518].
  * (16) ALICE collaboration, _D meson elliptic flow in non-central Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 2.76 TeV_, _Phys. Rev. Lett._ 111 (2013) 102301 [1305.2707].
  * (17) K. Zhou, N. Xu, Z. Xu and P. Zhuang, _Medium effects on charmonium production at ultrarelativistic energies available at the CERN Large Hadron Collider_ , _Phys. Rev. C_ 89 (2014) 054911 [1401.5845].
  * (18) X. Du and R. Rapp, _Sequential Regeneration of Charmonia in Heavy-Ion Collisions_ , _Nucl. Phys. A_ 943 (2015) 147 [1504.00670].
  * (19) M. He, B. Wu and R. Rapp, _Collectivity of J/ $\psi$ Mesons in Heavy-Ion Collisions_, _Phys. Rev. Lett._ 128 (2022) 162301 [2111.13528].
  * (20) A. Andronic, P. Braun-Munzinger, K. Redlich and J. Stachel, _Statistical hadronization of heavy quarks in ultra-relativistic nucleus-nucleus collisions_ , _Nucl. Phys. A_ 789 (2007) 334 [nucl-th/0611023].
  * (21) ALICE collaboration, _Prompt D 0, D+, and D∗+ production in Pb–Pb collisions at $\sqrt{s_{\mathrm{NN}}}$ = 5.02 TeV_, _JHEP_ 01 (2022) 174 [2110.09420].
  * (22) ALICE collaboration, _Measurements of inclusive J/ $\psi$ production at midrapidity and forward rapidity in Pb$-$Pb collisions at $\sqrt{s_{\mathrm{NN}}}$ = 5.02 TeV_, 2303.13361.
  * (23) M.I. Gorenstein, A. Kostyuk, H. Stoecker and W. Greiner, _Statistical coalescence model with exact charm conservation_ , _Phys. Lett. B_ 509 (2001) 277 [hep-ph/0010148].
  * (24) P. Braun-Munzinger, K. Redlich and J. Stachel, _Particle production in heavy ion collisions, Invited review for Quark Gluon Plasma 3, eds. R. C. Hwa and X.-N. Wang, World Scientific Publishing_ , nucl-th/0304013.
  * (25) F. Cooper and G. Frye, _Comment on the Single Particle Distribution in the Hydrodynamic and Statistical Thermodynamic Models of Multiparticle Production_ , _Phys. Rev. D_ 10 (1974) 186.
  * (26) A. Kirchner, E. Grossi and S. Floerchinger, _Cooper-Frye spectra of hadrons with viscous corrections including feed down from resonance decays_ , 2308.10616.
  * (27) B. Schenke, C. Shen and P. Tribedy, _Running the gamut of high energy nuclear collisions_ , _Phys. Rev. C_ 102 (2020) 044905 [2005.14682].
  * (28) B. Schenke, P. Tribedy and R. Venugopalan, _Event-by-event gluon multiplicity, energy density, and eccentricities in ultrarelativistic heavy-ion collisions_ , _Phys. Rev. C_ 86 (2012) 034908 [1206.6805].
  * (29) B. Schenke, P. Tribedy and R. Venugopalan, _Fluctuating Glasma initial conditions and flow in heavy ion collisions_ , _Phys. Rev. Lett._ 108 (2012) 252301 [1202.6646].
  * (30) D. Devetak, A. Dubla, S. Floerchinger, E. Grossi, S. Masciocchi, A. Mazeliauskas et al., _Global fluid fits to identified particle transverse momentum spectra from heavy-ion collisions at the Large Hadron Collider_ , _JHEP_ 06 (2020) 044 [1909.10485].
  * (31) J.S. Moreland, J.E. Bernhard and S.A. Bass, _Alternative ansatz to wounded nucleon and binary collision scaling in high-energy nuclear collisions_ , _Phys. Rev. C_ 92 (2015) 011901 [1412.4708].
  * (32) F. Capellino, A. Dubla, S. Floerchinger, E. Grossi, A. Kirchner and S. Masciocchi, _Fluid-dynamics of charm quarks in the quark-gluon plasma_ , 2307.14449.
  * (33) A. Andronic, P. Braun-Munzinger, K. Redlich and J. Stachel, _Statistical Hadronization of $b$-quarks in Pb–Pb Collisions at LHC Energy: A Case for Partial Equilibration of $b$-quarks?_, _Acta Phys. Polon. Supp._ 16 (2023) 107 [2209.14562].
  * (34) M. Cacciari, M. Greco and P. Nason, _The $p_{\text{T}}$ spectrum in heavy flavor hadroproduction_, _JHEP_ 05 (1998) 007 [hep-ph/9803400].
  * (35) C. Bierlich et al., _A comprehensive guide to the physics and usage of PYTHIA 8.3_ , 2203.11601.
  * (36) ALICE collaboration, _Production of charged pions, kaons, and (anti-)protons in Pb-Pb and inelastic $pp$ collisions at $\sqrt{s_{NN}}$ = 5.02 TeV_, _Phys. Rev. C_ 101 (2020) 044907 [1910.07678].
  * (37) ALICE collaboration, _$\psi(2S)$ suppression in Pb-Pb collisions at the LHC_, 2210.08893.
  * (38) B. Alver and G. Roland, _Collision geometry fluctuations and triangular flow in heavy-ion collisions_ , _Phys. Rev. C_ 81 (2010) 054905 [1003.0194].
  * (39) ATLAS collaboration, _Measurement of azimuthal anisotropy of muons from charm and bottom hadrons in Pb+Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV with the ATLAS detector_, _Phys. Lett. B_ 807 (2020) 135595 [2003.03565].
  * (40) A. Mazeliauskas, S. Floerchinger, E. Grossi and D. Teaney, “FastReso – program for computing irreducible components of the particle distribution from direct resonance decays.” GitHub repository, 2019.
  * (41) W. Florkowski, _Phenomenology of Ultra-Relativistic Heavy-Ion Collisions_ , World Scientific Publishing Co. (3, 2010).
  * (42) ALICE collaboration, _Constraining hadronization mechanisms with $\Lambda$c+/D0 production ratios in Pb–Pb collisions at sNN=5.02 TeV_, _Phys. Lett. B_ 839 (2023) 137796 [2112.08156].
  * (43) ALICE collaboration, _Charm-quark fragmentation fractions and production cross section at midrapidity in pp collisions at the LHC_ , _Phys. Rev. D_ 105 (2022) L011103 [2105.06335].
  * (44) M. He and R. Rapp, _Charm-Baryon Production in Proton-Proton Collisions_ , _Phys. Lett. B_ 795 (2019) 117 [1902.08889].
  * (45) A. Bazavov et al., _The melting and abundance of open charm hadrons_ , _Phys. Lett. B_ 737 (2014) 210 [1404.4043].
  * (46) HotQCD collaboration, _Heavy Quark Diffusion from 2+1 Flavor Lattice QCD with 320 MeV Pion Mass_ , _Phys. Rev. Lett._ 130 (2023) 231902 [2302.08501].
  * (47) ALICE collaboration, _Letter of intent for ALICE 3: A next-generation heavy-ion experiment at the LHC_ , 2211.02491.