# Hyper-order baryon number fluctuations at finite temperature and density

Wei-jie Fu  School of Physics, Dalian University of Technology, Dalian,
116024, P.R. China  Xiaofeng Luo  Key Laboratory of Quark & Lepton Physics
(MOE) and Institute of Particle Physics, Central China Normal University,
Wuhan 430079, China  Jan M. Pawlowski  Institut für Theoretische Physik,
Universität Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany  ExtreMe
Matter Institute EMMI, GSI, Planckstraße 1, D-64291 Darmstadt, Germany  Fabian
Rennecke  Brookhaven National Laboratory, Upton, NY 11973, USA  Rui Wen
School of Physics, Dalian University of Technology, Dalian, 116024, P.R. China
Shi Yin  School of Physics, Dalian University of Technology, Dalian, 116024,
P.R. China

###### Abstract

Fluctuations of conserved charges are sensitive to the QCD phase transition
and a possible critical endpoint in the phase diagram at finite density. In
this work, we compute the baryon number fluctuations up to tenth order at
finite temperature and density. This is done in a QCD-assisted effective
theory that accurately captures the quantum- and in-medium effects of QCD at
low energies. A direct computation at finite density allows us to assess the
applicability of expansions around vanishing density. By using different
freeze-out scenarios in heavy-ion collisions, we translate these results into
baryon number fluctuations as a function of collision energy. We show that a
non-monotonic energy dependence of baryon number fluctuations can arise in the
non-critical crossover region of the phase diagram. Our results compare well
with recent experimental measurements of the kurtosis and the sixth-order
cumulant of the net-proton distribution from the STAR collaboration. They
indicate that the experimentally observed non-monotonic energy dependence of
fourth-order net-proton fluctuations is highly non-trivial. It could be an
experimental signature of an increasingly sharp chiral crossover and may
indicate a QCD critical point. The physics implications and necessary upgrades
of our analysis are discussed in detail.

##  I Introduction

Some of the most challenging questions of heavy-ion physics are related to the
transition from the early, non-equilibrium, state of quarks and gluons to the
final hadronic states after chemical freeze out, which is observed in
experiments. Unravelling this dynamics necessitates a thorough grasp on the
physics in the QCD phase structure close to the confinement-deconfinement and
chiral transitions. This regime is strongly correlated with highly non-trivial
dynamics. Understanding this part of the phase structure, including the
location and dynamics of a potential critical end point (CEP), plays a pivotal
role in understanding phases of strongly interacting nuclear matter under
extreme conditions. For works on the phase structure of QCD, covering
experiment and theory see, e.g., Stephanov (2006a); Friman  _et al._ (2011a);
Blaschke  _et al._ (2016); Luo and Xu (2017); Dainese  _et al._ (2019); Bzdak
_et al._ (2020); Fischer (2019); Fu  _et al._ (2020a); Bazavov  _et al._
(2020); Borsanyi  _et al._ (2020); Lü  _et al._ (2017); Yang  _et al._ (2013);
Sako (2019), where theory covers first principles functional approaches and
lattice simulations.

Fluctuations of conserved charges are very sensitive to the physics of the
strongly correlated regime that governs the transition from the quark-gluon
plasma (QGP) to the hadronic phase. They provide detailed information on the
underlying dynamics. This includes, but is not limited to, possible
experimental signatures of a CEP Luo and Xu (2017). It has for example been
proposed in Stephanov  _et al._ (1999); Stephanov (2009, 2011) that non-
monotonic variations of conserved charge fluctuations as functions of the beam
energy can arise from critical physics in the vicinity of a CEP. During the
last decade, significant fluctuation measurements have been performed in the
first phase of the Beam Energy Scan (BES-I) program at the Relativistic Heavy
Ion Collider (RHIC), involving various cumulants of net-proton, net-charge and
net-kaon multiplicity distributions Adamczyk  _et al._ (2014a, b); Luo (2015);
Adamczyk  _et al._ (2018); Adam  _et al._ (2019). Remarkably, very recently
the STAR collaboration has reported the first evidence of a non-monotonic
variation in the kurtosis (multiplied by the variance) of the net-proton
number distribution as a function of the collision energy with $3.1\,\sigma$
significance for central collisions Adam  _et al._ (2020). The measurements
have been extended to the sixth-order cumulants of net-proton and net-charge
distributions, for preliminary results from STAR see Nonaka (2020); Pandav
(2020).

Recent first-principle QCD calculations at finite temperature and density,
within both the functional renormalisation group (fRG) and Dyson-Schwinger
equations (DSE), show that the transition from the QGP to the hadronic phase
is a crossover which becomes sharper with increasing baryon chemical
potential, $\mu_{B}$, for $\mu_{B}/T\lesssim 4$ Fischer (2019); Fu  _et al._
(2020a); Braun  _et al._ (2020a); Isserstedt  _et al._ (2019); Gao and
Pawlowski (2020a, b). Beyond this region, a CEP might occur, but quantitative
reliability of the theory computations cannot be guaranteed within the present
approximations Fischer (2019); Fu  _et al._ (2020a); Braun  _et al._ (2020a);
Gao and Pawlowski (2020b). In addition, critical physics may only be
observable in a very small region around the CEP, see e.g. Schaefer and
Wambach (2007). Since the available RHIC data is limited to $\mu_{B}/T\lesssim
3$, it is important to understand how conserved charge fluctuations are
affected by the increasingly sharp crossover away from a regime with critical
scaling.

To address these open questions related to the physics of strong correlations
in the QCD phase diagram, we study in detail the $T$\- and
$\mu_{B}$-dependence of net-baryon number fluctuations in the range
$\mu_{B}/T\lesssim 3$. We present results for fluctuations up to tenth order
(where we refer to everything above fourth order as _hyper-order_), including
comparisons to available results from RHIC Adam  _et al._ (2020); Nonaka
(2020); Pandav (2020) and predictions for the beam energy dependence of
fluctuations where no experimental results are available yet.

To facilitate the comparison between theory and experiment, $T$ and $\mu_{B}$
can be mapped onto the beam energy per nucleon, $\sqrt{s_{\rm NN}}$, via
phenomenological freeze-out curves Braun-Munzinger  _et al._ (2004a). While
these curves are expected to be close to the QCD crossover at large beam
energies (corresponding to small $\mu_{B}$) Braun-Munzinger  _et al._ (2004b),
they may move away from the transition region at low energies (i.e. larger
$\mu_{B}$) Floerchinger and Wetterich (2012). This can affect the beam energy
dependence of particle number fluctuations, and requires a detailed
understanding of the physics also outside the critical region.

Aside from their phenomenological relevance, net-baryon number fluctuations at
finite $\mu_{B}$ can also be used to assess the reliability of extrapolations
of thermodynamic quantities to finite $\mu_{B}$ based on a Taylor expansion at
$\mu_{B}=0$. Such a strategy is commonly used in lattice QCD simulations,
where a sign problem prevents direct simulations at finite $\mu_{B}$, see e.g.
Bazavov  _et al._ (2012); Borsanyi  _et al._ (2013, 2014); Bazavov  _et al._
(2017a, b); Borsanyi  _et al._ (2018); Bazavov  _et al._ (2020). By comparing
the results of direct computations at finite $\mu_{B}$ to the ones obtained
from extrapolations at $\mu_{B}=0$, we study the range of validity of a Taylor
expansion at a given order self-consistently. Understanding the limitations of
such an extrapolation is also relevant for phenomenologically constructed
equations of state, as, e.g., in Parotto  _et al._ (2020), where the non-
critical physics at finite $\mu_{B}$ crucially rely on this extrapolation.

All this is addressed within a QCD-assisted low-energy effective field theory
(LEFT) which is described in detail in the next section. We use first-
principles QCD-results on the $T$-dependence of the kurtosis and the
$\mu_{B}$-dependence of the chiral phase boundary to map the in-medium scales
of the LEFT onto QCD. This improves the reliability of our predictions, in
particular at finite $\mu_{B}$. Non-perturbative quantum-, thermal- and
density fluctuations are taken into account with the functional
renormalisation group (fRG). This work therefore is a significant upgrade of
previous work in Fu and Pawlowski (2015, 2016); Fu  _et al._ (2016), where
net-baryon number fluctuations up to fourth order have been studied. The
present QCD-assisted LEFT approach has various advantages. Most importantly,
it is directly embedded in QCD as the relevant low-energy degrees of freedom
emerge dynamically from systematically integrating-out the fast partonic modes
of QCD Mitter  _et al._ (2015); Braun  _et al._ (2016); Rennecke (2015a);
Cyrol  _et al._ (2018a); Fu  _et al._ (2020a). In addition, this approach
allows us to capture both critical and non-critical effects in the QCD phase
diagram. This entails in particular that our results agree with the results of
lattice QCD at small $\mu_{B}$ and show the correct universal behaviour of QCD
in the vicinity of the CEP, i.e. $3d$ Ising universality.

Concerning the existence and location of the latter we add that the first-
principles results in Fischer (2019); Fu  _et al._ (2020a); Braun  _et al._
(2020a); Isserstedt  _et al._ (2019); Gao and Pawlowski (2020a, b) include a
CEP in a region of $450\,\mathrm{MeV}\lesssim\mu_{B}\lesssim
650\,\mathrm{MeV}$ and therefore outside the regime of quantitative
reliability of these computations. This suggests that the experimental
detection of a CEP requires explorations of the high-$\mu_{B}$ in the region
with $\mu_{B}/T_{c}\gtrsim 4$. Moreover, the direct experimental measurement
of the CEP may be very challenging as it requires very high statistics, and
predictions that signal critical dynamics can be further complicated by non-
equilibrium effects. In the present work we shall therefore also outline how
the location of a CEP could be constrained based on data in the crossover
region, without the necessity of observing critical scaling. On the
theoretical side this asks for first principles QCD studies for
$\mu_{B}/T\gtrsim 4$. In turn, a first experimental step towards this goal is
the solidification of experimental observation of the non-monotonic energy
dependence of fourth-order net-proton fluctuations. Both is safely beyond the
scope of the present work.

This paper is organised as follows: In Section II we give a brief introduction
to the fRG-approach to QCD and low-energy effective theories, including their
mutual relationship. Thermodynamics and the hyper-order baryon number
fluctuations are discussed in Section III. In Section IV, we first introduce a
systematic scale-matching procedure between QCD and the low-energy effective
theory. We then present our numerical results and compare them to lattice QCD
simulations and experimental measurements. A summary with conclusions is given
in Section V. Technical details regarding the flow equations are presented in
the appendices.

##  II QCD and emergent low-energy effective theories

Figure 1: Diagrammatic representation of the QCD flow equation. The lines
stand for the full propagators of gluon, ghost, quark, and mesons,
respectively. The arrows in quark and meson lines indicate the quark number
(baryon number) flow. The crossed circles represent the infrared regulators.

At low momentum scales the quark-gluon dynamics of QCD successively decouple
due to the QCD mass gap and spontaneous chiral symmetry breaking. This
decoupling also applies to most dynamical (hadronic) low energy degrees of
freedom at even lower energies, finally leaving us with dynamical pions and
hence with chiral perturbation theory. Indeed, this successive decoupling is
at the root of the success of chiral effective field theory.

The functional renormalisation group approach to QCD with its successive
integrating-out of momentum modes is ideally suited to follow and study this
decoupling. Diagrammatically, this is already seen within the flow equation
for the QCD effective action, depicted in Figure 1. The different lines stand
for the full non-perturbative propagators of gluons, ghosts, quarks and
emergent low-energy degrees of freedom (hadrons in our case), where the loop
momentum $q$ is restricted by the infrared cutoff scale $k$, $q^{2}\lesssim
k^{2}$. In this setup, emergent bound states can be incorporated
systematically by dynamical hadronisation Gies and Wetterich (2002, 2004);
Pawlowski (2007); Floerchinger and Wetterich (2009). For quantitative QCD
applications in the vacuum see Braun  _et al._ (2016); Rennecke (2015a);
Mitter  _et al._ (2015); Cyrol  _et al._ (2018a), for further conceptual
developments and the application to the QCD phase structure important for the
present work see Fu  _et al._ (2020a). The decoupling is apparent in this
framework as the propagators carry the mass gaps $m_{\textrm{gap}}$ of gluons
and quarks and for cutoff scales $k\ll m_{\textrm{gap}}$ of a given field the
respective loop tends towards zero.

More importantly, in this way the emergent low-energy effective theory is
naturally embedded in QCD, and its ultraviolet parameters (at $\Lambda\lesssim
1$ GeV) as well as further input may be directly computed from QCD, leading to
QCD-assisted low-energy effective theories. We emphasise that this procedure
does not lead to a unique LEFT. The dynamical degrees of freedom of QCD-
assisted LEFTs at $\Lambda\lesssim 1$ GeV depend on the dynamical
hadronisation procedure applied within QCD-flows. This setup and the QCD-
embedding entail that, provided the relevant quantum, thermal and density
fluctuations of low energy QCD are taken into account in the QCD-assisted LEFT
at hand, all QCD-assisted LEFTs encode the same physics, namely that of low
energy QCD.

This leads to an equivalence relation of QCD-assisted Polyakov-loop–enhanced
NJL-type LEFTs (PNJL), Polyakov-loop–enhanced QM LEFTs (PQM) and variations
including higher meson multiplets and/or diquarks and baryons. We emphasise
again that this equivalence relation only holds if low energy quantum, thermal
and density fluctuations are taken into account. For more details see in
particular Fu  _et al._ (2020a), and the recent review Dupuis  _et al._
(2020). Most prominently this embedding has been used for determining the
temperature-dependence of the Polyakov loop potential, see Haas  _et al._
(2013); Herbst  _et al._ (2014). This setup was then applied to the
computation of fluctuations in Fu and Pawlowski (2016, 2015); Fu  _et al._
(2016); Wen  _et al._ (2019); Yin  _et al._ (2019).

In summary this entails, that for sufficiently small momenta $k$, temperatures
$T$, and also density or quark chemical potential $\mu_{q}$, the gluon (and
ghost) loop in Figure 1 decouple from the dynamics, and only provide a non-
trivial glue background at finite temperature and chemical potential. The
latter is taken into account with the Polyakov loop potential discussed in
detail below.

In the present work, we build upon previous investigations of the skewness and
kurtosis of baryon number distributions Fu and Pawlowski (2015, 2016); Fu  _et
al._ (2016), and baryon-strangeness correlations Fu  _et al._ (2020b, 2019)
within QCD-assisted LEFTs with the fRG. The present LEFT is an upgrade of
those used in the works above, and includes the quantum, thermal and density
dynamics of quarks, pions and the sigma mode in a Polyakov loop background. It
is a QCD-assisted PQM. As argued above, for low enough chemical potential,
this model sufficiently close to QCD, and leads to results that are
independent of the LEFT at hand.

For further investigations of fluctuation observables within the fRG-approach
to low-energy effective theories see e.g. Skokov  _et al._ (2010, 2011);
Friman  _et al._ (2011b); Morita  _et al._ (2015); Almasi  _et al._ (2017a),
the Dyson-Schwinger approach has been used in e.g. Xin  _et al._ (2014);
Isserstedt  _et al._ (2019), for mean-field investigations see e.g. Fu  _et
al._ (2010); Fu and Wu (2010); Karsch  _et al._ (2011); Schaefer and Wagner
(2012); Li  _et al._ (2019a). These functional works can be further adjusted
and benchmarked with results from lattice QCD simulations Bazavov  _et al._
(2012); Borsanyi  _et al._ (2013, 2014); Bazavov  _et al._ (2017a, b);
Borsanyi  _et al._ (2018); Bazavov  _et al._ (2020), at high temperatures,
$T\gtrsim T_{c}$, and vanishing $\mu_{B}$. In turn, at finite $\mu_{B}$, and
in particular for $\mu_{B}/T_{c}\gtrsim 3$, lattice simulations are obstructed
by the sign problem.

###  II.1 2-flavour setup

For the physics of fluctuations we are interested in scales below
approximately $1$ GeV. We restrict ourselves to $k\lesssim 700$ MeV and
temperatures and quark chemical potentials $T,\mu_{q}\lesssim 200$ MeV. In
this regime the only relevant quarks are the light quarks $q=(u,d)$ and the
strange quark $s$. The latter, while changing the momentum-scale running of
the correlation functions, has subleading effects on the form of the
fluctuations. Hence, the effect of the momentum-scale running induced by
strange fluctuations will be mimicked here by an appropriate scale-matching
detailed in Section II.2.

We also include the lowest lying hadronic resonances, the pion
${\bm{\pi}}=(\pi^{\pm},\pi^{0})$, and, for symmetry reasons, the scalar
resonance $\sigma$ as effective low energy degrees of freedom. Within QCD
flows these fields are emergent low energy degrees of freedom at cutoff scales
$k\simeq 1$ GeV, that are taken care of with dynamical hadronisation in e.g.
Fu  _et al._ (2020a). At the present low energy scales $k\leq 700$ MeV, they
are fully dynamical, and hence are part of the effective action at the initial
cutoff scale. The other members of the lowest lying multiplet as well as
further hadronic resonances produce rather subleading contributions to the
offshell dynamics and hence are dropped. The mesonic fields are stored in an
O(4) scalar field $\phi=(\sigma,{\bm{\pi}})$ with the corresponding chiral
invariant $\rho=\phi^{2}/2$.

Quantum, thermal and density fluctuations with scales $k\lesssim\Lambda=700$
MeV are taken into account within the fRG, whose dynamics is now reduced to
the last two loops in Figure 1. The respective effective action of QCD in the
low energy regime is approximated by

$\displaystyle\Gamma_{k}=$
$\displaystyle\int_{x}\bigg{\\{}Z_{q,k}\bar{q}\Big{[}\gamma_{\mu}\partial_{\mu}-\gamma_{0}(\mu+igA_{0})\Big{]}q+\frac{1}{2}Z_{\phi,k}(\partial_{\mu}\phi)^{2}$
$\displaystyle\hskip
14.22636pt+h_{k}\,\bar{q}\,\left(\tau^{0}\sigma+\bm{\tau}\cdot\bm{\pi}\right)\,q+V_{k}(\rho,A_{0})-c\sigma\bigg{\\}}\,,$
(1)

with $\int_{x}=\int_{0}^{1/T}dx_{0}\int d^{3}x$ and
$\tau=1/2(\mathbb{1},i\gamma_{5}\bm{\sigma})$. We assume isospin symmetry and
the corresponding chemical potential flavour-matrix is given by
$\mu=\text{diag}(\mu_{q},\mu_{q})=\frac{1}{3}\text{diag}(\mu_{B},\mu_{B})$.
$Z_{q,k}$ and $Z_{\phi,k}$ are the wave function renormalisations for the
light quarks and the meson respectively. Further running couplings considered
are the Yukawa coupling $h_{k}$, the scattering between quarks and mesons, as
well as the effective potential $V_{k}(\rho,A_{0})$, that describes the multi-
scattering of mesons in the non-trivial glue background present at finite
temperature and chemical potential.

The flow equation for the effective action Equation 1, and that for
$V_{k},h_{k},Z_{\phi,q}$ is described in Appendix A and Appendix B. The
initial condition for $\Gamma_{k}$ at the initial cutoff scale $k=700$ MeV is
described in Appendix C.

The potential $V_{k}(\rho,A_{0})$ has contributions
$V_{\mathrm{glue},k}(A_{0})$ from offshell glue fluctuations (first two
diagrams in Figure 1), and contributions $V_{\mathrm{mat},k}(\rho,A_{0})$ from
the quark loop (third diagram in Figure 1). This leads us to

$\displaystyle V_{k}(\rho,A_{0})=$ $\displaystyle
V_{\mathrm{glue},k}(A_{0})+V_{\mathrm{mat},k}(\rho,A_{0})\,,$ (2)

The first contribution is typically reformulated in terms of the Polyakov loop
$L(A_{0})$, while the latter is directly computed from the present low energy
flow. This allows us to trade the $A_{0}$-dependence for that of the traced
Polyakov loop, $L(A_{0}),\bar{L}(A_{0})$, see Appendix D, Equation 49, leading
us to to the final form of our potential,

$\displaystyle V_{k}(\rho,L,\bar{L})=V_{k}(\rho,A_{0})\,.$ (3)

More details can be found in Appendix D.

In conclusion, the QCD-assisted LEFT described above and used in the present
work, is a PQM-type model, e.g. Schaefer  _et al._ (2007, 2010); Skokov  _et
al._ (2010); Herbst  _et al._ (2011); Skokov  _et al._ (2011); Karsch  _et
al._ (2011); Schaefer and Wagner (2012); Morita  _et al._ (2011); Skokov  _et
al._ (2012); Mintz  _et al._ (2013); Haas  _et al._ (2013); Herbst  _et al._
(2013, 2014); Fu and Pawlowski (2016, 2015); Fu  _et al._ (2016); Stiele and
Schaffner-Bielich (2016); Sun  _et al._ (2018); Fu  _et al._ (2020b, 2019);
Wen  _et al._ (2019); Yin  _et al._ (2019); Hansen  _et al._ (2020). Quantum,
thermal and density fluctuations below $\Lambda=700$ MeV are taken account
with the functional renormalisation group, and the setup is well-embedded in
functional QCD. As argued above, within the present, and analogous, elaborate
approximation, the respective results (for fluctuation observables) for all
QCD-assisted LEFTs match those of QCD for sufficiently low density. Therefore,
we will refer to this model from now on as generic QCD-assisted LEFT.

###  II.2 $2+1$-favour scale-matching in $2$-flavour QCD

The current QCD-assisted LEFT setup enables us to compute thermodynamic
observables and in particular hyper-order baryon number fluctuations. However,
as already briefly discussed in Section II.1, we have dropped the dynamics of
the strange quark. While we expect sub-dominant effects on hyper-fluctuations,
the $s$-quark influences the momentum running of the correlations in the
ultraviolet.

Importantly, in Fu  _et al._ (2020a) it has been observed on the basis of
genuine $N_{f}=2$ and $N_{f}=2+1$ flavour computations in QCD, that the latter
effect is well approximated by a respective universal scale-matching of the
2-flavour results even in QCD. Such a scale-matching has already led to a
quantitative agreement of thermodynamics and kurtosis within the current LEFT
setup with lattice results, see Fu and Pawlowski (2016, 2015); Fu  _et al._
(2016). Thus, the present scale-matching entails that we use information on
the $T$\- and $\mu_{B}$-dependence of well-determined quantities in QCD. This
leads to an improved reliability of our results of finite $T$ and $\mu_{B}$,
as in-medium effects in the QCD-assisted LEFT are directly connected to in-
medium effects in QCD.

####  II.2.1 2- to 2+1-flavour scale-matching in QCD

Given its relevance for the predictive power of the present LEFT within a QCD
scale-matching procedure we briefly recall the respective results in Fu  _et
al._ (2020a): There, the phase boundaries of 2- and 2+1-flavour QCD have been
computed within the fRG approach. These results allow us to evaluate the
reliability of even linear scale-matchings of temperatures and chemical
potentials in 2- and 2+1-flavour QCD introduced by

$\displaystyle T^{(N_{f}=2)}$ $\displaystyle=c_{{}_{T}}\,T^{(N_{f}=2+1)}\,,$
$\displaystyle\mu_{B}^{(N_{f}=2)}$
$\displaystyle=c_{\mu_{B}}\,\mu_{B}^{(N_{f}=2+1)}\,.$ (4)

With such a linear scale-matching the scaling factors $c_{{}_{T}},c_{\mu_{B}}$
can be determined by evaluating the relations at a specific temperature and
chemical potential.

For the thermal scale-matching we naturally take $(T,\mu_{B})=(T_{c},0)$, the
crossover temperature at vanishing chemical potential. In Fu  _et al._ (2020a)
the crossover temperatures have been determined with thermal susceptibilities
of the renormalised light chiral condensate. Then, the linear rescaling of the
2-flavour chiral crossover temperature to the 2+1-flavour crossover
temperature is done with

$\displaystyle
T_{c}^{(N_{f}=2)}=c_{T}^{\textrm{QCD}}\,T_{c}^{(N_{f}=2+1)}\,,\qquad
c_{T}^{\textrm{QCD}}=1.1\,.$ (5)

For the matching of the chemical potentials we use the curvature $\kappa$ of
the phase boundary at vanishing $\mu_{B}=0$,

$\displaystyle\frac{T_{c}(\mu_{B})}{T_{c}}$
$\displaystyle=1-\kappa\left(\frac{\mu_{B}}{T_{c}}\right)^{2}+\lambda\left(\frac{\mu_{B}}{T_{c}}\right)^{4}+\cdots\,.$
(6)

Adjusting the 2-flavour curvature $-\kappa\,\mu_{B}^{2}/T_{c}^{2}$ to the
2+1-flavour one leads us to the relation

$\displaystyle c^{\textrm{QCD}}_{\mu_{B}}$
$\displaystyle=c^{\textrm{QCD}}_{{}_{T}}\left(\frac{\kappa^{(N_{f}=2+1)}}{\kappa^{(N_{f}=2)}}\right)^{1/2}\,,\qquad
c_{\mu_{B}}^{\textrm{QCD}}=0.99\,.$ (7)

The value $c_{\mu_{B}}^{\textrm{QCD}}\approx 1$ entails that the change in the
curvature coefficient $\kappa$ is balanced by that of the temperature.

Figure 2: Phase boundaries of 2- and 2+1-flavour QCD from Fu  _et al._ (2020a)
with a 2+1-flavour scale-matching of the 2-flavour data at the crossover
temperature and $\mu_{B}=0$. The bands denote the width of the chiral
crossover. The scale-matched 2-flavour phase boundary agrees quantitatively
with the genuine 2+1-flavour one including the location of the critical end
point. The dashed line at $\mu_{B}/T=4$ constitutes the reliability bound of
the computations in Fu  _et al._ (2020a) based on the potential emergence of
new degrees of freedom discussed in Fu  _et al._ (2020a); Braun  _et al._
(2020a); Fischer (2019). The dashed lines at $\mu_{B}/T=2,3$ are reliability
estimates of lattice results as well as old ones from functional approaches.

The fourth-order expansion coefficient $\lambda$ is found to be very small in
both functional, Fu  _et al._ (2020a); Gao and Pawlowski (2020b, a) as well as
lattice computations, Bazavov  _et al._ (2019); Borsanyi  _et al._ (2020).
Moreover, the results for the phase boundary at finite chemical potential in
Fu  _et al._ (2020a); Gao and Pawlowski (2020b, a); Fischer (2019) reveal that
the phase boundary is still described well by the leading order expansion with
$\mu_{B}^{2}$-terms. We estimate, that this prediction is quantitatively
reliable within $\mu_{B}/T\lesssim 4$, using results from Fu  _et al._
(2020a); Gao and Pawlowski (2020b, a); Fischer (2019); Braun  _et al._
(2020a). This covers the regime studied in the present work.

Applying the two scale-matching relations in Equation 4 with the coefficients
Equation 5 and Equation 7 to the 2- and 2+1-flavour data of the QCD phase
boundary in Fu  _et al._ (2020a) leads us to Figure 2. In conclusion, this
impressive agreement provides non-trivial support for the scale-matching
procedure in QCD.

####  II.2.2 2- to 2+1-flavour scale-matching in LEFTs

The convincing quantitative accuracy of the linear scale-matching analysis
presented for QCD in the last section also sustains its use in the LEFT within
the present work. Note however, that we cannot simply take over the above QCD-
relations for the present LEFT, which lacks the backcoupling of the glue-
dynamics on both large temperature and chemical potential physics. Still, the
dominance of the leading order term $-\kappa\mu_{B}^{2}/T_{c}^{2}$ in the
model reflects the same property in QCD. This allows us to employ a respective
linear scale-matching for $\mu_{B}/T\lesssim 4$ as studied in the present
work.

Figure 3: Matching of the temperature-scale in the QCD-assisted 2-flavour LEFT
with $R_{42}^{B}$ in Equation 14, using the 2+1-flavour lattice results of
Borsanyi  _et al._ (2018). This leads to $c_{{}_{T}}=1.247(12)$ in Equation 5.
The $T/T_{c}$-dependence of $R_{42}^{B}$ is a prediction of the QCD-assisted
LEFT, and agrees quantitatively with the lattice results.

Analogously to QCD we choose the chiral crossover temperature at vanishing
chemical potential, $(T,\mu_{B})=(T_{c},0)$ for fixing the scale factor
$c_{{}_{T}}$. Moreover, in the present work we are interested in fluctuations
of conserved charges. Hence, instead of the renormalised condensate we use the
kurtosis of baryon number fluctuations, or rather
$R_{42}^{B}=\chi_{4}^{B}/\chi_{2}^{B}$, for the definition see Equation 13 and
Equation 14 with Equation 16, Equation 18. This leads us to the following
determination of $c_{T}$: While the temperature-dependence of $R^{B}_{42}$ is
a prediction of the LEFT, its absolute temperature has to be adjusted. This is
done by minimising the $\chi^{2}$ of the difference between the lattice result
and the LEFT-prediction as a function of the rescaled absolute temperature
$c_{T}T_{c}$, leading us to

$\displaystyle c_{{}_{T}}=1.247(12)\,,$ (8)

The respective result for $R^{B}_{42}$ is shown in Figure 3 in comparison to
the lattice result from Borsanyi  _et al._ (2018). The two curves match
quantitatively supporting the predictive power of the LEFT.

For the scale-matching of $\mu_{B}$ with the curvature
$-\kappa\mu_{B}^{2}/T_{c}^{2}$ we have a plethora of results from state of the
art functional approaches: $\kappa=0.0142(2)$ in Fu  _et al._ (2020a),
$\kappa=0.0150(7)$ in Gao and Pawlowski (2020a) and $\kappa=0.0147(5)$ in Gao
and Pawlowski (2020b), the very recent update of Gao and Pawlowski (2020a).
Lattice results are provided with $\kappa=0.015(4)$ in Bazavov  _et al._
(2019), $\kappa=0.0149(21)$ in Bellwied  _et al._ (2015), $\kappa=0.0153(18)$
in Borsanyi  _et al._ (2020). Both, functional and lattice results agree
within the respective (statistical and systematic) errors with $\kappa\approx
0.015$.

Having adjusted the temperature with results from the WB-collaboration,
Bellwied  _et al._ (2015), we use $\kappa=0.0153(18)$ from Borsanyi  _et al._
(2020) for internal consistency. Note that the results presented here do only
change marginally if using one $\kappa$ in the range $\kappa=(0.0142-0.0153)$.
Within the current LEFT we obtain $\kappa_{{}_{\textrm{LEFT}}}=0.0193$. In
comparison, $\kappa_{{}_{\textrm{LEFT}}}$ is larger than the 2-flavour QCD
result in Fu  _et al._ (2020a) with $\kappa=0.0179(8)$. This reflects the lack
of glue-dynamics in the LEFT. We use this in the relation Equation 7 instead
of $\kappa^{(N_{f}=2)}$, and arrive at

$\displaystyle
c_{\mu_{B}}=c_{{}_{T}}\left(\frac{\kappa^{N_{f}=(2+1)}}{\kappa_{{}_{\textrm{LEFT}}}}\right)^{1/2}=1.110(66)\,,$
(9)

with the LEFT-$c_{{}_{T}}$ from Equation 8.

In summary, as our first step towards a quantitative prediction for hyper-
order baryon number fluctuations, in this work we do not deal with the strange
quark as a dynamical degree of freedom for the moment, but rather take into
account its effect on the modification of the momentum-scale running via an
appropriate scale-matching as shown in Equation 4. The validity of the scale-
matching relations between 2- and 2+1-flavour QCD has been well verified in
this section by means of results from the first-principle functional QCD in Fu
_et al._ (2020a). These relations were applied to the present QCD-assisted
LEFT. The scale-matching was done with two observables relevant for the
fluctuation physics studied here: $R^{B}_{42}$ as a function of $T$ and the
curvature of the phase boundary $\kappa$, both at vanishing chemical
potential. This led us to the coefficients Equation 8 and Equation 9 in
Equation 4. The errors in these coefficients determine the errors of our
results in Section IV.

##  III Thermodynamics and Hyper-order baryon number fluctuations

The thermodynamic potential in the LEFT at finite temperature and baryon
chemical potential is readily obtained from the effective action in Equation
1, or rather from its integrated flow: we evaluate the effective action on the
solution of the quantum equations of motion (EoMs). In the present work we
consider only homogeneous (constant) solutions,
$(\sigma_{\textrm{EoM}},A_{0,\textrm{EoM}})$ with

$\displaystyle\frac{\partial V(\rho,L,\bar{L})}{\partial\sigma}=\frac{\partial
V(\rho,L,\bar{L})}{\partial L}=\frac{\partial
V(\rho,L,\bar{L})}{\partial\bar{L}}=0\,,$ (10)

while the quark fields vanish on the EoMs, $q,\bar{q}=0$. We also note that
the assumption of homogeneous solutions has to be taken with a grain of salt
for larger chemical potentials with $\mu_{B}/T\gtrsim 4$, see Fu  _et al._
(2020a). Such a scenario has been investigated in LEFTs, see e.g. the review
Buballa and Carignano (2015) and references therein.

With these preparations we are led to the grand potential
$\Omega[T,\mu_{B}]=V_{k=0}(\rho,L,\bar{L})$, the effective potential,
evaluated at vanishing cutoff scale $k=0$. It reads

$\displaystyle\Omega[T,\mu_{B}]=$ $\displaystyle
V_{\mathrm{glue}}(L,\bar{L})+V_{\mathrm{mat}}(\rho,L,\bar{L})-c\sigma\,,$ (11)

where the gluonic background field $A_{0}$ in Equation 2 has been reformulated
in terms of the Polyakov loop $L$ and its complex conjugate $\bar{L}$. As
mentioned before, the matter sector of the effective potential is integrated
out towards the IR limit $k=0$, for details see Appendix B. In turn, the glue
sector is independent of $k$, see Appendix D. The pressure of the system
follows directly from the thermodynamic potential,

$\displaystyle p=$ $\displaystyle-\Omega[T,\mu_{B}]\,.$ (12)

The generalised susceptibilities of the baryon number $\chi^{B}_{n}$ are
defined through the $n$-th order derivatives of the pressure w.r.t. the baryon
chemical potential, to wit,

$\displaystyle\chi_{n}^{B}$
$\displaystyle=\frac{\partial^{n}}{\partial(\mu_{B}/T)^{n}}\frac{p}{T^{4}}\,.$
(13)

To remove the explicit volume dependence, it is advantageous to consider the
ratio between the $n$\- and $m$-th order susceptibilities, defined by,

$\displaystyle R_{nm}^{B}$
$\displaystyle=\frac{\chi_{n}^{B}}{\chi_{m}^{B}}\,.$ (14)

The generalised susceptibilities are related to various cumulants of the
baryon number distribution, which can be measured in heavy-ion collision
experiments through the cumulants of its proxy, i.e., the net proton
distribution, see, e.g., Luo and Xu (2017) for details. For the lowest four
orders we get,

$\displaystyle\chi^{B}_{1}=$ $\displaystyle\frac{1}{VT^{3}}\braket{N_{B}}\,,$
(15) $\displaystyle\chi^{B}_{2}=$
$\displaystyle\frac{1}{VT^{3}}\braket{(\delta N_{B})^{2}}\,,$ (16)
$\displaystyle\chi^{B}_{3}=$ $\displaystyle\frac{1}{VT^{3}}\braket{(\delta
N_{B})^{3}}\,,$ (17) $\displaystyle\chi^{B}_{4}=$
$\displaystyle\frac{1}{VT^{3}}\Big{(}\braket{(\delta
N_{B})^{4}}-3\braket{(\delta N_{B})^{2}}^{2}\Big{)}\,,$ (18)

with $\braket{\cdots}$ denoting the ensemble average and $\delta
N_{B}=N_{B}-\braket{N_{B}}$. Thus the mean value of the net baryon number of
the system is given by $M=VT^{3}\chi_{1}^{B}$, the variance
$\sigma^{2}=VT^{3}\chi_{2}^{B}$, skewness
$S=\chi_{3}^{B}/(\chi_{2}^{B}\sigma)$, and the kurtosis
$\kappa=\chi_{4}^{B}/(\chi_{2}^{B}\sigma^{2})$, respectively.

Figure 4: $R^{B}_{42}=\chi^{B}_{4}/\chi^{B}_{2}$ (left panel),
$R^{B}_{62}=\chi^{B}_{6}/\chi^{B}_{2}$ (middle panel), and
$R^{B}_{82}=\chi^{B}_{8}/\chi^{B}_{2}$ (right panel) as functions of the
temperature at vanishing baryon chemical potential ($\mu_{B}=0$). Results from
the QCD-assisted LEFT are compared with lattice results from the HotQCD
collaboration Bazavov  _et al._ (2017a, b, 2020) and the Wuppertal-Budapest
collaboration (WB) Borsanyi  _et al._ (2018). The inset in the plot of
$R^{B}_{82}$ shows its zoomed-out view. Our results agree quantitatively with
the WB-results , and are qualitatively compatible with the HotQCD results. We
also compare to the predictions of a hadron resonance gas (HRG) Braun-
Munzinger  _et al._ (2004a), which predicts only a very mild increase of
$R_{n2}^{B}$ from unity with increasing $T$.  Figure 5:
$R^{B}_{10\,2}=\chi^{B}_{10}/\chi^{B}_{2}$ as a function of the temperature
with $\mu_{B}=0$ from the QCD-assisted LEFT.

In this work the emphasis is, however, put on the baryon number fluctuations
of orders higher than the fourth, i.e. $\chi_{n>4}^{B}$, which are named
hyper-order baryon number fluctuations. As the low-order ones, the hyper-order
susceptibilities are also connected to their respective cumulants, and their
relations, taking the fifth through eighth ones for instance, are given as
follows,

$\displaystyle\chi^{B}_{5}=$
$\displaystyle\frac{1}{VT^{3}}\Big{(}\braket{(\delta
N_{B})^{5}}-10\braket{(\delta N_{B})^{2}}\braket{(\delta
N_{B})^{3}}\Big{)}\,,$ (19) $\displaystyle\chi^{B}_{6}=$
$\displaystyle\frac{1}{VT^{3}}\Big{(}\braket{(\delta
N_{B})^{6}}-15\braket{(\delta N_{B})^{4}}\braket{(\delta N_{B})^{2}}$
$\displaystyle-10\braket{(\delta N_{B})^{3}}^{2}+30\braket{(\delta
N_{B})^{2}}^{3}\Big{)}\,,$ (20) $\displaystyle\chi^{B}_{7}=$
$\displaystyle\frac{1}{VT^{3}}\Big{(}\braket{(\delta
N_{B})^{7}}-21\braket{(\delta N_{B})^{5}}\braket{(\delta N_{B})^{2}}$
$\displaystyle-35\braket{(\delta N_{B})^{4}}\braket{(\delta N_{B})^{3}}$
$\displaystyle+210\braket{(\delta N_{B})^{3}}\braket{(\delta
N_{B})^{2}}^{2}\Big{)}\,,$ (21) $\displaystyle\chi^{B}_{8}=$
$\displaystyle\frac{1}{VT^{3}}\Big{(}\braket{(\delta
N_{B})^{8}}-28\braket{(\delta N_{B})^{6}}\braket{(\delta N_{B})^{2}}$
$\displaystyle-56\braket{(\delta N_{B})^{5}}\braket{(\delta
N_{B})^{3}}-35\braket{(\delta N_{B})^{4}}^{2}$
$\displaystyle+420\braket{(\delta N_{B})^{4}}\braket{(\delta N_{B})^{2}}^{2}$
$\displaystyle+560\braket{(\delta N_{B})^{3}}^{2}\braket{(\delta
N_{B})^{2}}-630\braket{(\delta N_{B})^{2}}^{4}\Big{)}\,.$ (22)

Different aspects of hyper-order fluctuations have been studied in mean-field
approximations in the past, see e.g. Wagner  _et al._ (2010); Karsch  _et al._
(2011); Schaefer and Wagner (2012). However, due to the decisive role of non-
perturbative quantum fluctuations for these quantities, a treatment beyond
mean-field, as in the present work, is necessary for their accurate
description. A first study in this direction, discussing hyper-order
fluctuations up to $\chi_{8}$ within a PQM model with the fRG at small
$\mu_{B}/T$ can be found in Friman  _et al._ (2011b).

Figure 6: $R^{B}_{42}$ (left panel), $R^{B}_{62}$ (middle panel), and
$R^{B}_{82}$ (right panel) as functions of the temperature at several values
of $\mu_{B}$. Insets in each plot show their respective zoomed-out view.

##  IV Numerical results and discussions

In this section we present and discuss our numerical results for hyper-order
fluctuations on the freeze-out curve. At vanishing chemical potential the
lower orders are compared to results from lattice calculations. We then
discuss the implications of our predictions for the hyper-order baryon number
fluctuations for decreasing collision energies (increasing chemical potential)
for heavy-ion collision experiments.

###  IV.1 Hyper-order baryon number fluctuations at vanishing density:
benchmarks and predictions

We start our discussion of the numerical results in our QCD-assisted low-
energy effective theory with benchmark results at vanishing chemical
potential, $\mu_{B}=0$. We have already seen in Section II.2 that the fourth
order fluctuations $R^{B}_{42}$, Equation 14, agrees quantitatively with the
respective lattice result, see Figure 3 and Figure 4, left panel. We emphasise
again that the thermal dependence of $R^{B}_{42}$ is a prediction of the
present LEFT. Now we also compare the temperature dependence of the hyper-
fluctuations $R^{B}_{62}$ and $R^{B}_{82}$ with the corresponding lattice
results in the middle and right panels of Figure 4, respectively. We depict
both our numerical results and lattice results from the HotQCD collaboration,
Bazavov  _et al._ (2017a, b, 2020), and the Wuppertal-Budapest collaboration,
Borsanyi  _et al._ (2018). Note that lattice results in Figure 4 by the
Wuppertal-Budapest collaboration in Borsanyi  _et al._ (2018), and
$R^{B}_{62}$ and $R^{B}_{82}$ by the HotQCD collaboration in Bazavov  _et al._
(2017a) are not continuum extrapolated.

With the increase of the order of fluctuations, the uncertainties of lattice
results increase significantly. Moreover, the eighth-order fluctuations
$R^{B}_{82}$ obtained by the two collaborations show a significant
quantitative difference, although their form is qualitatively consistent with
each other.

The hyper-order baryon number fluctuations computed in the current setup are
in qualitative agreement with both lattice results. However, our results
single out the lattice results of the Wuppertal-Budapest collaboration, with
which we observe quantitative agreement. This situation is very reminiscent of
the pressure prediction in Herbst  _et al._ (2014): similarly to the current
situation with lattice predictions of hyper-fluctuations, the pressure
predictions of the lattice collaborations had not converged yet. A less
advanced version of the current QCD-assisted LEFT framework then predicted the
correct pressure result. We have also computed the hyper-order fluctuations
within a simple hadron resonance gas model Braun-Munzinger  _et al._ (2004a).
Essentially, they are all constant with only a very minor monotonic increase
with $T$ for $T\lesssim 140$ MeV, starting from unity at $T=0$. This is in
quantitative agreement with our findings at low temperatures. In summary, the
current setup passes all benchmark tests quantitatively and provides the full
temperature-dependence of hyper-fluctuations. We have also computed even
higher order baryon number fluctuations. In Figure 5 we show our result for
the temperature-dependence of the tenth order ratio
$R^{B}_{10\,2}=\chi^{B}_{10}/\chi^{B}_{2}$ at vanishing chemical potential,
$\mu_{B}=0$. So far no lattice results for the tenth-order fluctuation are
available, and the dependence of $R^{B}_{10\,2}$ on the temperature in Figure
5 is a prediction by the current QCD-assisted LEFT and awaits confirmation by
other calculations, e.g., lattice QCD, in the future.

Figure 7: $R^{B}_{42}$ (left panel) and $R^{B}_{62}$ (right panel) as
functions of $\mu_{B}/T$ with $T=155$ MeV and $T=160$ MeV from the eighth-
order Taylor expansion in $(\mu_{B}/T)^{2}$ around vanishing $\mu_{B}$ shown
in Equation 24. Our results from the QCD-assisted LEFT are compared to those
from lattice QCD by the HotQCD collaboration Bazavov  _et al._ (2020) and the
Wuppertal-Budapest collaboration Borsanyi  _et al._ (2018). We show the
comparison to HotQCD in the inlays, as these results deviate considerably from
both ours and the WB results.

###  IV.2 Hyper-order baryon number fluctuations at finite density

With successfully passing the benchmark tests, we proceed with the results for
baryon number fluctuations at finite chemical potential. This will allow us to
finally compare the theoretical predictions on the freeze-out curve with the
experimental measurements in Section IV.4.

Equally relevant is the self-consistent evaluation of the reliability of a
Taylor expansion in baryon-chemical potential that underlies the extension of
lattice results at vanishing chemical potential to $\mu_{B}\neq 0$. This is
particularly important for predictions of the location of the critical end
point based on such an expansion. Here we can investigate the reliability
range of the Taylor expansion around $\mu_{B}=0$ by comparison to our direct
computation at finite $\mu_{B}$.

First we investigate the temperature-dependence of the baryon number
fluctuations for different chemical potentials. This also allows us to discuss
the reliability bounds of the current LEFT-setup for increasing chemical
potential. In Figure 6 we show the temperature-dependence of the ratios
$R^{B}_{42}$, $R^{B}_{62}$ and $R^{B}_{82}$ for chemical potentials
$\mu_{B}=0,100,160,200,300,400$ MeV. First, we note that at small temperatures
the thermodynamic properties of the QCD medium are well described by a dilute
gas of hadrons, where the net-baryon number follows a Skellam distribution.
Thus, all ratios approach unity at vanishing temperature. At very large
temperature the system is governed by asymptotically free quarks, where
fluctuations approach to the trivial Stefan-Boltzmann limit, and $R^{B}_{n2}$
goes to zero for all $n>4$ at large $T$. Consequently, the non-trivial
behaviour of the fluctuations between these two limiting cases shown in Figure
6 is directly related to the crossover from the hadronic- to the quark-gluon
regime of QCD. The magnitude, but also the error of the fluctuations grow with
increasing chemical potential. Both effects are more pronounced for higher
order fluctuations. The increase in magnitude is directly linked to the
sharpening of the chiral crossover with increasing chemical potential, cf.
Figure 2. We expect that the current LEFT-setup is gradually loosing its
predictive power for fluctuations on the freeze-out curve due to the rapid
increase of the computational error for higher-order fluctuations at large
baryon chemical potential, e.g., $R^{B}_{82}$ with $\mu_{B}\gtrsim 200$ MeV.
All results of the subsequent investigations have to be evaluated with this
estimate on our systematic error.

Figure 8: Comparison between the direct calculation of baryon number
fluctuations $R^{B}_{42}$ (upper panels) and $R^{B}_{62}$ (lower panels) via
Equation 13 at finite $\mu_{B}$ and the Taylor expansion up to
$\chi^{B}_{8}(0)$ in Equation 24 (left panels) and to $\chi^{B}_{10}(0)$
(right panels). Both calculations are performed within the QCD-assisted LEFT
used in the present work. $R^{B}_{42}$, $R^{B}_{62}$ are shown as functions of
$\mu_{B}/T$ with $T=155$ MeV and $T=160$ MeV.

For the evaluation of the reliability regime of the Taylor expansion about
vanishing chemical potential we consider the Taylor expansion of the pressure
in Equation 12 in powers of $\hat{\mu}_{B}\equiv\mu_{B}/T$ around
$\hat{\mu}_{B}=0$. This leads us to

$\displaystyle\frac{p(\mu_{B})}{T^{4}}$
$\displaystyle=\frac{p(0)}{T^{4}}+\sum_{i=1}^{\infty}\frac{\chi^{B}_{2i}(0)}{(2i)!}\,\hat{\mu}_{B}^{2i}\,,$
(23)

with the expansion coefficients $\chi^{B}_{2i}(0)=\chi^{B}_{2i}(\mu_{B}=0)$,
the hyper-order fluctuations of the baryon charge. In Equation 23 we have
suppressed the temperature-dependence of all functions for the sake of
readability. Truncating the Taylor expansion in Equation 23 at the eighth
order, $\hat{\mu}_{B}^{8}$, and employing Equation 13, we obtain the expanded
baryon number fluctuations,

$\displaystyle\chi^{B}_{2}(\mu_{B})\simeq$
$\displaystyle\chi^{B}_{2}(0)+\frac{\chi^{B}_{4}(0)}{2!}\hat{\mu}_{B}^{2}+\frac{\chi^{B}_{6}(0)}{4!}\hat{\mu}_{B}^{4}+\frac{\chi^{B}_{8}(0)}{6!}\hat{\mu}_{B}^{6}\,,$
$\displaystyle\chi^{B}_{4}(\mu_{B})\simeq$
$\displaystyle\chi^{B}_{4}(0)+\frac{\chi^{B}_{6}(0)}{2!}\hat{\mu}_{B}^{2}+\frac{\chi^{B}_{8}(0)}{4!}\hat{\mu}_{B}^{4}\,,$
$\displaystyle\chi^{B}_{6}(\mu_{B})\simeq$
$\displaystyle\chi^{B}_{6}(0)+\frac{\chi^{B}_{8}(0)}{2!}\hat{\mu}_{B}^{2}\,.$
(24)

In Figure 7 we show the ratios $R^{B}_{42}=\chi^{B}_{4}/\chi^{B}_{2}$ and
$R^{B}_{62}=\chi^{B}_{6}/\chi^{B}_{2}$ based on the Taylor expansion for two
fixed temperatures: $T=155$ MeV (close to the crossover temperature $T_{c}$ at
$\mu_{B}=0$) and $T=160$ MeV (slightly above $T_{c}$). As an input we use
$\chi^{B}_{2i}(0)$ ($i\\!=\\!1,\,2,\,3,\,4$) from the current setup as well as
from the lattice (HotQCD collaboration Bazavov  _et al._ (2020) and Wuppertal-
Budapest collaboration Borsanyi  _et al._ (2018)), depicted already in Figure
4. As expected, the LEFT-results for the $\mu_{B}$-dependence of $R^{B}_{42}$
and $R^{B}_{62}$ agrees qualitatively with both lattice results. Moreover, it
agrees quantitatively with the Wuppertal-Budapest result. Note that
constraints, e.g., strangeness neutrality or a fixed ratio of the electric
charge to the baryon number density, are not implemented in all the results in
Figure 7. For more details about effects of these constraints on the
fluctuations and correlations of conserved charges, see the relevant
discussions in, e.g., Bazavov  _et al._ (2017a, 2020) in lattice QCD and Fu
_et al._ (2020b, 2019) in fRG.

Since we are not restricted by a sign problem within the fRG approach, the
$\chi^{B}_{n}(\mu_{B})$’s in Equation 13 can also be computed directly for the
current QCD-assisted LEFT without resorting to a Taylor expansion. By
comparing this to the results of the Taylor expansion, we can study its range
of validity. The results are presented in the left panel of Figure 8, again
for $T=155$ MeV and $T=160$ MeV and with the Taylor expansion up to eights
order in $\hat{\mu}_{B}$.

We observe that the result for $R^{B}_{42}$ from the Taylor expansion in
Equation 24 agrees quantitatively with that from the full calculation for
$\mu_{B}/T\lesssim 1.2$ for $T=155$ MeV and $\mu_{B}/T\lesssim 1.5$ for
$T=160$ MeV. Not surprisingly, this reliability regime is reduced
significantly for the hyper-order fluctuation $R^{B}_{62}$ to
$\mu_{B}/T\lesssim 1.2$ for $T=160$ MeV and $\mu_{B}/T\lesssim 0.8$ for
$T=155$ MeV. For larger $\mu_{B}/T$, the fluctuations show a non-trivial
oscillatory behaviour that cannot be captured by a (low-order) Taylor
expansion. We emphasise that this is not an artefact of our model, but rather
a generic, physical feature of these fluctuation observables. It reflects the
increasingly pronounced non-monotonic temperature dependence of $R^{B}_{n2}$
due to long-range correlations in the crossover region, as seen in Figure 6.
In particular, $R^{B}_{n2}$ develops distinctive areas around the crossover at
larger $\mu_{B}$ where its value varies significantly, even including sign
changes. By following a line of fixed $T$ close to $T_{c}$ and increasing
$\mu_{B}$ in the phase diagram, these areas are crossed eventually, leading to
the oscillatory behaviour seen in Figure 8. This is also evident in the right
plot of Figure 11, where we show the magnitude of $R^{B}_{42}$ in the phase
diagram. Since these strong fluctuations only occur at larger $\mu_{B}$, the
resulting characteristic qualitative features cannot be captured by a (low-
order) Taylor expansion at $\mu_{B}=0$; it is bound to fail at the onset of
the oscillatory behaviour, i.e. around $\mu_{B}/T\gtrsim 1$.

It is also interesting to evaluate to what extend a higher-order Taylor
expansion can improve its reliability. We therefore include our prediction for
$R^{B}_{10\,2}$ from Figure 5, hence extending the Taylor expansion in
Equation 24 to the tenth order. The resulting comparison is shown in right
panel of Figure 8. Interestingly, this does not change the compatibility
regime for $T=160$ MeV significantly. In turn, there are significant changes
for $T=155$ MeV. While the deviations for $R_{42}^{B}$ start to grow at
roughly the same $\mu_{B}/T$ as for the eights-order expansion, that is
$\mu_{B}/T\lesssim 1.2$, the result for $R_{42}^{B}$ stays compatible with the
full result for larger values. For $R_{62}^{B}$ the reliability regime is
essentially doubled: it rises from $\mu_{B}/T\lesssim 0.8$ to
$\mu_{B}/T\lesssim 1.5$.

The analysis above suggests the following picture: We have a temperature-
dependent reliability range of the Taylor expansion,

$\displaystyle T=155\,\textrm{MeV}:\quad$
$\displaystyle[\mu_{B}/T]_{\textrm{Max}}\approx 1.5\,,$ $\displaystyle
T=160\,\textrm{MeV}:\quad$ $\displaystyle[\mu_{B}/T]_{\textrm{Max}}\approx
1.2\,.$ (25)

Moreover, the results have already converged for $R_{42}^{B},R_{62}^{B}$ for
$T=160$ MeV as well as for $R^{B}_{42}$ for $T=155$ MeV within the eighth
order and for $\mu_{B}/T\lesssim[\mu_{B}/T]_{\textrm{Max}}(T)$. In turn,
convergence for $R^{B}_{62}$ for $\mu_{B}/T\lesssim[\mu_{B}/T]_{\textrm{Max}}$
requires the tenth order Taylor expansion for $T=155$ MeV. Note that the
values of $[\mu_{B}/T]_{\textrm{Max}}$ in Equation 25 might be increased a bit
when terms of orders higher than the tenth one are included in the Taylor
expansion. However, as we have discussed above, the full results of the baryon
number fluctuations show a non-trivial oscillatory behaviour, which is
generic, and stems from the fact that the system crosses over different phases
with the increase of $\mu_{B}$, cf. the right panel of Figure 11. For a
discussion of the radius of convergence based on mean-field theory we refer to
Karsch  _et al._ (2011). In general, it is given by the distance to the
nearest singularity of the equation of state in the complex $\mu_{B}/T$ plane.
Hence, possible candidates for this singularity are the CEP, the Roberge–Weiss
endpoint at imaginary $\mu_{B}$ Roberge and Weiss (1986), or the Yang–Lee edge
singularity in the complex plane Yang and Lee (1952). Evidently, the distance
to the CEP is far larger than the radius estimated here. In turn, the closest
endpoint at imaginary chemical potential is at $|\mu_{B}/T|\leq\pi$. For
physical quark masses it is most probably close to $|\mu_{B}/T|=\pi$. For a
discussion within QCD-flows see Braun  _et al._ (2011a), for lattice results
see e.g. Philipsen and Pinke (2014). Hence, a particularly intriguing option
is the Yang–Lee edge singularity, for a discussion see e.g. Stephanov (2006b);
Mukherjee and Skokov (2019). However, while the location of the edge
singularity has been determined for critical $O(N)$ theories Connelly  _et
al._ (2020), it is still unknown for QCD.

This interpretation also implies that the results from the Taylor expansion
fail to agree even qualitatively with the correct $\mu_{B}/T$-dependence for
$\mu_{B}/T_{c}\gtrsim[\mu_{B}/T]_{\textrm{Max}}(T)$, see Figure 8.
Interestingly, $[\mu_{B}/T]_{\textrm{Max}}(T)$ seems to grow for smaller
temperatures. Whether or not this holds true requires a more systematic study,
which will be considered elsewhere. In conclusion, the extrapolation of
fluctuations of conserved charges in the vicinity of the chiral crossover
within a Taylor expansion loses its predictive power for $\mu_{B}\gtrsim 200$
MeV, at least at tenth order.

Figure 9: $R^{B}_{31}$ (left panel), $R^{B}_{51}$ (middle panel), and
$R^{B}_{71}$ (right panel) as functions of the temperature at several values
of $\mu_{B}$. Insets in each plot show their respective zoomed-out view.

In Figure 9 we show our full results for the temperature-dependence of
$R^{B}_{31}$, $R^{B}_{51}$, and $R^{B}_{71}$ with different values of baryon
chemical potential. A further relevant odd fluctuation observable is
$R^{B}_{32}$, depicted in Figure 10. Its experimental analogue, the proton
number fluctuation $R^{p}_{32}$ has been already measured in Au+Au central
(0-5%) collisions at STAR, a comparison will be presented and discussed in
Section IV.5.

Figure 10: Baryon number fluctuation $R^{B}_{32}$ as a function of the
temperature at several values of $\mu_{B}$.

###  IV.3 Determination of the freeze-out curve

Figure 11: Left panel: chemical freeze-out temperature and baryon chemical
potential in the $T\\!-\\!\mu_{B}$ plane. The blue pentagons and red circles
show the freeze-out data from Andronic et al. Andronic  _et al._ (2018) and
STAR experiment Adamczyk  _et al._ (2017), respectively. The blue dashed line
represents the parametrisation of blue pentagons through Equation 26a and
Equation 26b. The red solid and green dotted lines show the parametrisation of
the STAR data based on all the seven data points, and only the four data
points in the middle region ($100\,\mathrm{MeV}\lesssim\mu_{B}\lesssim
300\,\mathrm{MeV}$), respectively. The grey squares are obtained by
interpolating the blue pentagons. The inlay zooms in the low-$\mu_{B}$ region.
Right panel: Baryon number fluctuations $R^{B}_{42}$ in the $T\\!-\\!\mu_{B}$
plane. The freeze-out curve is the STAR Fit II. The dashed line at
$\mu_{B}/T=4$ constitutes the reliability bound of the computations in Fu  _et
al._ (2020a) based on the potential emergence of new degrees of freedom
discussed in Fu  _et al._ (2020a); Braun  _et al._ (2020a); Fischer (2019).
The dashed lines at $\mu_{B}/T=2,3$ are reliability estimates of lattice
results as well as old ones from functional approaches, see also Figure 2.

The quantitatively successful benchmark tests analysed in Section IV.1, and
the evaluation of baryon number fluctuations at finite chemical potential in
Section IV.2 allow us to discuss our main goal: the comparison of theoretical
predictions on the baryon number fluctuations with experimental measurements.

A direct comparison between theory and experiment is a very challenging task.
This is due to the fact that experimental data are affected by many factors.
First, this concerns the acceptance of the detector such as the transverse
momentum $p_{T}$ range, rapidity window and the centrality dependence, e.g.
Adamczyk  _et al._ (2014a); Luo (2015); Bzdak  _et al._ (2017); He and Luo
(2017); Adam  _et al._ (2020); Abdallah  _et al._ (2021), see Luo and Xu
(2017); Adamczyk  _et al._ (2017) for more details. Second, the physics setup
used in theory and experiment may differ by the presence of volume
fluctuations, e.g. Luo  _et al._ (2013); Chatterjee  _et al._ (2020a, b),
finite volume effects on the location of the chiral phase boundary, e.g. Braun
_et al._ (2011b, 2012); Tripolt  _et al._ (2014); Almasi  _et al._ (2017b);
Klein (2017); Li  _et al._ (2019b); Liu  _et al._ (2020); Wan  _et al._
(2020), the question of global baryon number conservation, e.g. He  _et al._
(2016); Braun-Munzinger  _et al._ (2017); Vovchenko  _et al._ (2020), the
inclusion of resonance decays, e.g. Nahrgang  _et al._ (2015); Zhang  _et al._
(2020), and others.

All these different effects and experimental restrictions give rise to non-
critical contributions to fluctuation observables in experiments, and pinning
down their contributions plays a pivotal role in identifying the critical
signals in the BES experiment. Additionally, due to critical slowing down,
non-equilibrium effects become important in the vicinity of the CEP Berdnikov
and Rajagopal (2000), which necessitates a theoretical description of the
dynamics of critical fluctuations. For more details about recent progress on
the dynamics of critical fluctuations in QCD, see Bluhm  _et al._ (2020a) and
references therein. We emphasise, however, that the present results of QCD-
assisted LEFT model are well outside the critical region. Therefore they are
not subject to critical scaling in the vicinity of the CEP.

In this work we will not take into account the non-critical and dynamical
effects discussed above. Instead, we assume that the measured cumulants of the
net-proton multiplicity distribution at a given collision energy are in one-
to-one correspondence to the calculated fluctuations in Equation 13 with
single values for $T$ and $\mu_{B}$ (with other collision parameters e.g., the
centrality and rapidity range fixed). Then, it is suggestive to identify the
values of $T$ and $\mu_{B}$ with the ones when the chemical freeze-out occurs,
viz. $T_{{}_{\textrm{CF}}}$ and ${\mu_{B}}_{{}_{\textrm{CF}}}$. Such an
approach for the comparison is usually employed in fluctuation studies of
equilibrium QCD matter within functional methods or lattice simulations, see
e.g. Fu and Pawlowski (2016); Fu  _et al._ (2016); Almasi  _et al._ (2017a);
Isserstedt  _et al._ (2019); Bazavov  _et al._ (2020).

We adopt the freeze-out temperatures and baryon chemical potentials from
Andronic  _et al._ (2018) and from the STAR experiment Adamczyk  _et al._
(2017), which are shown in the left panel of Figure 11 by the blue pentagons
and red circles, respectively. They are both obtained from the analysis of
hadron yields in the statistical hadron resonance gas model, see the
aforementioned references for more details. The freeze-out data in Andronic
_et al._ (2018) has also been parametrised as functions of the collision
energy as follows

$\displaystyle{\mu_{B}}_{{}_{\textrm{CF}}}$
$\displaystyle=\frac{a}{1+0.288\sqrt{s_{\mathrm{NN}}}},$ (26a) with $a=1307.5$
MeV, and $\displaystyle T_{{}_{\textrm{CF}}}$
$\displaystyle=\frac{T^{(0)}_{{}_{\textrm{CF}}}}{1+\exp\big{(}2.60-\ln(\sqrt{s_{\mathrm{NN}}})/0.45\big{)}},$
(26b)

with $T^{(0)}_{{}_{\textrm{CF}}}=158.4$ MeV. This parametrisation is depicted
with the blue dashed line in the left panel of Figure 11. We use the same
parametrisation functions in Equation 26 to fit the freeze-out data in STAR
experiment, i.e., the red circle points. For this fit we invoke two
procedures, called STAR Fit I & II in the following:

For the first one, STAR Fit I, we simply take all 7 data points. The
corresponding freeze-out curve is depicted by the red solid line in the left
panel of Figure 11.

Figure 12: QCD-assisted LEFT (fRG-LEFT): Baryon number fluctuations
$R^{B}_{42}$ (top), $R^{B}_{62}$ (middle), and $R^{B}_{82}$ (bottom) as
functions of the collision energy. Left panels: the freeze-out points are
those from Andronic et al. Andronic  _et al._ (2018) (grey) and the STAR
experiment Adamczyk  _et al._ (2017) (red). Right panels: The freeze-out
curve, STAR Fit II, is obtained from the freeze-out parameters of the STAR
experiment Adamczyk  _et al._ (2017). The theoretical error bands show a
highly correlated error, and should be interpreted as a family of curves with
the same qualitative behaviour as the central curve. For more explanations see
Section IV.3 with Figure 11.
STAR data: $R^{p}_{42}$ (top) is the kurtosis of the net-proton distributions
measured in Au+Au central (0-5%) collisions Adam  _et al._ (2020).
$R^{p}_{62}$ (middle) is the result on the six-order cumulant of the net-
proton distribution at $\sqrt{s_{\mathrm{NN}}}$=200 GeV, 54.4 GeV and 27 GeV
with centrality 0-10% Abdallah  _et al._ (2021).

For the construction of the second one, STAR Fit II, we shall argue that some
of the data points are potentially flawed, or rather await a physics
explanation, and should be dropped accordingly in a fit based on Equation 26.
Accordingly, we drop the first two data points at small chemical potential as
well as the last one at the largest available chemical potential $\mu_{B}\sim
400$. From general considerations we do not expect the freeze-out curve to
rise with increasing chemical potential. Moreover, the physically motivated
fit formula does not describe sign-changes of the curvature of the freeze-out
curve. For a respective discussion and possible explanation for the only
apparent rise see Bluhm  _et al._ (2020b). The last data point is also not
well-described by the fitting procedure described here. This may indicate the
onset of a regime with different physics/phases. In this case, Equation 26
would not be an appropriate fit function. It may also indicate the onset of a
regime of rapidly worsening systematics. In this case more points are needed
in this regime.

The freeze-out line of STAR Fit II is depicted by the green dotted line in the
left panel of Figure 11. In comparison to STAR Fit I, STAR Fit II is located
at slightly lower temperatures, which is more pronounced when $\mu_{B}\gtrsim
200$ MeV. In the right panel of Figure 11, we show the baryon number
fluctuation $R^{B}_{42}$ in the $T\\!-\\!\mu_{B}$ plane. It can be observed
that a narrow blue band, indicating the regime of negative $R^{B}_{42}$,
develops around the crossover starting at $\mu_{B}\sim 250$ MeV. The freeze-
out curve STAR Fit II is approaching towards the boundary of the blue region
firstly at small $\mu_{B}$, and then deviates a bit from it at large
$\mu_{B}$. We emphasise that the large chemical potential region, and in
particular asymptotically large $\mu_{B}\gtrsim 500$ MeV, is beyond of the
reliability bound of the current computation, $\mu_{B}/T=4$, see Figure 2. For
a detailed discussion see Section II.2.

###  IV.4 Hyper-order baryon number fluctuations on the freeze-out curve

The determination of the freeze-out curve completes our setup, which enables
us to compute hyper-order baryon number fluctuations $R_{nm}^{B}$ along the
freeze-out line within the QCD-assisted LEFT. These results are then used to
compare with the experimental measurements of cumulants $R_{nm}^{p}$ of the
net-proton distribution from STAR experiment.

Before we discuss the numerical results, we also emphasise once more, that it
follows from the analysis of Section IV.2, that the simple extrapolation with
the Taylor expansion about $\mu_{B}=0$ lacks predictive power for
$\mu_{B}\gtrsim 250$ MeV, that is $\sqrt{s_{\rm NN}}\lesssim 15$ GeV, see
Equation 25. Moreover, it even lacks predictive power for the qualitative
behaviour.

In the left panel of Figure 12 we show the $\sqrt{s_{\rm NN}}$\- or chemical
potential dependence of the baryon number fluctuations $R^{B}_{42}$,
$R^{B}_{62}$, and $R^{B}_{82}$ for the freeze-out lines from Andronic et al.
Andronic  _et al._ (2018) and STAR Fit I. The freeze-out line from Andronic et
al. is obtained from an interpolation of the freeze-out data, the grey squares
in Figure 11.

In the right panel of Figure 12 we show the same observables for the freeze-
out line of STAR Fit II. As discussed in Section IV.3, we have singled out the
results for this freeze-out curve as the best-informed computation.

In both panels of Figure 12 we also show the experimental measurement of
cumulants of the net-proton distributions in the beam energy scan experiments
from the STAR collaboration. The fourth-order fluctuations, $R^{p}_{42}$, of
the net-proton multiplicity distributions are measured in Au+Au collisions
with centrality 0-5%, transverse momentum range
$0.4<p_{T}\,(\mathrm{GeV}/c)\,<2.0$, and rapidity $|y|<0.5$, cf. Adam  _et
al._ (2020) for more details. Moreover, results for the sixth-order cumulant
of the net-proton distribution, $R^{p}_{62}$, are also presented in the middle
plot of Figure 12, which are obtained at three values of the collision energy,
i.e., $\sqrt{s_{NN}}$=200 GeV, 54.4 GeV and 27 GeV with centrality 0-10%
Abdallah  _et al._ (2021).

The theoretical results for the fourth-order fluctuations $R^{B}_{42}$ from
the present QCD-assisted LEFT for all freeze-out curves considered are
compatible with the respective experimental measurement of the
$\kappa\sigma^{2}$ of net-proton distributions in 0-5% central Au+Au
collisions. In particular, the theoretical results feature a non-monotonic
$\sqrt{s_{\rm NN}}$-dependence: $R^{B}_{42}$ first decreases with decreasing
beam energy and then increases. The details of this behaviour, in particular
how pronounced it is, is highly sensitive to the precise location of the
freeze-out. For example, the increase at small $\sqrt{s_{\rm NN}}$ is larger
for smaller freeze-out temperatures. Thus, the weak increase for STAR Fit 1
originates in the slightly larger freeze-out temperature of this freeze-out
fit. This shows that even small variations in the freeze-out temperature have
a substantial effect on the fluctuations in this regime. The underlying reason
is that the freeze-out happens in or close to the crossover region, where the
fluctuations vary significantly, see Figure 6. Importantly, this regime cannot
be accessed within the extrapolation of the Taylor expansion at least within
the current order.

This entails that extrapolations based on a Taylor expansion are bound to fail
to describe the data in this regime reliably. Consequently this calls for
qualitatively improved direct theoretical computations at small beam-energies.
This is work in progress and we hope to report on the respective results soon.

Our results for the hyper-order fluctuations $R^{B}_{62}$ and $R^{B}_{82}$ are
shown in the middle and bottom panel of Figure 12. For small chemical
potentials or large collision energies both fluctuation observables are
negative. Moreover, $R^{B}_{62}$ decreases with decreasing
$\sqrt{s_{\mathrm{NN}}}$, while $R^{B}_{82}$ increases. The occurrence of non-
monotonicities of $R^{B}_{62}$ and $R^{B}_{82}$ at lower beam energies cannot
be shown within the accuracy limits of the current study.

For $\sqrt{s_{\mathrm{NN}}}$=200 GeV, 54.4 GeV and 27 GeV we can compare our
results for $R^{B}_{62}$ to STAR data within 0-10% centrality Abdallah  _et
al._ (2021). One observes that our results are in agreement with the
experimental data within errors at $\sqrt{s_{\mathrm{NN}}}$=200 GeV and 54.4
GeV, though the central value of STAR data at $\sqrt{s_{\mathrm{NN}}}$=54.4
GeV is positive. Both the theory and experiment show negative values for the
sixth-order fluctuations at the collision energy $\sqrt{s_{\mathrm{NN}}}$=27
GeV.

Another interesting property of the current LEFT setting is that the non-
monotonic behaviour of our results for $R^{B}_{n2}$ at large $\mu_{B}$ in
Figure 12 does not arise from critical physics: in the LEFT used here, the CEP
is at significantly larger $\mu_{B}\gtrsim 700$ MeV. Moreover, it is well
established that the critical region is only very small. It is already small
within mean-field approximations of low-energy effective theories, and
additionally shrinks considerably if quantum, thermal and density equilibrium
fluctuations are taken into account, see Schaefer and Wambach (2007).
Moreover, this does not change if transport processes are taken into account,
see Bluhm  _et al._ (2019).

In the present LEFT the increasing trend at large $\mu_{B}$ region originates
from two effects: First, fluctuations are enhanced since the chiral crossover
becomes sharper with increasing $\mu_{B}$. This leads to a stronger non-
monotonic behaviour of $R^{B}_{42}$ as a function of $T$, see Figure 6. This
sharpening is also present in the vicinity of a CEP. Second, the freeze-out
temperature is shifted away from the pseudo-critical temperature towards small
beam-energies, thereby probing different regimes of the cumulants. However, it
should be noted that the uncertainty of our results increases significantly in
the low energy region. These uncertainties include an estimate for the
systematic error of the QCD-assisted LEFT-approach. This systematic error
stems from the uncertainty in the matching of the in-medium scales of the LEFT
and QCD, encoded in the coefficients $c_{{}_{T}}$ and $c_{\mu}$ in Equation 4.
Moreover, for larger chemical potential the current LEFT lacks the back-
reaction of the $\mu_{B}$-dependence of the glue-dynamics. While inherently
small, it might still play a r$\hat{\textrm{o}}$le. Furthermore, it is found
that $R^{B}_{42}$ can also be suppressed in the regime of low collision
energies due to the effect of global baryon number conservation, cf. Braun-
Munzinger  _et al._ (2017, 2021); Vovchenko  _et al._ (2020) which will be
included in our future work.

###  IV.5 Search for the CEP

Figure 13: Baryon number fluctuation $R^{B}_{32}$ as a function of the
collision energy in comparison to STAR-data for $R^{p}_{32}$ (0-5%) centrality
Adam  _et al._ (2020). Left panel: the freeze-out points are those from
Andronic et al. Andronic  _et al._ (2018) (grey) and the STAR experiment
Adamczyk  _et al._ (2017) (red). Right panel: The freeze-out curve, STAR Fit
II, is obtained from the freeze-out parameters of the STAR experiment Adamczyk
_et al._ (2017). The theoretical error bands show a highly correlated error,
and should be interpreted as a family of curves with the same qualitative
behaviour as the central curve. For more explanations see Section IV.3 with
Figure 11.  Figure 14: Baryon number fluctuations $R^{B}_{31}$ (top),
$R^{B}_{51}$ (middle), and $R^{B}_{71}$ (bottom) as functions of the collision
energy. The freeze-out curve, STAR Fit II, is obtained from the freeze-out
parameters of the STAR experiment Adamczyk  _et al._ (2017). The theoretical
error bands show a highly correlated error, and should be interpreted as a
family of curves with the same qualitative behaviour as the central curve. For
more explanations see Section IV.3 with Figure 11.

The analysis in the previous sections entails that in the present QCD-assisted
LEFT the non-monotonic behaviour of baryon-number fluctuations is triggered by
the sharpening of the chiral crossover. This is highly non-trivial, since it
is evident, e.g., from Figure 6 that neither a system only in the hadronic-
nor only in the QGP phase could produce the beam-energy dependencies shown in
Figure 12.

In conclusion, the agreement of $R^{B}_{42}$, computed in the QCD-assisted
LEFT and the measured beam-energy dependence of $R^{p}_{42}$ shows, that the
latter could be a signature for the presence of a sharpening crossover between
these two phases. Whether or not it also signals the onset of the critical
region will be subject of a future improved study. In the context of this
latter study we also emphasise that the non-universal properties of the LEFT
such as the existence and location of the CEP may not quantitatively agree
with QCD as the latter regime lies outside the LEFT-regime with quantitative
reliability. Still, the present LEFT probably has the same qualitative non-
universal properties at large chemical potential, and it certainly has the
same universal ones.

A non-monotonic energy dependence for the fluctuations is a highly relevant
experimental observation, since this behaviour has been proposed as an
experimental signature of a CEP Stephanov  _et al._ (1999); Stephanov (2011).
The present analysis based on QCD-assisted LEFT model demonstrates that the
non-monotonic behaviour of fluctuations can serve as an indication of a CEP,
but is not necessarily a smoking gun signature for it. The latter requires the
extraction of critical scaling, or similar definite signatures such as the
detection of a first order regime for large $\mu_{B}$, etc. .

Still, the non-monotonic behaviour observed in both theory and experiments is
a clear signature for interesting strongly correlated physics, whose
uncovering requires joint and intensified effort of both, theory and
experiment. Of course, whether or not these properties carry over completely
to QCD remains to be seen.

Note also that the non-monotonic regime is far away from that covered by a
simple extrapolation of the Taylor expansion at $\mu_{B}=0$. It might be
covered by a resummation of the latter, which can already be investigated
within the present QCD-assisted LEFT. Constraints on such a resummation should
also make use of odd hyper-order fluctuations at finite chemical potential,
that are readily computed in the present setup:

A prominent and relevant example is $R_{32}$, already measured in the STAR
experiment. In Figure 13 we show our predictions for $R^{B}_{32}$ computed in
the current QCD-assisted LEFT on the different freeze-out curves defined in
Section IV.3, see in particular Figure 11. In this section it also has been
argued, that our best-informed freeze-out curve is given by STAR Fit II. The
respective results are shown in the right panel of Figure 13 in comparison
with the STAR data for $R^{p}_{32}$ (0-5% centrality). Indeed, these results
show the best compatibility with the experimental data. Moreover, within the
respective systematic and statistical errors the theoretical results with the
freeze-out curve STAR Fit II and the experimental data agree down to collision
energies of $\sqrt{s_{\mathrm{NN}}}\approx$ 14.5 GeV or $\mu_{B}\approx 250$
MeV.

Interestingly, below $\sqrt{s_{\mathrm{NN}}}\approx$ 14.5 GeV the experimental
data show a plateau, which is not present in the theoretical prediction. While
this is in the large chemical potential regime, in which the LEFT gradually
looses its predictive power, also the respective functional first principles
QCD computation in Fischer (2019); Fu  _et al._ (2020a); Gao and Pawlowski
(2020a, b), based on a grand potential, do not show any sign of new physics in
this regime. This suggests that for $\sqrt{s_{\mathrm{NN}}}\lesssim$ 14.5 GeV
at least one of the implicit assumptions underlying the identification of
$R^{B}_{nm}$ with $R^{p}_{nm}$ within a grand canonical ensemble with variable
baryon charge (density) for given beam energies breaks down. As discussed
before, this asks for a re-assessment of the identification of baryon and
proton number fluctuations, finite volume effects and finite volume
fluctuations, the determination of the freeze-out curve for smaller collision
energies, the evaluation of non-equilibrium effects such as transport, and
finally the use of the grand potential in the theory computations. While
highly relevant and interesting, this goes far beyond the scope of the present
work and we defer any further investigation to future work.

The above example of $R_{32}$ demonstrates very impressively, that the odd
(hyper-order) fluctuation observables encode highly relevant physics
information which may be difficult or even impossible to extract from the even
orders. As a first step in this direction, finally aiming at a resummation of
the $\mu_{B}$-expansion that allows us to go beyond the validity regime of the
Taylor expansion, we also have computed the fluctuation observables
$R_{31},R_{51},R_{71}$ on the freeze-out curve STAR Fit II in Figure 14. An
experimental confirmation of the respective predictions at least for the lower
orders would be highly desirable.

The discussion in this section leaves us with the highly exciting possibility
of unravelling the location and properties of a potential CEP within a
combined experiment-theory analysis: First principle QCD at finite density
should provide us with a prediction for the location of the CEP in terms of
hyper-order fluctuations, $\textrm{Loc}_{\textrm{\tiny{CEP}}}(R_{nm})$. This
would allow us to use the experimental data on hyper-order fluctuation
observables $R_{nm}^{p}$ as input. We emphasise that this prediction does not
necessitate the observation of critical behaviour in the $R_{nm}$, but
utilises the details of the non-monotonicity of the $R_{nm}$.

In summary, such an analysis does explicitly not rely on the universal
property of critical scaling measured in the $R_{nm}$. Indeed, it uses the
non-universal properties of the $R_{nm}$ to predict the non-universal location
of the CEP, and hence is far more robust. We hope to report on this matter in
the near future.

##  V Conclusions

In this work we have computed baryon number fluctuations up to tenth order
with a QCD-assisted low-energy effective theory. This LEFT incorporates
quantum, thermal and density fluctuations from momentum scales less than 700
MeV within the functional renormalisation group approach, and is embedded in
QCD, for details see Section II. The quantitative predictability has been
benchmarked with a comparison of baryon number fluctuations at $\mu_{B}=0$ up
to the eighth order from the lattice, see Section IV.1. Our results are in
quantitative agreement with that from the Wuppertal-Budapest collaboration,
and are compatible with that of the HotQCD collaboration, as shown in Figure
4.

Our direct computation at finite $\mu_{B}$, presented in Section IV.2, has
allowed us to assess the range of validity of the Taylor expansion of the free
energy of QCD around $\mu_{B}=0$. Such an expansion is commonly used to
extrapolate lattice results to finite density. We have shown that the
expansion up to tenth order in $\mu_{B}/T$ is only valid for
$\mu_{B}/T\lesssim 1.5$ in the chiral crossover regime, see Figure 8. Beyond
this range, the Taylor expansion, at least to this order, fails to even
capture the qualitative behaviour of the fourth- and sixth-order baryon number
fluctuations. Thus, results for fluctuations at the freeze-out curve based on
a Taylor expansion around $\mu_{B}=0$ should be interpreted with great caution
for $\mu_{B}/T\gtrsim 1.5$, as relevant physical effects might not be captured
by this extrapolation.

The main goal of the current work was the computation of baryon number
fluctuations and in particular hyper-order fluctuations along the freeze-out
curve at collision energies $\sqrt{s}\gtrsim 7.7$ GeV. The respective results
are discussed in Section IV.4. They have been compared to experimental data of
net-proton number cumulants from STAR for different estimates for the freeze-
out curve, see Figure 12. Our result for the kurtosis, $R^{B}_{42}$, is in
good agreement with the experimental data for collision energies
$\sqrt{s}\gtrsim 7.7$ GeV. In particular the increasing trend at lower beam
energies $\sqrt{s}\lesssim 19.6$ GeV is captured well. This non-monotonicity
is also present in the hyper-order fluctuations $R^{B}_{62},R^{B}_{82}$. We
also note that a comprehensive comparison for the higher order cumulants is
not possible due to the lack of experimental data. Accordingly, our results in
Figure 12 are predictions that await experimental verification.

We have also investigated the twofold origin of the non-monotonicity for
$\sqrt{s}\lesssim 19.6$ GeV in the present LEFT: First, for increasing
chemical potential the chiral crossover gets sharper. Secondly, for smaller
beam energies the freeze-out temperature may move away from the pseudo-
critical temperatures. In the current setup both phenomena happen far away
from a potential critical end point in the LEFT located at
$\sqrt{s}_{\textrm{CEP}}\lesssim 3$ GeV. The latter regime is also safely
outside the reliability regime of the current setup, which gradually looses
reliability for $\sqrt{s}\lesssim 27$ GeV. However, its qualitative features
may well be present in QCD. The current LEFT-results and its upgrades towards
first principle QCD can be compared with the future experimental results of
the high statistic data taken from the second phase of RHIC beam energy scan
(BES-II, 2019-2021). From the year of 2018 to 2020, the STAR experiment has
collected high statistics data of Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}$
= 9.2, 11.5, 14.6, 19.6 and 27 GeV in the collider mode, and
$\sqrt{s_{\mathrm{NN}}}$ = 3.0 – 7.7 GeV in the fixed target mode. These data
give us access to the QCD phase structure for baryon chemical potential up to
$\mu_{B}$ $\approx$ 720 MeV.

Related further steps in a comprehensive understanding of the physics of
fluctuations in a heavy ion collision has been undertaken in Section IV.5,
where we have presented results for odd fluctuations observables
$R^{B}_{31},R^{B}_{51},R^{B}_{71}$ as well as $R^{B}_{32}$. While the former
observables have not been measured yet, $R^{B}_{32}$ agrees quantitatively
within the systematic and statistical error with the STAR measurement for
$R^{p}_{32}$ for collision energies $\sqrt{s_{\mathrm{NN}}}\gtrsim$ 14.5 GeV.
For smaller energies, the experimental data show a plateau, that may indicate
the loss of one or several of the underlying assumption in the identification
of theoretical equilibrium computations of baryon number fluctuations
$R^{B}_{nm}$ in a grand canonical ensemble with the experimental results for
proton number fluctuations on the freeze-out curve.

In summary, the non-monotonicities of hyper-order fluctuations, observed both
in experiment and theory, are important signatures for interesting physics in
the border regime between quark-gluon plasma and the hadron phase. This of
course can include a potential CEP, and in any case deserves further
investigation from both experiment and theory. In particular, we envisage that
experimental data of fluctuation observables and their dependence on collision
energy allow us to constrain the onset regime of this strongly correlated
physics/CEP. Importantly, such a prediction does not rely on the observation
of critical scaling in the hyper-order fluctuations, but is far more robust,
for more details see Section IV.5. We hope to report on this in the near
future.

###### Acknowledgements.

We thank Jens Braun, Rob Pisarski, Bernd-Jochen Schaefer and Nu Xu for
discussions. The work was supported by the National Key Research and
Development Program of China (Grant No. 2020YFE0202002 and 2018YFE0205201) and
the National Natural Science Foundation of China under (Grant No. 11775041,
11828501, 11890711 and 11861131009). The work is also supported by EMMI, and
the BMBF grant 05P18VHFCA. It is part of and supported by the DFG
Collaborative Research Centre SFB 1225 (ISOQUANT) and the DFG under Germany’s
Excellence Strategy EXC - 2181/1 - 390900948 (the Heidelberg Excellence
Cluster STRUCTURES).

##  Appendix A The fRG-approach to QCD & LEFTs

The functional renormalisation group or flow equation for QCD provides the
evolution of its effective action $\Gamma_{k}$ with an infrared cutoff scale
$k$. Here we use the setup with dynamical hadronisation, Gies and Wetterich
(2002, 2004); Pawlowski (2007); Braun (2009); Floerchinger and Wetterich
(2009); Fu  _et al._ (2020a). The formulation used here has been developed in
Mitter  _et al._ (2015); Braun  _et al._ (2016); Rennecke (2015a); Cyrol  _et
al._ (2018a); Fu  _et al._ (2020a). Its current form has been described and
further developed in Fu  _et al._ (2020a), and for further details we refer to
this work. The flow equation of the QCD effective action reads,

$\displaystyle\partial_{t}\Gamma_{k}[\Phi]=$
$\displaystyle\frac{1}{2}\mathrm{Tr}\Big{(}G_{AA,k}\partial_{t}R_{A,k}\Big{)}-\mathrm{Tr}\Big{(}G_{c\bar{c},k}\partial_{t}R_{c,k}\Big{)}$
$\displaystyle-\mathrm{Tr}\Big{(}G_{q\bar{q},k}\partial_{t}R_{q,k}\Big{)}+\frac{1}{2}\mathrm{Tr}\Big{(}G_{\phi\phi,k}\partial_{t}R_{\phi,k}\Big{)}\,.$
(27)

In Equation 27, the $\Phi=(A,c,\bar{c},q,\bar{q},\phi)$ is a superfield that
comprises all fields. This also includes hadronic (composite) low energy
degrees of freedom introduced by dynamical hadronisation. The $G$’s and $R$’s
are the propagators and regulators of the different fields, respectively.
Diagrammatically it is depicted in Figure 1. For more works on QCD-flows at
finite temperature and density see Braun  _et al._ (2010); Braun (2009); Braun
_et al._ (2011a); Mitter  _et al._ (2015); Braun  _et al._ (2016); Rennecke
(2015a); Cyrol  _et al._ (2016, 2018a, 2018b); Fu  _et al._ (2020a); Braun
_et al._ (2020b, c), for reviews on QCD and LEFTs for QCD see Litim and
Pawlowski (1998); Berges  _et al._ (2002); Pawlowski (2007); Schaefer and
Wambach (2008); Gies (2012); Rosten (2012); Braun (2012); Pawlowski (2014);
Dupuis  _et al._ (2020).

For scales $k\lesssim 1$ GeV, the gluon decouples from the system due to its
confinement-related mass gap. For these momentum scales, the (off-shell)
dynamics of QCD is dominated by quarks and the emergent composite hadronic
degrees of freedom. In particular, the lowest lying meson multiplet, and
specifically the $\pi$ meson is driving the dynamics. The pion is the pseudo-
Goldstone boson of strong chiral symmetry breaking, and hence is the lightest
hadron with a mass $\sim 140$ MeV in the vacuum.

Consequently, in this regime with $k\lesssim 1$ GeV, the flow equation of the
QCD effective action in Equation 27 is reduced to,

$\displaystyle\partial_{t}\Gamma_{k}[\Phi]=$
$\displaystyle-\mathrm{Tr}\Big{(}G_{q\bar{q},k}\partial_{t}R_{q,k}\Big{)}+\frac{1}{2}\mathrm{Tr}\Big{(}G_{\phi\phi,k}\partial_{t}R_{\phi,k}\Big{)}\,,$
(28)

where $R_{q,k}$ and $R_{\phi,k}$ are the regulators for the quark and meson
fields, respectively. The full propagators in Equation 28 read,

$\displaystyle G_{q\bar{q}/\phi\phi,k}=$
$\displaystyle\left(\frac{1}{\Gamma^{(2)}_{k}[\Phi]+R_{k}}\right)_{q\bar{q}/\phi\phi}\,,$
(29)

with $\Gamma^{(2)}_{k}[\Phi]=\delta^{2}\Gamma_{k}[\Phi]/\delta\Phi^{2}$. In
this work we employ $3d$-flat or Litim regulators Litim (2000, 2001); Litim
and Pawlowski (2006),

$\displaystyle R_{\phi,k}(q_{0},\bm{q})$
$\displaystyle=Z_{\phi,k}\bm{q}^{2}r_{B}(\bm{q}^{2}/k^{2})\,,$ $\displaystyle
R_{q,k}(q_{0},\bm{q})$
$\displaystyle=Z_{q,k}i\bm{\gamma}\cdot\bm{q}r_{F}(\bm{q}^{2}/k^{2})\,,$ (30)

with

$\displaystyle r_{B}(x)$
$\displaystyle=\left(\frac{1}{x}-1\right)\Theta(1-x)\,,$ $\displaystyle
r_{F}(x)$ $\displaystyle=\left(\frac{1}{\sqrt{x}}-1\right)\Theta(1-x)\,,$ (31)

where $\Theta(x)$ denotes the Heaviside step function. Inserting the effective
action Equation 1 into the flow equation Equation 28, we arrive at

$\displaystyle\partial_{t}V_{\mathrm{mat},k}(\rho)=$
$\displaystyle\frac{k^{4}}{4\pi^{2}}\bigg{[}\big{(}N^{2}_{f}-1\big{)}l^{(B,4)}_{0}(\tilde{m}^{2}_{\pi,k},\eta_{\phi,k};T)$
$\displaystyle+l^{(B,4)}_{0}(\tilde{m}^{2}_{\sigma,k},\eta_{\phi,k};T)$
$\displaystyle-4N_{c}N_{f}l^{(F,4)}_{0}(\tilde{m}^{2}_{q,k},\eta_{q,k};T,\mu)\bigg{]}\,,$
(32)

where the threshold functions $l^{(B/F,4)}_{0}$ as well as other threshold
functions used in the following can be found in e.g., Fu  _et al._ (2020a);
Yin  _et al._ (2019). The dimensionless renormalised quark and meson masses
read,

$\displaystyle\tilde{m}^{2}_{q,k}=$
$\displaystyle\frac{h^{2}_{k}\rho}{2k^{2}Z^{2}_{q,k}}\,,\qquad\tilde{m}^{2}_{\pi,k}=\frac{V^{\prime}_{\mathrm{mat},k}(\rho)}{k^{2}Z_{\phi,k}}\,,$
(33) $\displaystyle\tilde{m}^{2}_{\sigma,k}=$
$\displaystyle\frac{V^{\prime}_{\mathrm{mat},k}(\rho)+2\rho
V^{\prime\prime}_{\mathrm{mat},k}(\rho)}{k^{2}Z_{\phi,k}}\,.$ (34)

The anomalous dimensions for the quark and meson fields in Equation 32 are
defined as

$\displaystyle\eta_{q,k}$
$\displaystyle=-\frac{\partial_{t}Z_{q,k}}{Z_{q,k}}\,,\qquad\eta_{\phi,k}=-\frac{\partial_{t}Z_{\phi,k}}{Z_{\phi,k}}\,,$
(35)

respectively. Accordingly, the flow equation for the mesonic anomalous
dimension is obtained from the (spatial) momentum derivative w.r.t.
$\bm{p}^{2}$ of the pion two-point function, to wit,

$\displaystyle\eta_{\phi,k}$
$\displaystyle=-\frac{1}{3Z_{\phi,k}}\delta_{ij}\frac{\partial}{\partial\bm{p}^{2}}\frac{\delta^{2}\partial_{t}\Gamma_{k}}{\delta\pi_{i}(-p)\delta\pi_{j}(p)}\Bigg{|}_{\begin{subarray}{c}p_{0}=0\\\
\bm{p}=0\end{subarray}}\,.$ (36)

The approximation Equation 1 to the effective action together with Equation 36
are based on two approximations: Firstly, in Equation 1 we have dropped the
field-dependence of $Z_{\phi}$, which would lead to different $Z_{\pi}$ and
$Z_{\sigma}$. In Equation 36 we have identified $Z_{\phi}=Z_{\pi}$, and hence
also $Z_{\sigma}=Z_{\pi}$. This is motivated by the fact that the meson
dynamics are only dominant in the broken regime where the three pions are far
lighter than the single sigma mode, which quickly decouples. Hence, the three
pions drive the dynamics.

Furthermore, in Equation 1 we do not distinguish between spatial and temporal
components of $Z_{\phi}$. For finite temperature and density, the Euclidean
$\mathrm{O}(4)$ rotation symmetry is broken, as the heat bath of density
singles out a rest frame. This entails, that $\eta_{\phi,k}$ splits into
$\eta_{\phi,k}^{\perp}$ and $\eta_{\phi,k}^{\parallel}$, the components
transverse and longitudinal to the heat bath/density. We have used the
approximation $\eta_{\phi,k}=\eta_{\phi,k}^{\perp}$ as we have three spatial
directions. The influence of the splitting of $\eta_{\phi,k}$ on the
thermodynamics and baryon number fluctuations has been investigated in detail
e.g. in Yin  _et al._ (2019). There it has been found that the impact is
small, supporting the reliability of the present approximation.

Similarly, the quark anomalous dimension is obtained by projecting the
relevant flow onto the vector channel of the 1PI quark–anti-quark correlation
function,

$\displaystyle\eta_{q}=$ $\displaystyle\frac{1}{4Z_{q,k}}$
$\displaystyle\times\mathrm{Re}\left[\frac{\partial}{\partial\bm{p}^{2}}\mathrm{tr}\left(i\bm{\gamma}\cdot\bm{p}\left(-\frac{\delta^{2}}{\delta\bar{q}(p)\delta
q(p)}\partial_{t}\Gamma_{k}\right)\right)\right]\Bigg{|}_{\begin{subarray}{c}p_{0,ex}\\\
\bm{p}=0\end{subarray}}\,.$ (37)

In Equation 37, the spatial momentum is set to zero, ${\bm{p}}=0$ as in the
mesonic case: vanishing momentum is most relevant to the flow of effective
potential in Equation 32. Note, that the lowest fermionic Matsubara frequency
is non-vanishing. We use $p_{0,ex}\neq 0$, its value is further described in
Appendix B, based on Fu and Pawlowski (2015); Fu  _et al._ (2016).

As is implicit in Equation 37, the flow of the quark two-point function is
complex-valued at non-vanishing chemical potential. This originates in the
Silver-Blaze property of QCD at $T=0$. For quark correlation functions this
entails that they are functions of $p_{0}-i\mu_{q}$ already before the onset
of the baryon density, for a discussion in the present fRG-approach see Khan
_et al._ (2015); Fu and Pawlowski (2015); Fu  _et al._ (2016). In turn, the
couplings are still real (i.e. real functions of the complex variable
$p_{0}-i\mu_{q}$) in particular below the density onset. Hence, couplings
(i.e. expansion coefficients in a Taylor expansion in momenta) are real. This
is readily seen in a resummation of the external frequency of the quark
propagator Fu  _et al._ (2016). Without resummation they are obtained from a
projection on the real part of the flow, see Equation 37.

This projection is also used for the Yukawa coupling. Within the present
approximation, the flow equation of the (real) Yukawa coupling is given by,

$\displaystyle\partial_{t}h_{k}$
$\displaystyle=\frac{1}{2\sigma}\mathrm{Re}\left[\mathrm{tr}\left(-\frac{\delta^{2}}{\delta\bar{q}(p)\delta
q(p)}\partial_{t}\Gamma_{k}\right)\right]\Bigg{|}_{\begin{subarray}{c}p_{0,ex}\\\
\bm{p}=0\end{subarray}}\,.$ (38)

The explicit expressions for the meson and quark anomalous dimensions, as well
as the flow of the Yukawa coupling can be found in Appendix B.

##  Appendix B Flow equations for $V_{k}(\rho)$, $h_{k}$, and $\eta_{\phi,q}$

The flow equation for the effective potential is given in Equation 32. To
resolve its field dependence, we use a Taylor expansion about a $k$-dependent
$\rho$-value $\kappa_{k}$,

$\displaystyle V_{\mathrm{mat},k}(\rho)$
$\displaystyle=\sum_{n=0}^{N_{v}}\frac{\lambda_{n,k}}{n!}(\rho-\kappa_{k})^{n}\,,$
(39)

with the running expansion coefficients $\lambda_{n,k}$. Here, $N_{v}$ is the
maximal order of Taylor expansion included in the numerics. $N_{v}=5$ is
adopted in this work, which is large enough to guarantee the convergence of
expansion, for more details, see e.g., Pawlowski and Rennecke (2014); Yin  _et
al._ (2019). It is more convenient to rewrite Equation 39 by means of the
renormalised variables, i.e.,

$\displaystyle\bar{V}_{\mathrm{mat},k}(\bar{\rho})$
$\displaystyle=\sum_{n=0}^{N_{v}}\frac{\bar{\lambda}_{n,k}}{n!}(\bar{\rho}-\bar{\kappa}_{k})^{n}\,,$
(40)

with $\bar{V}_{\mathrm{mat},k}(\bar{\rho})=V_{\mathrm{mat},k}(\rho)$,
$\bar{\rho}=Z_{\phi,k}\rho$, $\bar{\kappa}_{k}=Z_{\phi,k}\kappa_{k}$, and
$\bar{\lambda}_{n,k}=\lambda_{n,k}/(Z_{\phi,k})^{n}$. Inserting Equation 40
into the l.h.s. of Equation 32 leads us to,

$\displaystyle\partial^{n}_{\bar{\rho}}\left(\partial_{t}\big{|}_{\rho}\bar{V}_{\mathrm{mat},k}(\bar{\rho})\right)\Big{|}_{\bar{\rho}=\bar{\kappa}_{k}}$
$\displaystyle=$
$\displaystyle(\partial_{t}-n\eta_{\phi,k})\bar{\lambda}_{n,k}-(\partial_{t}\bar{\kappa}_{k}+\eta_{\phi,k}\bar{\kappa}_{k})\bar{\lambda}_{n+1,k}\,.$
(41)

In the present work, we use the EoM of $\rho$ as our expansion point. With
Equation 1 this yields,

$\displaystyle\frac{\partial}{\partial\bar{\rho}}\Big{(}\bar{V}_{\mathrm{mat},k}(\bar{\rho})-\bar{c}_{k}\bar{\sigma}\Big{)}\bigg{|}_{\bar{\rho}=\bar{\kappa}_{k}}$
$\displaystyle=0\,,$ (42)

with $\bar{\sigma}=Z_{\phi,k}^{1/2}\sigma$ and
$\bar{c}_{k}=Z_{\phi,k}^{-1/2}c$, with a cutoff-independent $c$. Another
commonly used expansion point is a fixed expansion point,
$\partial_{t}\kappa_{k}=0$. For further details on these two different
expansion approaches, and their respective convergence properties see
Pawlowski and Rennecke (2014); Braun  _et al._ (2016); Rennecke (2015b); Fu
and Pawlowski (2015); Rennecke and Schaefer (2017); Yin  _et al._ (2019).

From Equation 41 and Equation 42 we get the flow equation for the expansion
point,

$\displaystyle\partial_{t}\bar{\kappa}_{k}=$
$\displaystyle\,-\frac{\bar{c}_{k}^{2}}{\bar{\lambda}_{1,k}^{3}+\bar{c}_{k}^{2}\bar{\lambda}_{2,k}}\Bigg{[}\partial_{\bar{\rho}}\left(\partial_{t}\big{|}_{\rho}\bar{V}_{\mathrm{mat},k}(\bar{\rho})\right)\Big{|}_{\bar{\rho}=\bar{\kappa}_{k}}$
$\displaystyle\hskip
42.67912pt+\eta_{\phi,k}\left(\frac{\bar{\lambda}_{1,k}}{2}+\bar{\kappa}_{k}\bar{\lambda}_{2,k}\right)\Bigg{]}\,.$
(43)

The meson anomalous dimension in Equation 36 reads,

$\displaystyle\eta_{\phi,k}=$
$\displaystyle\frac{1}{6\pi^{2}}\Bigg{\\{}\frac{4}{k^{2}}\bar{\kappa}_{k}(\bar{V}^{\prime\prime}_{k}(\bar{\kappa}_{k}))^{2}\mathcal{BB}_{(2,2)}(\tilde{m}^{2}_{\pi,k},\tilde{m}^{2}_{\sigma,k};T)$
$\displaystyle\hskip
17.07182pt+N_{c}\bar{h}^{2}_{k}\bigg{[}\mathcal{F}_{(2)}(\tilde{m}^{2}_{q,k};T,\mu)(2\eta_{q,k}-3)$
$\displaystyle\hskip
17.07182pt-4(\eta_{q,k}-2)\mathcal{F}_{(3)}(\tilde{m}^{2}_{q,k};T,\mu)\bigg{]}\Bigg{\\}}\,,$
(44)

The quark anomalous dimension in Equation 37 reads,

$\displaystyle\eta_{q,k}=$
$\displaystyle\frac{1}{24\pi^{2}N_{f}}(4-\eta_{\phi,k})\bar{h}^{2}_{k}$
$\displaystyle\times\bigg{\\{}(N^{2}_{f}-1)\mathcal{FB}_{(1,2)}(\tilde{m}^{2}_{q,k},\tilde{m}^{2}_{\pi,k};T,\mu,p_{0,ex})$
$\displaystyle+\mathcal{FB}_{(1,2)}(\tilde{m}^{2}_{q,k},\tilde{m}^{2}_{\sigma,k};T,\mu,p_{0,ex})\bigg{\\}}\,.$
(45)

In the threshold function $\mathcal{FB}$’s we have employed $p_{0,ex}=\pi T$
for the finite temperature sector and $p_{0,ex}=\pi T\exp\\{-k/(\pi T)\\}$ for
the vacuum sector. This choice guarantees a consistent temperature dependence
for all $k$, which is particularly relevant for the thermodynamics in the low
temperature region Fu and Pawlowski (2015). This can be resolved by means of a
full frequency summation of the quark external leg Fu  _et al._ (2016), and
the present procedure mimics this physical behaviour.

The flow of the Yukawa coupling in Equation 38 reads,

$\displaystyle\partial_{t}\bar{h}_{k}=$
$\displaystyle\left(\frac{1}{2}\eta_{\phi,k}+\eta_{q,k}\right)\bar{h}_{k}(\bar{\rho})$
$\displaystyle+\frac{\bar{h}^{3}_{k}}{4\pi^{2}N_{f}}\bigg{[}L^{(4)}_{(1,1)}(\tilde{m}^{2}_{q,k},\tilde{m}^{2}_{\sigma,k},\eta_{q,k},\eta_{\phi,k};T,\mu,p_{0,ex})$
$\displaystyle-(N^{2}_{f}-1)L^{(4)}_{(1,1)}(\tilde{m}^{2}_{q,k},\tilde{m}^{2}_{\pi,k},\eta_{q,k},\eta_{\phi,k};T,\mu,p_{0,ex})\bigg{]}\,.$
(46)

Explicit expressions of all the threshold functions mentioned above, such as
$\mathcal{BB}$, $\mathcal{F}$’s, $\mathcal{FB}$’s, and $L$ can be found in
e.g., Fu  _et al._ (2020a); Yin  _et al._ (2019).

In summary, the flow equations Equation 32, Equation 41, Equation 43, Equation
46, supplemented with Equation 44 and Equation 45, constitute a closed set of
ordinary differential equations, which is evolved from the UV cutoff
$k=\Lambda$ to the IR limit $k=0$.

##  Appendix C Initial conditions

To solve the flow equation, we need to specify initial conditions. To this
end, we choose initial values at a scale $k=\Lambda$ such that known
observables of QCD in the vacuum at $k=0$, such as the pion mass and decay
constant, are reproduced. The effective potential at the UV cutoff reads,

$\displaystyle
V_{\mathrm{mat},\Lambda}(\rho)=\frac{\lambda_{\Lambda}}{2}\rho^{2}+\nu_{\Lambda}\rho\,.$
(47)

We initialise the flows at $\Lambda=700$ MeV. The input parameters of the QCD-
assisted LEFT are given in Table 1.

Observables | Value | Parameter in $\Gamma_{\Lambda}$, $\Lambda=700\,$MeV
---|---|---
|  $m_{\pi,\textrm{pol}}$ [MeV]
---
$m_{\sigma}$ [MeV]
$f_{\pi}$ [MeV]
| 136
---
479
92
| $\lambda_{\Lambda}$ | = | 11
---|---|---
$\nu_{k=\Lambda}$ | = | $(0.830\,\mathrm{GeV})^{2}$
$c_{\sigma}$ | = | 2.82$\times 10^{-3}\,\text{GeV}^{3}$
|  |
|  $\bar{m}_{l}$ [MeV]
---
300 |  | $h_{\Lambda}$ | = | 10.18
---|---|---
Table 1: Observables and related initial values for the LEFT couplings at the
initial cutoff scale $\Lambda=700$ MeV. The parameters are fixed with the pion
pole mass $m_{\pi,\textrm{pol}}$, the mass of the sigma resonance,
$m_{\sigma}$, and the pion decay constant,
$f_{\pi}\approx\bar{\sigma}_{\textrm{EoM}}$, the expectation value of the
sigma-field. The input parameters are those of the initial effective potential
$V_{\Lambda}$: the meson self-coupling $\lambda_{\Lambda}$ and the meson mass
parameter $\nu_{\Lambda}$. The pion mass is tuned by $c_{\sigma}$, the
parameter of explicit chiral symmetry breaking. Finally, the constituent quark
mass is fixed via the initial value for the Yukawa coupling, $h_{\Lambda}$.

##  Appendix D Glue potential

The dynamics of the glue sector in QCD is partly imprinted in the glue
potential $V_{\mathrm{glue},k}(A_{0})$, see Equation 2. This has been
discussed in Section II. This potential is only needed for the determination
of the expectation value of the Polyakov loop. Its inherent glue correlation
functions are gapped and their backcoupling is suppressed for $k\lesssim 1$
GeV. Accordingly, we can simply drop the scale dependence of the glue
potential for the present purposes. This leads us to,

$\displaystyle V_{\mathrm{glue}}(L,\bar{L})$
$\displaystyle=V_{\mathrm{glue},k=0}(A_{0})=T^{4}\bar{V}_{\mathrm{glue}}(L,\bar{L})\,.$
(48)

In Equation 48 we have introduced a dimensionless glue potential
$\bar{V}_{\mathrm{glue}}$. Its dependence on the temporal gluon background
field, $A_{0}$, is encoded in the traced Polyakov loop $L[A_{0}]$ and its
conjugate $\bar{L}[A_{0}]$,

$\displaystyle L(\bm{x})$
$\displaystyle=\frac{1}{N_{c}}\left\langle\operatorname{Tr}\,{\cal
P}(\bm{x})\right\rangle\,,\quad\bar{L}(\bm{x})=\frac{1}{N_{c}}\langle\operatorname{Tr}\,{\cal
P}^{\dagger}(\bm{x})\rangle\,,$ (49)

with

$\displaystyle{\cal P}(\bm{x})$
$\displaystyle=\mathcal{P}\exp\Big{(}ig\int_{0}^{\beta}d\tau\hat{A}_{0}(\bm{x},\tau)\Big{)}\,,$
(50)

where $\mathcal{P}$ on the r.h.s. stands for path ordering.

| 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
$a_{i}$ | -44.14 | 151.4 | -90.0677 | 2.77173 | 3.56403
$b_{i}$ | -0.32665 | -82.9823 | 3.0 | 5.85559 |
$c_{i}$ | -50.7961 | 114.038 | -89.4596 | 3.08718 | 6.72812
$d_{i}$ | 27.0885 | -56.0859 | 71.2225 | 2.9715 | 6.61433
Table 2: Values of the parameters for the glue potential in Equation 53 and
Equation 54.

In this work we adopt the parametrisation of the glue potential in Lo  _et
al._ (2013), which reads

$\displaystyle V_{\text{glue}}(L,\bar{L})=$
$\displaystyle-\frac{a(T)}{2}\bar{L}L+b(T)\ln M_{H}(L,\bar{L})$
$\displaystyle+\frac{c(T)}{2}(L^{3}+\bar{L}^{3})+d(T)(\bar{L}L)^{2}\,,$ (51)

with the $\mathrm{SU}(N_{c})$ Haar measure

$\displaystyle M_{H}(L,\bar{L})$
$\displaystyle=1-6\bar{L}L+4(L^{3}+\bar{L}^{3})-3(\bar{L}L)^{2}\,.$ (52)

Both, the parametrisation of glue potential in Equation 51, as well as
determination of relevant parameters in Table 2, is based on lattice results
of $\mathrm{SU}(3)$ Yang-Mills theory at finite temperature. This potential
does not only reproduce the lattice expectation value of the Polyakov loop and
the pressure, but also the correct quadratic fluctuations of the Polyakov
loop, Lo  _et al._ (2013). These fluctuations, and higher ones, are important
for the fluctuation observables discussed here Fu and Pawlowski (2015); Fu
_et al._ (2016). The coefficients in Equation 51 are temperature-dependent,

$\displaystyle x(T)$
$\displaystyle=\frac{x_{1}+x_{2}/(t_{r}+1)+x_{3}/(t_{r}+1)^{2}}{1+x_{4}/(t_{r}+1)+x_{5}/(t_{r}+1)^{2}}\,,$
(53)

with $x=a,c,d$, and

$\displaystyle b(T)$
$\displaystyle=b_{1}(t_{r}+1)^{-b_{4}}\left(1-e^{b_{2}/(t_{r}+1)^{b_{3}}}\right)\,.$
(54)

In Equation 53 and Equation 54 we have used the reduced temperature
$t_{r}=(T-T_{c})/T_{c}$. The parameter values are taken from Lo  _et al._
(2013), and are collected in Table 2 for convenience.

The parameters in Table 2 are that of the glue potential in Yang-Mills theory.
It has been argued and shown in Pawlowski (2011); Haas  _et al._ (2013);
Herbst  _et al._ (2014) that unquenching effects in QCD are well captured by a
linear rescaling of the reduced temperature in the regime about $T_{c}$, very
similar to the rescaling discussed in Section II.2.2. This leads us to,

$\displaystyle(t_{r})_{\text{\tiny{YM}}}$
$\displaystyle\rightarrow\alpha\,(t_{r})_{\text{\tiny{glue}}}\,,$ (55)

with

$\displaystyle(t_{r})_{\text{\tiny{glue}}}$
$\displaystyle=(T-T_{c}^{\text{\tiny{glue}}})/T_{c}^{\text{\tiny{glue}}}\,.$
(56)

In the present work we have used $\alpha=0.75$ and
$T_{c}^{\text{\tiny{glue}}}=213\,{\rm MeV}$.

## References

  * Stephanov (2006a) M. A. Stephanov, “QCD phase diagram: An Overview,”  _Proceedings, 24th International Symposium on Lattice Field Theory (Lattice 2006): Tucson, USA, July 23-28, 2006_ , PoS LAT2006, 024 (2006a), arXiv:hep-lat/0701002 [hep-lat] .
  * Friman  _et al._ (2011a) Bengt Friman, Claudia Hohne, Jorn Knoll, Stefan Leupold, Jorgen Randrup, Ralf Rapp, and Peter Senger, “The CBM physics book: Compressed baryonic matter in laboratory experiments,” Lect. Notes Phys. 814, pp.1–980 (2011a).
  * Blaschke  _et al._ (2016) David Blaschke, Jörg Aichelin, Elena Bratkovskaya, Volker Friese, Marek Gazdzicki, Jørgen Randrup, Oleg Rogachevsky, Oleg Teryaev, and Viacheslav Toneev, “Topical issue on exploring strongly interacting matter at high densities - nica white paper,” The European Physical Journal A 52, 267 (2016).
  * Luo and Xu (2017) Xiaofeng Luo and Nu Xu, “Search for the QCD Critical Point with Fluctuations of Conserved Quantities in Relativistic Heavy-Ion Collisions at RHIC : An Overview,” Nucl. Sci. Tech. 28, 112 (2017), arXiv:1701.02105 [nucl-ex] .
  * Dainese  _et al._ (2019) A. Dainese _et al._ (QCD Working Group), “Physics Beyond Colliders: QCD Working Group Report,” (2019), arXiv:1901.04482 [hep-ex] .
  * Bzdak  _et al._ (2020) Adam Bzdak, Shinichi Esumi, Volker Koch, Jinfeng Liao, Mikhail Stephanov, and Nu Xu, “Mapping the Phases of Quantum Chromodynamics with Beam Energy Scan,” Phys. Rept. 853, 1–87 (2020), arXiv:1906.00936 [nucl-th] .
  * Fischer (2019) Christian S. Fischer, “QCD at finite temperature and chemical potential from Dyson-Schwinger equations,” Prog. Part. Nucl. Phys. 105, 1–60 (2019), arXiv:1810.12938 [hep-ph] .
  * Fu  _et al._ (2020a) Wei-jie Fu, Jan M. Pawlowski, and Fabian Rennecke, “QCD phase structure at finite temperature and density,” Phys. Rev. D 101, 054032 (2020a), arXiv:1909.02991 [hep-ph] .
  * Bazavov  _et al._ (2020) A. Bazavov _et al._ , “Skewness, kurtosis, and the fifth and sixth order cumulants of net baryon-number distributions from lattice QCD confront high-statistics STAR data,” Phys. Rev. D 101, 074502 (2020), arXiv:2001.08530 [hep-lat] .
  * Borsanyi  _et al._ (2020) Szabolcs Borsanyi, Zoltan Fodor, Jana N. Guenther, Ruben Kara, Sandor D. Katz, Paolo Parotto, Attila Pasztor, Claudia Ratti, and Kalman K. Szabo, “QCD Crossover at Finite Chemical Potential from Lattice Simulations,” Phys. Rev. Lett. 125, 052001 (2020), arXiv:2002.02821 [hep-lat] .
  * Lü  _et al._ (2017) LiMing Lü, Han Yi, ZhiGang Xiao, Ming Shao, Song Zhang, GuoQing Xiao, and Nu Xu, “Conceptual design of the HIRFL-CSR external-target experiment,” Sci. China Phys. Mech. Astron. 60, 012021 (2017).
  * Yang  _et al._ (2013) J. C. Yang _et al._ , “High Intensity heavy ion Accelerator Facility (HIAF) in China,”  _Proceedings, 16th International Conference on Electromagnetic Isotope Separators and Techniques Related to their Applications (EMIS 2012): Matsue, Japan, December 2-7, 2012_ , Nucl. Instrum. Meth. B317, 263–265 (2013).
  * Sako (2019) H. Sako (J-PARC-HI), “Studies of extremely dense matter in heavy-ion collisions at J-PARC,” Nucl. Phys. A 982, 959–962 (2019).
  * Stephanov  _et al._ (1999) Misha A. Stephanov, K. Rajagopal, and Edward V. Shuryak, “Event-by-event fluctuations in heavy ion collisions and the QCD critical point,” Phys. Rev. D60, 114028 (1999), arXiv:hep-ph/9903292 [hep-ph] .
  * Stephanov (2009) M.A. Stephanov, “Non-Gaussian fluctuations near the QCD critical point,” Phys. Rev. Lett. 102, 032301 (2009), arXiv:0809.3450 [hep-ph] .
  * Stephanov (2011) M.A. Stephanov, “On the sign of kurtosis near the QCD critical point,” Phys. Rev. Lett. 107, 052301 (2011), arXiv:1104.1627 [hep-ph] .
  * Adamczyk  _et al._ (2014a) L. Adamczyk _et al._ (STAR), “Energy Dependence of Moments of Net-proton Multiplicity Distributions at RHIC,” Phys. Rev. Lett. 112, 032302 (2014a), arXiv:1309.5681 [nucl-ex] .
  * Adamczyk  _et al._ (2014b) L. Adamczyk _et al._ (STAR), “Beam energy dependence of moments of the net-charge multiplicity distributions in Au+Au collisions at RHIC,” Phys. Rev. Lett. 113, 092301 (2014b), arXiv:1402.1558 [nucl-ex] .
  * Luo (2015) Xiaofeng Luo (STAR), “Energy Dependence of Moments of Net-Proton and Net-Charge Multiplicity Distributions at STAR,”  _Proceedings, 9th International Workshop on Critical Point and Onset of Deconfinement (CPOD 2014): Bielefeld, Germany, November 17-21, 2014_ , PoS CPOD2014, 019 (2015), arXiv:1503.02558 [nucl-ex] .
  * Adamczyk  _et al._ (2018) L. Adamczyk _et al._ (STAR), “Collision Energy Dependence of Moments of Net-Kaon Multiplicity Distributions at RHIC,” Phys. Lett. B785, 551–560 (2018), arXiv:1709.00773 [nucl-ex] .
  * Adam  _et al._ (2019) Jaroslav Adam _et al._ (STAR), “Collision-energy dependence of second-order off-diagonal and diagonal cumulants of net-charge, net-proton, and net-kaon multiplicity distributions in Au + Au collisions,” Phys. Rev. C 100, 014902 (2019), arXiv:1903.05370 [nucl-ex] .
  * Adam  _et al._ (2020) J. Adam _et al._ (STAR), “Net-proton number fluctuations and the Quantum Chromodynamics critical point,” (2020), arXiv:2001.02852 [nucl-ex] .
  * Nonaka (2020) Toshihiro Nonaka (STAR), “Measurement of the Sixth-Order Cumulant of Net-Proton Distributions in Au+Au Collisions from the STAR Experiment,” in  _28th International Conference on Ultrarelativistic Nucleus-Nucleus Collisions_ (2020) arXiv:2002.12505 [nucl-ex] .
  * Pandav (2020) Ashish Pandav (STAR), “Measurement of cumulants of conserved charge multiplicity distributions in Au+Au collisions from the STAR experiment,” (2020), arXiv:2003.12503 [nucl-ex] .
  * Braun  _et al._ (2020a) Jens Braun, Marc Leonhardt, and Martin Pospiech, “Fierz-complete NJL model study III: Emergence from quark-gluon dynamics,” Phys. Rev. D101, 036004 (2020a), arXiv:1909.06298 [hep-ph] .
  * Isserstedt  _et al._ (2019) Philipp Isserstedt, Michael Buballa, Christian S. Fischer, and Pascal J. Gunkel, “Baryon number fluctuations in the QCD phase diagram from Dyson-Schwinger equations,” Phys. Rev. D 100, 074011 (2019), arXiv:1906.11644 [hep-ph] .
  * Gao and Pawlowski (2020a) Fei Gao and Jan M. Pawlowski, “QCD phase structure from functional methods,” Phys. Rev. D 102, 034027 (2020a), arXiv:2002.07500 [hep-ph] .
  * Gao and Pawlowski (2020b) Fei Gao and Jan M. Pawlowski, “Chiral phase structure and critical end point in QCD,” (2020b), arXiv:2010.13705 [hep-ph] .
  * Schaefer and Wambach (2007) Bernd-Jochen Schaefer and Jochen Wambach, “Susceptibilities near the QCD (tri)critical point,” Phys. Rev. D75, 085015 (2007), arXiv:hep-ph/0603256 [hep-ph] .
  * Braun-Munzinger  _et al._ (2004a) Peter Braun-Munzinger, Krzysztof Redlich, and Johanna Stachel, “Particle production in heavy ion collisions,” in  _Quark-Gluon Plasma_, edited by Rudolph C. Hwa and Xin-Nian Wang (2004) pp. 491–599, arXiv:nucl-th/0304013 .
  * Braun-Munzinger  _et al._ (2004b) P. Braun-Munzinger, J. Stachel, and Christof Wetterich, “Chemical freezeout and the QCD phase transition temperature,” Phys. Lett. B 596, 61–69 (2004b), arXiv:nucl-th/0311005 .
  * Floerchinger and Wetterich (2012) Stefan Floerchinger and Christof Wetterich, “Chemical freeze-out in heavy ion collisions at large baryon densities,” Nucl. Phys. A890-891, 11–24 (2012), arXiv:1202.1671 [nucl-th] .
  * Bazavov  _et al._ (2012) A. Bazavov _et al._ , “Freeze-out Conditions in Heavy Ion Collisions from QCD Thermodynamics,” Phys. Rev. Lett. 109, 192302 (2012), arXiv:1208.1220 [hep-lat] .
  * Borsanyi  _et al._ (2013) S. Borsanyi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo, “Freeze-out parameters: lattice meets experiment,” Phys. Rev. Lett. 111, 062005 (2013), arXiv:1305.5161 [hep-lat] .
  * Borsanyi  _et al._ (2014) S. Borsanyi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo, “Freeze-out parameters from electric charge and baryon number fluctuations: is there consistency?” Phys. Rev. Lett. 113, 052301 (2014), arXiv:1403.4576 [hep-lat] .
  * Bazavov  _et al._ (2017a) A. Bazavov _et al._ , “The QCD Equation of State to $\mathcal{O}(\mu_{B}^{6})$ from Lattice QCD,” Phys. Rev. D 95, 054504 (2017a), arXiv:1701.04325 [hep-lat] .
  * Bazavov  _et al._ (2017b) A. Bazavov _et al._ (HotQCD), “Skewness and kurtosis of net baryon-number distributions at small values of the baryon chemical potential,” Phys. Rev. D96, 074510 (2017b), arXiv:1708.04897 [hep-lat] .
  * Borsanyi  _et al._ (2018) Szabolcs Borsanyi, Zoltan Fodor, Jana N. Guenther, Sandor K. Katz, Kalman K. Szabo, Attila Pasztor, Israel Portillo, and Claudia Ratti, “Higher order fluctuations and correlations of conserved charges from lattice QCD,” JHEP 10, 205 (2018), arXiv:1805.04445 [hep-lat] .
  * Parotto  _et al._ (2020) Paolo Parotto, Marcus Bluhm, Debora Mroczek, Marlene Nahrgang, Jacquelyn Noronha-Hostler, Krishna Rajagopal, Claudia Ratti, Thomas Schäfer, and Mikhail Stephanov, “QCD equation of state matched to lattice data and exhibiting a critical point singularity,” Phys. Rev. C 101, 034901 (2020), arXiv:1805.05249 [hep-ph] .
  * Fu and Pawlowski (2015) Wei-jie Fu and Jan M. Pawlowski, “Relevance of matter and glue dynamics for baryon number fluctuations,” Phys. Rev. D92, 116006 (2015), arXiv:1508.06504 [hep-ph] .
  * Fu and Pawlowski (2016) Wei-jie Fu and Jan M. Pawlowski, “Correlating the skewness and kurtosis of baryon number distributions,” Phys. Rev. D93, 091501 (2016), arXiv:1512.08461 [hep-ph] .
  * Fu  _et al._ (2016) Wei-jie Fu, Jan M. Pawlowski, Fabian Rennecke, and Bernd-Jochen Schaefer, “Baryon number fluctuations at finite temperature and density,” Phys. Rev. D 94, 116020 (2016), arXiv:1608.04302 [hep-ph] .
  * Mitter  _et al._ (2015) Mario Mitter, Jan M. Pawlowski, and Nils Strodthoff, “Chiral symmetry breaking in continuum QCD,” Phys. Rev. D91, 054035 (2015), arXiv:1411.7978 [hep-ph] .
  * Braun  _et al._ (2016) Jens Braun, Leonard Fister, Jan M. Pawlowski, and Fabian Rennecke, “From Quarks and Gluons to Hadrons: Chiral Symmetry Breaking in Dynamical QCD,” Phys. Rev. D94, 034016 (2016), arXiv:1412.1045 [hep-ph] .
  * Rennecke (2015a) Fabian Rennecke, “Vacuum structure of vector mesons in QCD,” Phys. Rev. D92, 076012 (2015a), arXiv:1504.03585 [hep-ph] .
  * Cyrol  _et al._ (2018a) Anton K. Cyrol, Mario Mitter, Jan M. Pawlowski, and Nils Strodthoff, “Nonperturbative quark, gluon, and meson correlators of unquenched QCD,” Phys. Rev. D97, 054006 (2018a), arXiv:1706.06326 [hep-ph] .
  * Gies and Wetterich (2002) Holger Gies and Christof Wetterich, “Renormalization flow of bound states,” Phys. Rev. D65, 065001 (2002), arXiv:hep-th/0107221 [hep-th] .
  * Gies and Wetterich (2004) Holger Gies and Christof Wetterich, “Universality of spontaneous chiral symmetry breaking in gauge theories,” Phys. Rev. D69, 025001 (2004), arXiv:hep-th/0209183 [hep-th] .
  * Pawlowski (2007) Jan M. Pawlowski, “Aspects of the functional renormalisation group,” Annals Phys. 322, 2831–2915 (2007), arXiv:hep-th/0512261 [hep-th] .
  * Floerchinger and Wetterich (2009) S. Floerchinger and C. Wetterich, “Exact flow equation for composite operators,” Phys. Lett. B680, 371–376 (2009), arXiv:0905.0915 [hep-th] .
  * Dupuis  _et al._ (2020) N. Dupuis, L. Canet, A. Eichhorn, W. Metzner, J.M. Pawlowski, M. Tissier, and N. Wschebor, “The nonperturbative functional renormalization group and its applications,” (2020), arXiv:2006.04853 [cond-mat.stat-mech] .
  * Haas  _et al._ (2013) Lisa M. Haas, Rainer Stiele, Jens Braun, Jan M. Pawlowski, and Jürgen Schaffner-Bielich, “Improved Polyakov-loop potential for effective models from functional calculations,” Phys. Rev. D87, 076004 (2013), arXiv:1302.1993 [hep-ph] .
  * Herbst  _et al._ (2014) Tina Katharina Herbst, Mario Mitter, Jan M. Pawlowski, Bernd-Jochen Schaefer, and Rainer Stiele, “Thermodynamics of QCD at vanishing density,” Phys. Lett. B731, 248–256 (2014), arXiv:1308.3621 [hep-ph] .
  * Wen  _et al._ (2019) Rui Wen, Chuang Huang, and Wei-Jie Fu, “Baryon number fluctuations in the 2+1 flavor low energy effective model,” Phys. Rev. D 99, 094019 (2019), arXiv:1809.04233 [hep-ph] .
  * Yin  _et al._ (2019) Shi Yin, Rui Wen, and Wei-jie Fu, “Mesonic dynamics and the QCD phase transition,” Phys. Rev. D 100, 094029 (2019), arXiv:1907.10262 [hep-ph] .
  * Fu  _et al._ (2020b) Wei-jie Fu, Jan M. Pawlowski, and Fabian Rennecke, “Strangeness Neutrality and QCD Thermodynamics,” SciPost Phys. Core 2, 002 (2020b), arXiv:1808.00410 [hep-ph] .
  * Fu  _et al._ (2019) Wei-jie Fu, Jan M. Pawlowski, and Fabian Rennecke, “Strangeness neutrality and baryon-strangeness correlations,” Phys. Rev. D 100, 111501 (2019), arXiv:1809.01594 [hep-ph] .
  * Skokov  _et al._ (2010) V. Skokov, B. Stokic, B. Friman, and K. Redlich, “Meson fluctuations and thermodynamics of the Polyakov loop extended quark-meson model,” Phys. Rev. C82, 015206 (2010), arXiv:1004.2665 [hep-ph] .
  * Skokov  _et al._ (2011) V. Skokov, B. Friman, and K. Redlich, “Quark number fluctuations in the Polyakov loop-extended quark-meson model at finite baryon density,” Phys. Rev. C83, 054904 (2011), arXiv:1008.4570 [hep-ph] .
  * Friman  _et al._ (2011b) B. Friman, F. Karsch, K. Redlich, and V. Skokov, “Fluctuations as probe of the QCD phase transition and freeze-out in heavy ion collisions at LHC and RHIC,” Eur. Phys. J. C 71, 1694 (2011b), arXiv:1103.3511 [hep-ph] .
  * Morita  _et al._ (2015) Kenji Morita, Bengt Friman, and Krzysztof Redlich, “Criticality of the net-baryon number probability distribution at finite density,” Phys. Lett. B741, 178–183 (2015), arXiv:1402.5982 [hep-ph] .
  * Almasi  _et al._ (2017a) Gabor Andras Almasi, Bengt Friman, and Krzysztof Redlich, “Baryon number fluctuations in chiral effective models and their phenomenological implications,” Phys. Rev. D 96, 014027 (2017a), arXiv:1703.05947 [hep-ph] .
  * Xin  _et al._ (2014) Xian-yin Xin, Si-xue Qin, and Yu-xin Liu, “Quark number fluctuations at finite temperature and finite chemical potential via the Dyson-Schwinger equation approach,” Phys. Rev. D 90, 076006 (2014).
  * Fu  _et al._ (2010) Wei-jie Fu, Yu-xin Liu, and Yue-Liang Wu, “Fluctuations and Correlations of Conserved Charges in QCD at Finite Temperature with Effective Models,” Phys. Rev. D 81, 014028 (2010), arXiv:0910.5783 [hep-ph] .
  * Fu and Wu (2010) Wei-jie Fu and Yue-liang Wu, “Fluctuations and Correlations of Conserved Charges near the QCD Critical Point,” Phys. Rev. D 82, 074013 (2010), arXiv:1008.3684 [hep-ph] .
  * Karsch  _et al._ (2011) Frithjof Karsch, Bernd-Jochen Schaefer, Mathias Wagner, and Jochen Wambach, “Towards finite density QCD with Taylor expansions,” Phys. Lett. B 698, 256–264 (2011), arXiv:1009.5211 [hep-ph] .
  * Schaefer and Wagner (2012) B.J. Schaefer and M. Wagner, “QCD critical region and higher moments for three flavor models,” Phys. Rev. D 85, 034027 (2012), arXiv:1111.6871 [hep-ph] .
  * Li  _et al._ (2019a) Zhibin Li, Kun Xu, Xinyang Wang, and Mei Huang, “The kurtosis of net baryon number fluctuations from a realistic Polyakov–Nambu–Jona-Lasinio model along the experimental freeze-out line,” Eur. Phys. J. C 79, 245 (2019a), arXiv:1801.09215 [hep-ph] .
  * Schaefer  _et al._ (2007) Bernd-Jochen Schaefer, Jan M. Pawlowski, and Jochen Wambach, “The Phase Structure of the Polyakov–Quark-Meson Model,” Phys. Rev. D76, 074023 (2007), arXiv:0704.3234 [hep-ph] .
  * Schaefer  _et al._ (2010) Bernd-Jochen Schaefer, Mathias Wagner, and Jochen Wambach, “Thermodynamics of (2+1)-flavor QCD: Confronting Models with Lattice Studies,” Phys. Rev. D 81, 074013 (2010), arXiv:0910.5628 [hep-ph] .
  * Herbst  _et al._ (2011) Tina Katharina Herbst, Jan M. Pawlowski, and Bernd-Jochen Schaefer, “The phase structure of the Polyakovquarkmeson model beyond mean field,” Phys. Lett. B696, 58–67 (2011), arXiv:1008.0081 [hep-ph] .
  * Morita  _et al._ (2011) Kenji Morita, Vladimir Skokov, Bengt Friman, and Krzysztof Redlich, “Role of mesonic fluctuations in the Polyakov loop extended quark-meson model at imaginary chemical potential,” Phys. Rev. D 84, 074020 (2011), arXiv:1108.0735 [hep-ph] .
  * Skokov  _et al._ (2012) V. Skokov, B. Friman, and K. Redlich, “Non-perturbative dynamics and charge fluctuations in effective chiral models,” Phys. Lett. B 708, 179–185 (2012), arXiv:1108.3231 [hep-ph] .
  * Mintz  _et al._ (2013) Bruno W. Mintz, Rainer Stiele, Rudnei O. Ramos, and Juergen Schaffner-Bielich, “Phase diagram and surface tension in the three-flavor Polyakov-quark-meson model,” Phys. Rev. D 87, 036004 (2013), arXiv:1212.1184 [hep-ph] .
  * Herbst  _et al._ (2013) Tina K. Herbst, Jan M. Pawlowski, and Bernd-Jochen Schaefer, “Phase structure and thermodynamics of QCD,” Phys. Rev. D88, 014007 (2013), arXiv:1302.1426 [hep-ph] .
  * Stiele and Schaffner-Bielich (2016) Rainer Stiele and Juergen Schaffner-Bielich, “Phase diagram and nucleation in the Polyakov-loop-extended Quark-Meson truncation of QCD with the unquenched Polyakov-loop potential,” Phys. Rev. D 93, 094014 (2016), arXiv:1601.05731 [hep-ph] .
  * Sun  _et al._ (2018) Ke-xin Sun, Rui Wen, and Wei-jie Fu, “Baryon number probability distribution at finite temperature,” Phys. Rev. D98, 074028 (2018), arXiv:1805.12025 [hep-ph] .
  * Hansen  _et al._ (2020) Hubert Hansen, Rainer Stiele, and Pedro Costa, “Quark and Polyakov-loop correlations in effective models at zero and nonvanishing density,” Phys. Rev. D 101, 094001 (2020), arXiv:1904.08965 [hep-ph] .
  * Bazavov  _et al._ (2019) A. Bazavov _et al._ (HotQCD), “Chiral crossover in QCD at zero and non-zero chemical potentials,” Phys. Lett. B795, 15–21 (2019), arXiv:1812.08235 [hep-lat] .
  * Bellwied  _et al._ (2015) R. Bellwied, S. Borsanyi, Z. Fodor, J. Guenther, S. D. Katz, C. Ratti, and K. K. Szabo, “The QCD phase diagram from analytic continuation,” Phys. Lett. B751, 559–564 (2015), arXiv:1507.07510 [hep-lat] .
  * Buballa and Carignano (2015) Michael Buballa and Stefano Carignano, “Inhomogeneous chiral condensates,” Prog. Part. Nucl. Phys. 81, 39–96 (2015), arXiv:1406.1367 [hep-ph] .
  * Wagner  _et al._ (2010) Mathias Wagner, Andrea Walther, and Bernd-Jochen Schaefer, “On the efficient computation of high-order derivatives for implicitly defined functions,” Comput. Phys. Commun. 181, 756–764 (2010), arXiv:0912.2208 [hep-ph] .
  * Roberge and Weiss (1986) Andre Roberge and Nathan Weiss, “Gauge Theories With Imaginary Chemical Potential and the Phases of QCD,” Nucl. Phys. B275, 734–745 (1986).
  * Yang and Lee (1952) Chen-Ning Yang and T.D. Lee, “Statistical theory of equations of state and phase transitions. 1. Theory of condensation,” Phys. Rev. 87, 404–409 (1952).
  * Braun  _et al._ (2011a) Jens Braun, Lisa M. Haas, Florian Marhauser, and Jan M. Pawlowski, “Phase Structure of Two-Flavor QCD at Finite Chemical Potential,” Phys. Rev. Lett. 106, 022002 (2011a), arXiv:0908.0008 [hep-ph] .
  * Philipsen and Pinke (2014) Owe Philipsen and Christopher Pinke, “Nature of the Roberge-Weiss transition in $N_{f}=2$ QCD with Wilson fermions,” Phys. Rev. D 89, 094504 (2014), arXiv:1402.0838 [hep-lat] .
  * Stephanov (2006b) M.A. Stephanov, “QCD critical point and complex chemical potential singularities,” Phys. Rev. D 73, 094508 (2006b), arXiv:hep-lat/0603014 .
  * Mukherjee and Skokov (2019) Swagato Mukherjee and Vladimir Skokov, “Universality driven analytic structure of QCD crossover: radius of convergence in baryon chemical potential,” (2019), arXiv:1909.04639 [hep-ph] .
  * Connelly  _et al._ (2020) Andrew Connelly, Gregory Johnson, Fabian Rennecke, and Vladimir Skokov, “Universal Location of the Yang-Lee Edge Singularity in $O(N)$ Theories,” Phys. Rev. Lett. 125, 191602 (2020), arXiv:2006.12541 [cond-mat.stat-mech] .
  * Andronic  _et al._ (2018) Anton Andronic, Peter Braun-Munzinger, Krzysztof Redlich, and Johanna Stachel, “Decoding the phase structure of QCD via particle production at high energy,” Nature 561, 321–330 (2018), arXiv:1710.09425 [nucl-th] .
  * Adamczyk  _et al._ (2017) L. Adamczyk _et al._ (STAR), “Bulk Properties of the Medium Produced in Relativistic Heavy-Ion Collisions from the Beam Energy Scan Program,” Phys. Rev. C96, 044904 (2017), arXiv:1701.07065 [nucl-ex] .
  * Bzdak  _et al._ (2017) Adam Bzdak, Volker Koch, and Nils Strodthoff, “Cumulants and correlation functions versus the QCD phase diagram,” Phys. Rev. C 95, 054906 (2017), arXiv:1607.07375 [nucl-th] .
  * He and Luo (2017) Shu He and Xiaofeng Luo, “Proton Cumulants and Correlation Functions in Au + Au Collisions at $\sqrt{s_{\mathrm{NN}}}$=7.7-200 GeV from UrQMD Model,” Phys. Lett. B 774, 623–629 (2017), arXiv:1704.00423 [nucl-ex] .
  * Abdallah  _et al._ (2021) Mohamed Abdallah _et al._ (STAR), “Measurement of the sixth-order cumulant of net-proton multiplicity distributions in Au+Au collisions at $\sqrt{s_{\rm NN}}=$ 27, 54.4, and 200 GeV at RHIC,” (2021), arXiv:2105.14698 [nucl-ex] .
  * Luo  _et al._ (2013) Xiaofeng Luo, Ji Xu, Bedangadas Mohanty, and Nu Xu, “Volume fluctuation and auto-correlation effects in the moment analysis of net-proton multiplicity distributions in heavy-ion collisions,” J. Phys. G 40, 105104 (2013), arXiv:1302.2332 [nucl-ex] .
  * Chatterjee  _et al._ (2020a) Arghya Chatterjee, Yu Zhang, Jingdong Zeng, Nihar Ranjan Sahoo, and Xiaofeng Luo, “Effect of centrality selection on higher-order cumulants of net-proton multiplicity distributions in relativistic heavy-ion collisions,” Phys. Rev. C 101, 034902 (2020a), arXiv:1910.08004 [nucl-ex] .
  * Chatterjee  _et al._ (2020b) Arghya Chatterjee, Yu Zhang, Hui Liu, Ruiqin Wang, Shu He, and Xiaofeng Luo, “Effects of centrality fluctuation and deuteron formation on proton number cumulant in Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}$ = 3 GeV from JAM model,” (2020b), arXiv:2009.03755 [nucl-ex] .
  * Braun  _et al._ (2011b) Jens Braun, Bertram Klein, and Piotr Piasecki, “On the scaling behavior of the chiral phase transition in QCD in finite and infinite volume,” Eur. Phys. J. C 71, 1576 (2011b), arXiv:1008.2155 [hep-ph] .
  * Braun  _et al._ (2012) Jens Braun, Bertram Klein, and Bernd-Jochen Schaefer, “On the Phase Structure of QCD in a Finite Volume,” Phys. Lett. B713, 216–223 (2012), arXiv:1110.0849 [hep-ph] .
  * Tripolt  _et al._ (2014) Ralf-Arno Tripolt, Jens Braun, Bertram Klein, and Bernd-Jochen Schaefer, “Effect of fluctuations on the QCD critical point in a finite volume,” Phys. Rev. D90, 054012 (2014), arXiv:1308.0164 [hep-ph] .
  * Almasi  _et al._ (2017b) Gabor Almasi, Robert Pisarski, and Vladimir Skokov, “Volume dependence of baryon number cumulants and their ratios,” Phys. Rev. D95, 056015 (2017b), arXiv:1612.04416 [hep-ph] .
  * Klein (2017) Bertram Klein, “Modeling Finite-Volume Effects and Chiral Symmetry Breaking in Two-Flavor QCD Thermodynamics,” Phys. Rept. 707-708, 1–51 (2017), arXiv:1710.05357 [hep-ph] .
  * Li  _et al._ (2019b) Bo-Lin Li, Zhu-Fang Cui, Bo-Wen Zhou, Sun An, Li-Ping Zhang, and Hong-Shi Zong, “Finite volume effects on the chiral phase transition from DysonSchwinger equations of QCD,” Nucl. Phys. B938, 298–306 (2019b), arXiv:1711.04914 [hep-ph] .
  * Liu  _et al._ (2020) Run-Lin Liu, Meng-Yun Lai, Chao Shi, and Hong-Shi Zong, “Finite volume effects on QCD susceptibilities with a chiral chemical potential,” Phys. Rev. D 102, 014014 (2020).
  * Wan  _et al._ (2020) Shen-Song Wan, Daize Li, Bonan Zhang, and Marco Ruggieri, “Finite Size Effects on the Chiral Phase Transition of Quantum Chromodynamics,” (2020), arXiv:2012.05734 [hep-ph] .
  * He  _et al._ (2016) Shu He, Xiaofeng Luo, Yasushi Nara, ShinIchi Esumi, and Nu Xu, “Effects of Nuclear Potential on the Cumulants of Net-Proton and Net-Baryon Multiplicity Distributions in Au+Au Collisions at $\sqrt{s_{\text{NN}}}=5\,\text{GeV}$,” Phys. Lett. B 762, 296–300 (2016), arXiv:1607.06376 [nucl-ex] .